The spin-only formula
$$\mu_\mathrm{so} = \mu_\mathrm{B} \sqrt{n(n+2)} = \mu_\mathrm{B}\cdot 2\sqrt{S(S+1)}$$
is usually a good first approximation to calculate the magnetic moment of transition metal complexes. (I am aware of its limitations.) The usual derivation of this formula involves the Zeeman splitting of the ground state into its separate $M_S$ states under the influence of a magnetic field, followed by evaluating the average magnetic moment using a Boltzmann distribution.
However, if I were to write that $\vec{\mu} = \gamma_\mathrm{e} g_\mathrm{e} \vec{S}$, then I could arrive at the simpler "derivation"
$$\begin{align} |\vec{\mu}| &= |\gamma_\mathrm{e}|\cdot g_\mathrm{e} \cdot |\vec{S}| \\ &= \frac{e}{2m_\mathrm{e}} \cdot 2 \cdot \sqrt{S(S+1)}\hbar \\ &= \mu_\mathrm{B}\cdot 2\sqrt{S(S+1)} \end{align}$$
Is this a coincidental result? I'm assuming there is something wrong with the second "derivation", or even worse, it is just an entirely unphysical interpretation, since I am using $S$ instead of $S_z$. What is the mistake?