Consider the following **xtmelogit**.

use https://stats.idre.ucla.edu/stat/data/hsbdemo, clear xtmelogit honors i.female read || cid:Refining starting values: Iteration 0: log likelihood = -82.156921 Iteration 1: log likelihood = -80.132249 Iteration 2: log likelihood = -79.765159 Performing gradient-based optimization: Iteration 0: log likelihood = -79.765159 Iteration 1: log likelihood = -79.721813 Iteration 2: log likelihood = -79.72173 Iteration 3: log likelihood = -79.72173 Mixed-effects logistic regression Number of obs = 200 Group variable: cid Number of groups = 20 Obs per group: min = 7 avg = 10.0 max = 12 Integration points = 7 Wald chi2(2) = 6.39 Log likelihood = -79.72173 Prob > chi2 = 0.0410 ------------------------------------------------------------------------------ honors | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- 1.female | 1.203685 .498816 2.41 0.016 .2260232 2.181346 read | .0426387 .0579243 0.74 0.462 -.0708908 .1561682 _cons | -4.827711 2.919208 -1.65 0.098 -10.54925 .8938321 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ cid: Identity | sd(_cons) | 2.097709 .9227794 .8857332 4.96807 ------------------------------------------------------------------------------ LR test vs. logistic regression: chibar2(01) = 11.44 Prob>=chibar2 = 0.0004

Say that we wanted to know the predicted probabilities for males and females at five
different values of **read**, 30, 40, 50, 60 and 70. We might try using the **margins**
command.

margins female, at(read=(30(10)70)) vsquishdefault prediction is a function of possibly stochastic quantities other than e(b) r(498);

That didn’t work because **margins** can’t compute predicted values for models that
have both fixed and random components.

We can, however, get the predicted probabilities from just the fixed part of the model, like this.

margins female, at(read=(30(10)70)) predict(mu fixedonly) vsquishAdjusted predictions Number of obs = 200 Expression : Predicted mean, fixed portion only, predict(mu fixedonly) 1._at : read = 30 2._at : read = 40 3._at : read = 50 4._at : read = 60 5._at : read = 70 ------------------------------------------------------------------------------ | Delta-method | Margin Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- _at#female | 1 0 | .027962 .0353643 0.79 0.429 -.0413509 .0972748 1 1 | .0874748 .1006056 0.87 0.385 -.1097085 .2846582 2 0 | .0422023 .0352079 1.20 0.231 -.0268039 .1112085 2 1 | .1280315 .0900429 1.42 0.155 -.0484494 .3045123 3 0 | .0632231 .0416911 1.52 0.129 -.01849 .1449363 3 1 | .1836082 .0928378 1.98 0.048 .0016494 .365567 4 0 | .0936902 .0807071 1.16 0.246 -.0644929 .2518733 4 1 | .2562211 .1691374 1.51 0.130 -.075282 .5877243 5 0 | .1366968 .1661537 0.82 0.411 -.1889584 .462352 5 1 | .3454012 .3085949 1.12 0.263 -.2594336 .9502361 ------------------------------------------------------------------------------

For many purposes these probabilities from the fixed effects only will be all that we
will need and these probabilities could be graphed using **marginsplot**. If we are specifically interested in the estimated of probabilities that
include both fixed and random effects we can make use of the **predict** command.

First, we will estimate the predicted probabilities from the fixed and random parts of the model directly.

predict mu1, mu tabstat mu1, by(female)Summary for variables: mu1 by categories of: female female | mean -------+---------- male | .1909927 female | .3182445 -------+---------- Total | .2603449 ------------------

The values produced by **tabstat** give us the predicted probabilities separately for males
and females while **read** is held at its observed values. We can estimate these
same values in two steps by estimating the linear predictor for the random and fixed
effects separately.

predict re*, reffects // linear predictor for the random effects predict xb, xb // linear predictor for the fixed effects gen mu2 = 1 /(1+exp(-1*(xb + re1))) // compute probabilities using both fixed and random components tabstat mu2, by(female)Summary for variables: mu2 by categories of: female female | mean -------+---------- male | .1909927 female | .3182445 -------+---------- Total | .2603449 ------------------

