The only trouble is that this is not a unique prescription. The point is to "gauge" the Pauli equation. The analysis by Levy-Leblond can be shortened a bit. The basic point, why he gets the correct gyro-factor of 2 is that he takes a first-order differential equation representation of the wave function's dynamics as Dirac did for his equation. This derivation can be much shortened by just taking the analysis of the unitary ray representations of the Galilei group.
Non-relativistically, i.e., from the analysis of the Galilei group you get
[tex]\hat{H}=\frac{\hat{\vec{p}}^2}{2m}[/tex]
as the free-particle Hamiltionian of particles of any spin. If you "gauge" this, using the principle of minimal coulpling, you don't get the correct coupling of the spin to the em. field.
Now you can use the following trick. You just write
[tex]\hat{H}=\frac{(\vec{\sigma} \cdot \hat{\vec{p}})\cdot (\vec{\sigma} \cdot \hat{\vec{p}})}{2m},[/tex]
which is the same as above, because the momentum components commute, and the Pauli matrices generate the Clifford algebra of 3D Euclidean space, i.e., fullfil the anticommutation relations
[tex]\{\sigma_j,\sigma_k \}=2 \delta_{jk}.[/tex]
If you now do the "gauging" via minimal coupling, you find
[tex]\displaystyle{\hat{H}_{\text{em}}=\frac{[\vec{\sigma} \cdot (\hat{\vec{p}}-q \vec{A}(t,\hat{\vec{x}}))] \cdot [\vec{\sigma} \cdot (\hat{\vec{p}}-q \vec{A}(t,\hat{\vec{x}}))]}{2m} + q \Phi(t,\hat{\vec{x}})},[/tex]
and this gives a different result than the naive "gauging", because the momenta and the vector potential components don't commute. Of course [itex]\Phi[/itex] and [itex]\vec{A}[/itex] are the scalar and vector potential of the electromagnetic field.
Multiplying out the above Hamiltonian you get the correct Pauli Hamiltonian for spin-1/2-particles:
[tex]\displaystyle{\hat{H}_{\text{em}}=\frac{[\hat{\vec{p}}-q \vec{A}(t,\hat{\vec{x}})]^2}{2m} -\frac{q}{2m} g_s \hat{\vec{S}} \cdot \vec{B}(t,\hat{\vec{x}})+q \Phi(t,\hat{\vec{x}}) \quad \text{with} \quad g_s=2.}[/tex]
Without the trick to introduce the vector potential into the non-naive Hamiltonian you wouldn't have gotten the coupling of the spin with the magnetic field [itex]\propto \vec{B} \cdot \hat{\vec{S}}[/itex] with [itex]\vec{B}=\vec{\nabla} \times \vec{A}[/itex] at all.
That this is the "correct" equation for an elementary particle can only be derived by taking the non-relativistic limit of the Dirac equation.
What's of course right is that spin is not a relativistic feature per se but it emerges as "naturally" from the analysis of the unitary ray representations of the Galilei group as it arises from the analogous analysis for the Poincare group.