(Questions summarized at the bottom of the post.) For discrete BJT circuits, a biasing scheme like the one below is often described which helps reduce the sensitivity of \$I_C\$ to variations in the BJT's \$\beta\$ and the divider resistances \$R_{1,2}\$. (This figure is from Razavi's Fundamentals of Microelectronics, 2nd edition, p. 186).
Analysis of the above circuit shows that if the voltage divider carries sufficient current compared to the base current, the influence of the relation \$I_C=\beta I_B\$ is reduced compared to \$I_C=I_S(e^{V_{BE}/V_T}-1)\$, and so changes in \$\beta\$ are not as significant as they would be if, for example, the biasing omitted \$R_2, R_E\$.
However, I'm wondering about the practical variations observed in \$\beta\$ compared to those in \$I_S\$, which basic device physics seems to suggest are highly related parameters. It seems that the same physical processes which lead to variations in \$\beta\$ would also lead to similar variation in \$I_S\$, which would then translate directly to variations in \$I_C\$ since \$I_C \propto I_S \$. From Chenming Hu's book Modern Semiconductor Devices for Integrated Circuits (publicly available online), we have that
$$ \beta = \frac{D_B W_E N_E n_{iB}^2}{D_E W_B N_B n_{iE}^2} $$
$$ I_S = A_E q \frac{D_B}{W_B} \frac{n_{iB}^2}{N_B} $$
(Here \$D\$ is the diffusion constant in cm^2/s in the base or emitter; \$W\$ is the width in cm of the base or emitter; \$N\$ is the doping concentration in the base or emitter in 1/cm^3; \$n_i\$ is the intrinsic carrier concentration in the base or emitter in 1/cm^3.) On a surface level it seems that variations in \$\beta\$ are likely to directly translate to variations in \$I_S\$.
To summarize, is it the case that in real BJTs, \$\beta\$ tends to vary much more than \$I_S\$? Is there information available on the statistical distribution in the variation of both parameters in some example devices? (I haven't been able to find info comparing the two.) Is there a physical reason for why one would vary much more than the other? If they both vary, does the biasing strategy above actually reduce variations in \$I_S\$ as well?