# how do we reproduce average marginal effects from the xtlogit model?

There are a couple of ways to do this, but essentially the problem boils down to the fact

$$\Pr(y_{it}=1 \vert x_{it})=\int\Lambda(u_i + x_{it}’\beta)\cdot \varphi(0,\sigma_u^2) du_i$$

where $\varphi()$ is the normal density. In your code, you are effectively setting the random effect $u_i$ to zero (which is what predict(pu0) does). This sets the RE to its average, which may not be what you had it mind. Of course, $u_i$ is not observed or even estimated by xtlogit, re, so if you want to replicate what predict(pr) does, you need to integrate the random effect out to get the unconditional probability using the estimated variance.

One way to do this in Stata is to use the user-written integrate command to do one dimensional numerical integration like this:

webuse union, clear
xtset idcode year
xtlogit union age, nolog
margins, at(age=(16)) predict(pr)
margins, dydx(*) at(age=16) predict(pr)
capture ssc install integrate
/* phat at age 16 */
integrate, f(invlogit(x - 3.0079682 + .01929225*16)*normalden(x,0,2.477537654945496)) l(-10) u(10) vectorise quadpts(1000)
/* ME at age 16 */
integrate, f(invlogit(x - 3.0079682 + .01929225*16)*(1-invlogit(x - 3.0079682 + .01929225*16))*(.01929225)*normalden(x,0,2.477537654945496)) l(-10) u(10) vectorise quadpts(1000)


You can probably get even better precision if you use the actual coefficients (like _b[_cons] and e(sigma_u)) instead of pasting the values in.

There might be a more efficient way to do this using Mata, Python, or even with a simulation approach, but I will leave that up to you to work out. You can also get posterior modal estimates of the REs using xtmelogit.