I can try to help you mainly about point 2 and 3. What I say comes mainly from some lectures I attended. In general you may be interested in Coleman aspects of symmetry for some very similar discussion about these topics. Especially for theory of scale invariance as I haven't told much about it. For point 1 I don't know a formal proof, only some obvious arguments so I'll avoid it.
Concerning your last two questions, I think that to understand the conceptual and practical relevance of the running coupling constant is important to understand how renormalization modifies the theory of scale invariance.
Basically, if you derive the Ward identity (WI) for scale invariance, you obtain and identity which shows that, under a scale transformation, a green function scales just as dimensional analysis tells you:
$$\mu\frac{\partial}{\partial\mu}\Gamma^{(n)}(p_1, ..., p_{n-1})=-(D-4)\Gamma^{(n)}(p_1, ..., p_{n-1}) \tag{1}\label{1}$$
where $\Gamma^{(n)}(p_1, ..., p_{n-1})$ is the 1PI green function obtained from a generic $G^{n}(x_1, ..., x_n)=\langle 0|T[\phi(x_1)...\phi(x_n)]|0\rangle$. There are n-1 moemntum as the delta of conservation of momentum was factorized away, while $D=n\text{dim}(\phi)$. You can check that by this conditions $\Gamma^{(n)}$ has actually dimension of mass $-(D-4)$ (I think Lewis and Ryder or Srednicky books have a table with dimensions of green functions which could be useful to prove it).
It is evident that (1) is a dimensional analysis statement because
$$\mu\frac{\partial}{\partial\mu}\Gamma^{(n)}(p_1, ..., p_{n-1})=\mu\frac{\partial}{\partial\mu}\Gamma^{(n)}(p_1/\mu, ..., p_{n-1}/\mu, \mu)=\mu\frac{\partial}{\partial\mu}\Gamma^{(n)}(\Pi_1, ..., \Pi_{n-1},\mu)=\mu\frac{\partial}{\partial\mu}\tilde{\Gamma}^{(n)}(\Pi_1, ..., \Pi_{n-1})\mu^{-(D-4)}=-(D-4)\Gamma^{(n)}(p_1, ..., p_{n-1}) $$
This WI is correct for bare green functions. Though renormalization, bringing in a new scale $\mu_R$, spoils it. (from now on you can assume I have identified the parameter of the scale transformation $\mu$ with the renormalization scale, which are both arbitrary)
For example consider $\phi^4$ ind 4 dimension, with $\lambda$ the adimensional coupling constant, and with null mass as we suppose we are in the high energy limit. Consider the two in two scattering amplitude $\Gamma^{(4)}(s, \lambda)=\lambda$.
It is obvious that $\Gamma^{(4)}$ at the leading order LO satisfies the WI, as $\mu\frac{\partial}{\partial\mu}\Gamma^{(n)}=\mu\frac{\partial}{\partial\mu}\lambda=0$. But yet it is easy to understand how after renormalization at next to leading order NLO $\mu\frac{\partial}{\partial\mu}\Gamma^{(n)}\neq 0$, as $\lambda$ takes a $\mu$ dependence even if is adimensional. The 1PI green function now can have both an explicit and implicit dependence on $\mu$, which enters this way $\Gamma^{(4)}(s,\lambda(\mu^2), \mu^2)=\Gamma^{(4)}(\frac{s}{\mu^2},\lambda(\mu^2))$. What must be true though is that the renormalization scale is unobservable in a physical quantity, so the $\Gamma^{(4)}$ must satisfy $\mu^2\frac{d}{d\mu^2}\Gamma^{(4)}(\frac{s}{\mu^2},\lambda(\mu^2))=0$. The last equation is actually the renormalization group equation for observables:
$$\frac{\partial\Gamma^{(4)}}{\partial ln(\frac{s}{\mu^2})}=\beta(\lambda(\mu^2))\frac{\partial\Gamma^{(4)}}{\partial\lambda(\mu^2)} \tag{2}\label{2}$$
(in general the RGE eq or Callan-Symanzik are the equations that corrects the (1) for generic green functions, which can have an anomalous dimension).
Now what you can say is that, given the definition of the running coupling constant (RCC) $\frac{d\lambda_R(s)}{dln(\frac{s}{\mu^2})}=\beta(\lambda_R)$ with $\lambda_R(\mu^2)=\lambda(\mu^2)$ as a inital condition, by Dini's theorem the RCC satisfy the RGE equation for physical observables, that is
$$\frac{\partial\lambda_R}{\partial ln(\frac{s}{\mu^2})}=\beta(\lambda(\mu^2))\frac{\partial\lambda_R}{\partial\lambda(\mu^2)} \tag{3}\label{3}$$
equations (2) and (3) allows you to prove that $$\Gamma^{(4)}(\frac{s}{\mu^2}, \lambda(\mu^2))=\Gamma^{(4)}(\lambda_R)\tag{4}\label{4}$$ that is the observable depends on its kinematics only through the RCC. What all this is telling us is that:
1)Renormalization breaks the IW (1) associated to scale invariance
2)The IW hence are modfied by renormalization, and are substituted by the RGE equation (in our case for physical observables) which differs by the naive one (1) for an anomaly factor proportional to the beta function.
3)The introduction of the renormalization scale allows to have a non trivial kinematic dependence, which at the bare level, due to dimensional analysis, wasn't possible. Due to the unobservability of the renormalization scale, this kinematic dependence by s can enter only in a certain way. So there is a kind of universal behaviour for the way the kinematic enters, which allows us to reabsorb it in the running coupling constant.
Another thing you can notice is that eq (4) tells you that the observable, if expandend around the RCC, is actually a series in $\lambda_R$ with adimensional coefficents. But this is exactly what we would say to hold $\lambda_\mu^2$ due the naive Ward identity $\mu^2\frac{\partial\Gamma^{(4)}(\lambda)}{\partial\mu^2}=0$.
This to me explains the conceptual relevance of the running coupling constant, in relation to the theory of scale invariance. How renormalization modifies (1), how (1) is corrected and how a "similar kind of scale invariance" (obviously scale invariance is still broken) consequence is obatined if the observale is rexpressed in terms of the RCC.
Now the more practical point of view of all this. Now you have seen that the RGE equation, as shown before, poses a great constraint on how the kinematic dependence can enter. What we can do is use this constraint to improve a fixed order calculation with a renormalization group improvement procedure. That is given a fixed order calculation you substitute the renormalized coupling with the RCC. This procedure allows you to include towers of logarithmic terms. This logarithmic terms to be obtained would require difficult calculation of Feynman diagrams at high perturbative orders. In this way you obtain them almost for free, because the beta function is a property of the theory, so once the $\beta^{(i)}$ are determined you know them for every process.
Take as an example the $\Gamma^{(4)}$ case, and imagine to do the RG improvement (RGI) at LO
$$\Gamma^{(4), LO}=\lambda(\mu^2)\rightarrow \Gamma^{(4), RGI}=\lambda_R$$
Suppose you have solved at the first order of the beta function $\beta(\lambda)=-\beta^{(0)}\lambda^2(1+O(\lambda))$ the differential equation for the RCC. You find that
$$\lambda_R(s)=\frac{\lambda_(\mu^2)}{1+\beta^{(0)}\lambda(\mu^2)ln(s/\mu^2)}\tag{5}\label{5}$$
and expanding the denominator, what you get is something like
$$\lambda_R=\lambda(\mu^2)+\sum_k\lambda^kln^{k-1}(s/\mu^2)$$
That is we obtain the LO preditcion from the amplitude, plus all the higher powers of logarithms at every fixed order N$^m$LO, which we call leading logarithms LL. So by substituting the renormalized CC with RCC you include all this LL terms. If you would have solved the eq for the RCC with $\beta^{(1)}$ you would find even the NLL terms that goes like $\sum_k\lambda^kln^{k-2}(s/\mu^2)$. It is evident that in the kinematic limit where $\lambda<<1$, in which a perturbative approach is reliable, the LL terms are dominant respect the NLL terms. Simply because the first goes as $\lambda(\lambda ln)^{m}$ and the latter as $\lambda^2(\lambda ln)^{m}$.
These logarithmic terms doesn't correspond to the whole fixed orders, are just some terms that would be found by an explicit calculation at the relative order. They can be more or less relevant to the magnitude of the whole respective fixed order, depending on the kinematic regions. But due to the RGI, once you have calculated the beta coefficents, which are universal, you can add them to your perturbative expansion.
There are situations in which they are more relevant for others. An example is an asymptotic free theory like QCD. In this case in the high energy limit the CC goes like $$\alpha_s \sim \frac{1}{ln(s/\mu^2)}$$. The LL terms then are exactly of the same order of the LO term because $\alpha_s ln(s/\mu^2) \sim 1$, they are not negligible. Actually they must be added to the perturbative series for this reason, but at the same time they kind of spoil the sense of the perturbative approach, as you get something that can goes like an infinite some of constants of the same order. So to maintain an adequate perturbative approach you must include this LL terms, but you have to choose the renormalization scale in a suitable way ( usually for what I know 2 or 0.5 times the scale of the process s) to avoid this instability. You can actually make this choice since chosing $\mu$ or $2\mu$ as your renormalization scale differs only by terms which are subleading, so in the case of LL these differences are NLL). You can see this by making this substituition $\mu \rightarrow 2\mu$ in eq (5). So in general it should be all consistent according to the accuracy of order N$^m$LO N$^m$LL you choose.
Hope this helps you.
