It seems that some people liked this question so I shall post my thoughts so far. I don't have a definitive answer, but I did get some interesting results.

Let $\rho_m(\boldsymbol r)$ and $\rho_e(\boldsymbol r)$ be the mass and charge densities of the electron. The $g$ factor is given by $$ g=\frac{m}{e}\frac{\int\mathrm d\boldsymbol r\ r^2\sin\theta\ \rho_m(\boldsymbol r)}{\int\mathrm d\boldsymbol r\ r^2\sin\theta\ \rho_e(\boldsymbol r)} \tag{1} $$

From this, it's easy to see that if $\rho_m\propto \rho_e$, we get $g=1$. This means that if we have a solid sphere with constant charge density and constant mass density, the $g$ factor is 1; al hollow sphere with surface charge has also $g=1$. If we want $g\neq 1$ we must take a charge density that is not proportional to the mass density.

The first model that comes to mind is to take a volume mass density and a surface charge density, that is, a filled sphere with its charge on the surface: \begin{align} \rho_m&=\frac{m}{V}\Theta(R-r)\\ \rho_e&=\frac{e}{S}\delta(r-R)\tag{2} \end{align} where $V=\frac{4}{3}\pi R^3$ and $S=4\pi R^2$. If we plug these functions into $(1)$ we get $g=5/3$ as already acknowledged by Ilja and Anubhav. This means that Arabatzis' and Pais' claims are inaccurate: this model does not predict $g=2$ but $g=1.67$ instead.

To go a step further, we may take the same model before, but with a different mass and charge radii, that is, \begin{align} \rho_m&=\frac{m}{V}\Theta(R_m-r)\\ \rho_e&=\frac{e}{S}\delta(R_e-r)\tag{3} \end{align} with $R_m\neq R_e$. In this case, we find $g=5R_e^2/3R_m^2$, which equals 2 if $R_e=1.095 R_m$. This model seems highly artificial though.

The next possible example could be to take exponential densities, which could be the result of some kind of screening at some fundamental level: \begin{align} \rho_m&\propto\exp\left[-\frac{r^2}{R_m^2}\right]\\ \rho_e&\propto\exp\left[-\frac{r}{R_e}\right]\tag{4} \end{align} from which we find $g=8R_e^2/R_m^2$; if we take $R_m=2R_e$ we get $g=2$. This is still very artificial but there might be some electrostatic model that is able to accommodate this.

Other possible models could consist of non-spherical densities, such as cylinders or string-like wires. I leave to the reader to explore this models. In any case, it is clear that the most natural models don't predict $g=2$, and it's not easy to find another one that fixes this while not getting too ad-hoc. But it is possible to write down exotic models with tunable parameters so as to get $g=2$, which means at least that $g=2$ is achievable at the classical level.

It seems that some people liked this question so I shall post my thoughts so far. I don't have a definitive answer, but I did get some interesting results.

Let $\rho_m(\boldsymbol r)$ and $\rho_e(\boldsymbol r)$ be the mass and charge densities of the electron. The $g$ factor is given by $$ g=\frac{m}{e}\frac{\int\mathrm d\boldsymbol r\ r^2\sin\theta\ \rho_m(\boldsymbol r)}{\int\mathrm d\boldsymbol r\ r^2\sin\theta\ \rho_e(\boldsymbol r)} \tag{1} $$

From this, it's easy to see that if $\rho_m\propto \rho_e$, we get $g=1$. This means that if we have a solid sphere with constant charge density and constant mass density, the $g$ factor is 1; al hollow sphere with surface charge has also $g=1$. If we want $g\neq 1$ we must take a charge density that is not proportional to the mass density.

The first model that comes to mind is to take a volume mass density and a surface charge density, that is, a filled sphere with its charge on the surface: \begin{align} \rho_m&=\frac{m}{V}\Theta(R-r)\\ \rho_e&=\frac{e}{S}\delta(r-R)\tag{2} \end{align} where $V=\frac{4}{3}\pi R^3$ and $S=4\pi R^2$. If we plug these functions into $(1)$ we get $g=5/3$ as already acknowledged by Ilja and Anubhav. This means that Arabatzis' and Pais' claims are inaccurate: this model does not predict $g=2$ but $g=1.67$ instead.

To go a step further, we may take the same model before, but with a different mass and charge radii, that is, \begin{align} \rho_m&=\frac{m}{V}\Theta(R_m-r)\\ \rho_e&=\frac{e}{S}\delta(R_e-r)\tag{3} \end{align} with $R_m\neq R_e$. In this case, we find $g=5R_e^2/3R_m^2$, which equals 2 if $R_e=1.095 R_m$. This model seems highly artificial though.

The next possible example could be to take exponential densities, which could be the result of some kind of screening at some fundamental level: \begin{align} \rho_m&\propto\exp\left[-\frac{r^2}{R_m^2}\right]\\ \rho_e&\propto\exp\left[-\frac{r}{R_e}\right]\tag{4} \end{align} from which we find $g=8R_e^2/R_m^2$; if we take $R_m=2R_e$ we get $g=2$. This is still very artificial but there might be some electrostatic model that is able to accommodate this.

Other possible models could consist of non-spherical densities, such as cylinders or string-like wires. I leave to the reader to explore this models. In any case, it is clear that the most natural models don't predict $g=2$, and it's not easy to find another one that fixes this while not getting too ad-hoc. But it is possible to write down exotic models with tunable parameters so as to get $g=2$, which means at least that $g=2$ is achievable at the classical level.