Setting of renormalization scale in field theory calculations

In dimensional regularization an arbitrary mass parameter $\mu$ must be introduced in going to $4-\epsilon$ dimensions. I am trying to understand to what extent this parameter can be eliminated from physical observables.

Since $\mu$ is arbitrary, physical quantities such as pole masses and scattering amplitudes must be independent of it. Nevertheless at any fixed order in perturbation theory these quantities contain residual $\mu$-dependence. One expects this dependence to decrease at higher orders in perturbation theory.

For concreteness, consider dimensional regularization with minimal subtraction of $\phi^4$ theory, which has bare Lagrangian

$\mathcal{L}_B = \frac{1}{2}(\partial \phi_B)^2 - \frac{1}{2}m_B^2 \phi_B^2 - \frac{\lambda}{4!}\phi_B^4$

Here are the 1-loop expressions for the physical mass $m_P$ and 4-point coupling $\lambda_P \equiv (\sqrt{Z})^4 \Gamma^{(4)}$ after minimal subtraction of poles and taking $\epsilon \to 0$:

$$\begin{equation} m_P^2 = m_R^2 \left\{1 + \frac{\lambda_R}{2(4\pi)^2}\left[\log\left(\frac{m_R^2}{4\pi\mu^2}\right)\right] + \gamma - 1 \right\} \end{equation}$$

$$\lambda_P = \lambda_R + \frac{3\lambda_R^2}{2(4\pi)^2}\left[\log\left(\frac{m_R^2}{4\pi\mu^2}\right) + \gamma - 2+\frac{1}{3}A\left(\frac{m_R^2}{s_E},\frac{m_R^2}{t_E},\frac{m_R^2}{u_E}\right) \right]$$

where $A\left(\frac{m_R^2}{s_E},\frac{m_R^2}{t_E},\frac{m_R^2}{u_E}\right) = \sum_{z_E = s_E,t_E,u_E} A\left(\frac{m_R^2}{z_E}\right)$ and $A(x) \equiv \sqrt{1+4x}\log\left(\frac{\sqrt{1+4x}+1}{\sqrt{1+4x}-1}\right)$.

Both of these quantities ($m_P$ and $\lambda_P$) are physically observable.

Suppose we conduct an experiment at a reference momentum $p_{E0}\equiv(s_{E0},t_{E0},u_{E0})$ and make measurements of the pole mass and 4-point coupling with the result $\lambda_{P0}, m_{P0}$. We now have a system of two equations in the three unknowns ($\lambda_R,m_R,\mu$). This means that in principle I can solve for $\lambda_R = \lambda_R(\mu)$ and $m_R = m_R(\mu)$.

I would now like to make a prediction for the 4-point amplitude at a different momentum $p_{E}' \neq p_{E0}$. Since I have two equations in three unknowns I need to guess a suitable value for $\mu$ (say $\mu' = \sqrt{s_{E}'}$) which allows me to fix $\lambda_R$ and $m_R$ and then calculate $\lambda_P (p_E')$. This procedure seems quite ad hoc to me because a different (arbitrary) choice of $\mu$ (e.g. $\mu'/2$ or $2\mu'$) will lead a different physical answer (albeit only logarithmically different).

From what I can gather from the literature, the problem of determining the renormalization scale I described above is a genuine problem in actual calculations of QCD (e.g. http://arxiv.org/abs/1302.0599) which leads theorists to introduce so-called "systematic uncertainties".

What concerns me is that I haven't been able to find any mention of this problem in any textbook on quantum field theory that deals with QED or QCD (anyone know of a reference?). Since this problem appears already in arguably the simplest QFT of $\lambda\phi^4$, I would expect it also to occur in e.g. 1-loop calculations of Bhabha scattering but I haven't been able to find any mention of it in this context.

Does anyone know how this problem is dealt with in real loop calculations (e.g. at LEP or LHC?)

Also, I would be interested to know if there is any analogue of this problem in condensed matter theory.