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`

.

CLICK HERE to find out more related problems solutions.