There is no way you can derive itthe spin interaction term $H_I$ from non-relativistic mechanics fromand using $\mathbf p,\mathbf r$ only. The spin is an intrinsic property of the electron and you have to postulate it. Here I propose three ways to convince yourself:
One option is to take it as it is. You can convince yourself in a heuristic manner that if you expand $(\mathbf p +e \mathbf A)^2$ with the gauge $\mathbf A = \mathbf B \times \mathbf r /2$ you get various terms, one of them is of the form:
$$H_L=\frac{\mu_{\rm B}}{\hbar}\mathbf B \cdot \mathbf L$$
you can say that the spin $\mathbf S$ is an angular momentum too so to obtain the true Hamiltonian, you have to do the following replacement:
$$H_L\to H_{L-S}=H_L+H_I=\frac{\mu_{\rm B}}{\hbar}\mathbf B \cdot (\mathbf L+ g \mathbf S)$$ where $g$ is some constant (the gyromagnetic factor).
Another possibility is to give up, go to relativistic Dirac equation and work your way to a non-relativistic equation, you will find that the true Hamiltonian (in the non-relativistic limit) is
$$H = \frac12 \frac{[2\mathbf S \cdot (\mathbf p + e\mathbf A)]^2}{2m\hbar}$$ if you expand that you get $H_I$ with $g=2$.
But there is a third possibility. If you have ever heuristically derived Dirac equation you'll know that it is derived by linearizing the relativistic dispersion relation $E^{\rm rel}= \sqrt{p^2c^2+m^2c^4}$ so that you only have linear derivatives ($p$ and not powers of it, $p^2,p^4,\cdots$). The third possibility consist on linearizing not the relativistic $E^{\rm rel}$ but linearizing your usual non-relativistic $E=p^2/2m$ instead. This is called Lévy-Leblond equation and leads to Pauli equation (Schrödinger+$H_I)$ when you try to solve it. Here is a link to some page that discusses it http://quantumchymist.blogspot.com/2014/04/is-spin-relativistic-effect-l-and-first.html but a formal derivation can be found in Greiner's Quantum Mechanics. Interestingly, it also provides $g=2$.
One option is to take it as it is, an interaction with an angular momentum. You can convince yourself that if you expand $(\mathbf p +e \mathbf A)^2$ with the gauge $\mathbf A = \mathbf B \times \mathbf r /2$ you get various terms, one of them is of the form:$$H_L=\frac{\mu_{\rm B}}{\hbar}\mathbf B \cdot \mathbf L\;.$$ As the spin $\mathbf S$ is also an angular momentum operator, to obtain the true Hamiltonian, you have to do the following replacement: $$H_L\to H_{L-S}=H_L+H_I=\frac{\mu_{\rm B}}{\hbar}\mathbf B \cdot (\mathbf L+ g \mathbf S)\;,$$ where $g$ is some constant (the gyromagnetic factor).
Another possibility is to give up the nonrelativistic origin. Take relativistic Dirac equation and work your way to a non-relativistic equation, you will find that the true Hamiltonian (in the non-relativistic limit) is $$H = \frac12 \frac{[2\mathbf S \cdot (\mathbf p + e\mathbf A)]^2}{2m\hbar}\,$$ if you expand that you get $H_I$ with $g=2$.
There is still a third possibility. If you have ever heuristically derived Dirac equation you'll know that it involves the linearization of the relativistic dispersion relation $E^{\rm rel}= \sqrt{p^2c^2+m^2c^4}$ imposing that you obtain an equation that only has linear derivatives ($p$ and not powers of it, $p^2,p^4,\cdots$). The alternative derivation of $H_I$ consists on linearizing, not the relativistic $E^{\rm rel}$, but your usual non-relativistic kinetic energy $E=p^2/2m$ instead. The resulting formula is called the Lévy-Leblond equation and naturally leads to Pauli equation (Schrödinger+$H_I)$ when you try to solve it. A step by step derivation can be found in Greiner's Quantum Mechanics. Interestingly, it also provides $g=2$.
Which is the best? To me all are equally ad hoc and heuristic,In my opinion the derivation should start from the most general Lagrangian and not a non-relativistic one. However, when one does not know the most general Lagrangian it is good to have different approaches to be convinced of the final result.