As you can see the predicted probabilities computed this way are the same as when
we predicted **mu1** directly.
The term **(xb + re1)** combines the fixed and random linear predictors while
** 1 /(1+exp(-1*(xb + re1)))** converts the predictions to the probability metric.

We will now use the same approach to fix **read** at its mean value while letting
**female** vary as observed.

sum read replace read = r(mean) predict xb3, xb gen mu3 = 1 /(1+exp(-1*(xb3 + re1))) tabstat mu3, by(female)Summary for variables: mu3 by categories of: female female | mean -------+---------- male | .1588983 female | .3037404 -------+---------- Total | .2378372 ------------------

We can use this same method to compute predicted probabilities for each gender at the
five values of **read** of interest by computing a separate **xb** for each combination of
values. We can continue to use the the linear random predictor, **re1**, because the observations
remain nested in the same **cid**. We will do the computations in a **forvalues** loop.

forvalues i=30(10)70 { quietly replace read = `i' quietly replace female = 0 // begin with the males predict xbm`i' , xb gen mum`i' = 1 /(1+exp(-1*(xbm`i' + re1))) quietly replace female = 1 // now for the females predict xbf`i' , xb gen muf`i' = 1 /(1+exp(-1*(xbf`i' + re1))) } sum mum* muf*Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- mum30 | 200 .0878644 .1277633 .0062536 .452949 mum40 | 200 .1181626 .161502 .0095469 .5591283 mum50 | 200 .1545762 .1959134 .0145493 .660161 mum60 | 200 .1969638 .2284821 .0221143 .7484569 mum70 | 200 .2450655 .2570279 .0334791 .8200648 -------------+-------------------------------------------------------- muf30 | 200 .189036 .2229538 .0205396 .7339823 muf40 | 200 .2361479 .2523496 .0311209 .8086573 muf50 | 200 .2887174 .2762624 .0468924 .8661915 muf60 | 200 .3463831 .2931488 .0700785 .9083859 muf70 | 200 .4087083 .3014779 .1034842 .9382238

The variables that begin **mum** give the predicted values for the males for each of
the five different values of **read**. The **muf** do the same thing for the
females. We will plot these values by collapsing the data to a single observation and
then reshaping the data long.

collapse (mean) mum* muf* gen tid = 1 reshape long mum muf, i(tid) j(read) twoway (line mum read)(line muf read), legend(order(1 "male" 2 "female")) /// ytitle(predicted probability) /// title("Predicted probabilities with both" "fixed and random effects") /// name(mu, replace)

With this approach you can replace the observed values with the mean values or any value that you wish. So, this approach worked okay for random intercepts, how would it work if we include a random coefficient (slope).

use https://stats.idre.ucla.edu/stat/data/hsbdemo, clear xtmelogit honors i.female read || cid: female, cov(unstr)Refining starting values: Iteration 0: log likelihood = -81.308446 Iteration 1: log likelihood = -79.509044 Iteration 2: log likelihood = -79.233888 Performing gradient-based optimization: Iteration 0: log likelihood = -79.233888 Iteration 1: log likelihood = -79.203192 Iteration 2: log likelihood = -79.20305 Iteration 3: log likelihood = -79.20305 Mixed-effects logistic regression Number of obs = 200 Group variable: cid Number of groups = 20 Obs per group: min = 7 avg = 10.0 max = 12 Integration points = 7 Wald chi2(2) = 2.07 Log likelihood = -79.20305 Prob > chi2 = 0.3544 ------------------------------------------------------------------------------ honors | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- 1.female | 1.107571 .9144827 1.21 0.226 -.6847819 2.899924 read | .0548958 .0614466 0.89 0.372 -.0655374 .1753289 _cons | -5.459728 3.237309 -1.69 0.092 -11.80474 .8852801 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ cid: Unstructured | sd(female) | 1.270624 .877583 .3281837 4.919457 sd(_cons) | 2.104286 1.168272 .7088079 6.247133 corr(female,_cons) | -.1370656 1.016493 -.9741795 .9555909 ------------------------------------------------------------------------------ LR test vs. logistic regression: chi2(3) = 12.48 Prob > chi2 = 0.0059 Note: LR test is conservative and provided only for reference.predict mu1, mu tabstat mu1, by(female)Summary for variables: mu1 by categories of: female female | mean -------+---------- male | .1884487 female | .3188527 -------+---------- Total | .2595188 ------------------

To reproduce the predicted probabilities that include both fixed and random effects, we need to obtain the linear predictor of the random effects as follows:

predict re*, reffects // obtain the random effects des re*storage display value variable name type format label variable label -------------------------------------------------------------------------------------------- read byte %9.0g reading score re1 float %9.0g random effects for cid: female re2 float %9.0g random effects for cid: _cons

As you can see, **re1** is the linear predictor for the random coefficient while **re2**
linear the linear predictor for the random intercept.

Now we need to compute the interaction of **female** with **re1**, obtain the
linear predictor for the fixed effect and compute the probability.

gen fre1=female*re1 predict xb2, xb gen mu2 = 1 /(1+exp(-1*(xb2 + fre1 + re2))) tabstat mu2, by(female)Summary for variables: mu2 by categories of: female female | mean -------+---------- male | .1884487 female | .3188527 -------+---------- Total | .2595188 ------------------

Once again **mu1** and **mu2** are the same.

Next, with this knowledge, we are going to compute the predicted probabilities for
**read** for the values
30, 40 50 60 and 70 for both males and females, just like we did in the first section.
This time we need to compute **fre1** for both males and female. For male the
value of **female*re1** is zero. While for female the value of **female*re1** is
equal to **re1**.

forvalues i=30(10)70 { quietly replace read = `i' quietly replace female = 0 quietly replace fre1 = 0 // this value is 0 for males predict xbm`i' , xb gen mum`i' = 1 /(1+exp(-1*(xbm`i' +fre1 + re2))) quietly replace female = 1 quietly replace fre1 = re1 // this value is re1 for females predict xbf`i' , xb gen muf`i' = 1 /(1+exp(-1*(xbf`i' + fre1 + re2))) } sum mum* muf*Variable | Obs Mean Std. Dev. Min Max -------------+-------------------------------------------------------- mum30 | 200 .067521 .1073255 .0062084 .4335063 mum40 | 200 .1006205 .1453719 .010701 .5698889 mum50 | 200 .1441323 .1850342 .0183843 .6964301 mum60 | 200 .1988144 .2226292 .031409 .7988807 mum70 | 200 .2647876 .2545738 .0531617 .8730579 -------------+-------------------------------------------------------- muf30 | 200 .1651364 .2082666 .0140494 .5906218 muf40 | 200 .2210162 .252848 .0240784 .7141231 muf50 | 200 .2841479 .2869894 .040969 .8122126 muf60 | 200 .3544923 .3071029 .0688718 .8821978 muf70 | 200 .432452 .3113492 .1135287 .9283999

As before we will **collapse** the data then **reshape** them prior to plotting.

collapse (mean) mum* muf* gen tid = 1 reshape long mum muf, i(tid) j(read) twoway (line mum read)(line muf read), legend(order(1 "male" 2 "female")) /// ytitle(predicted probability) /// title("Predicted probabilities with both" "fixed and random effects") /// name(mu2, replace)

The two graphs that we generated for this page look very similar but inspection of the tables for the predicted probabilities show that adding the random coefficient to the model changes the predicted probabilities. Let’s check the two models using a likelihood-ratio test.

use https://stats.idre.ucla.edu/stat/data/hsbdemo, clear quietly xtmelogit honors i.female read || cid: female, cov(unstr) estimates store m1 quietly xtmelogit honors i.female read || cid: estimates store m2 lrtest m1 m2Likelihood-ratio test LR chi2(2) = 1.04 (Assumption: m1 nested in m2) Prob > chi2 = 0.5953 Note: The reported degrees of freedom assumes the null hypothesis is not on the boundary of the parameter space. If this is not true, then the reported test is conservative.

So, in this instance, the model with the random coefficient is not significantly better than the model with just the random intercept.

**Acknowledgements:** Thanks to Jurjen Iedema of the Netherlands Institute for Social Research / SCP
for suggesting this approach.