Let’s start off by down loading some data and running a binary mixed model with random intercepts and slopes.
use https://stats.idre.ucla.edu/stat/data/hsbdemo, clear
meqrlogit honors read i.female || cid: read, var covar(unstr)
Refining starting values:
Iteration 0: log likelihood = -119.30049 (not concave)
Iteration 1: log likelihood = -105.82055 (not concave)
Iteration 2: log likelihood = -88.907174
Performing gradient-based optimization:
Iteration 0: log likelihood = -88.907174 (not concave)
Iteration 1: log likelihood = -84.544443 (not concave)
[output omitted]
Iteration 12: log likelihood = -78.268094
Iteration 13: log likelihood = -78.268094
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) = 8.29
Log likelihood = -78.268094 Prob > chi2 = 0.0159
------------------------------------------------------------------------------
honors | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
read | .1084052 .0671505 1.61 0.106 -.0232074 .2400177
1.female | 1.309964 .5183077 2.53 0.011 .2941 2.325829
_cons | -8.18226 3.894754 -2.10 0.036 -15.81584 -.548682
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
-----------------------------+------------------------------------------------
cid: Unstructured |
var(read) | .0239163 .0244789 .0032171 .1777943
var(_cons) | 102.4491 96.65305 16.12367 650.9578
cov(read,_cons) | -1.55161 1.532405 -4.555068 1.451848
------------------------------------------------------------------------------
LR test vs. logistic regression: chi2(3) = 14.35 Prob > chi2 = 0.0025
Note: LR test is conservative and provided only for reference.
Let’s try to use the margins command to get the predictive means (probabilities) for males and females holding the reading score at 50.
margins female, at(read=50) predict(mu) prediction is a function of possibly stochastic quantities other than e(b)
It doesn’t work because the model has random effects in addition to the fixed effects. To get margins to work we need to include the predict(mu fixed) option.
margins female, at(read=50) predict(mu fixed)
Adjusted predictions Number of obs = 200
Expression : Predicted mean, fixed portion only, predict(mu fixed)
at : read = 50
------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
female |
0 | .0594124 .0466851 1.27 0.203 -.0320887 .1509134
1 | .1896881 .1177394 1.61 0.107 -.0410769 .4204531
------------------------------------------------------------------------------
So, if margins won’t compute predictive margins with random effects we will have to compute them manually. We will begin with the easier task of computing predicted probabilities that include both the fixed and random effects. Of course, there is an option in predict that will do this. We will use predict, mu to check the results of our computation.
predict mu, mu // mu contains both fixed effects and random effects
Now we will replicate mu using the predicted random effects for both the intercept and the slope. These are the Best Linear Unbiased Predictors (BLUPs) also called empirical Bayes predictors.
predict re*, reffects
. describe re1 re2
storage display value
variable name type format label variable label
--------------------------------------------------------------------------------------------------
re1 float %9.0g random effects for cid: read
re2 float %9.0g random effects for cid: _cons
We will use these BLUPs along with the coefficients from our model. The re2 values get added to the constant from the model and the re1 values are added to the coefficient for read. Since the variable rxb is in the log odds metric we need to exponentiate it to obtain predicted probabilities.
generate rxb = _b[1.female]*female + (_b[read]+re1)*read + (_b[_cons]+re2)
generate mu2 = exp(rxb)/(1+exp(rxb))
list mu mu2 in 1/20
+---------------------+
| mu mu2 |
|---------------------|
1. | .0068868 .0068868 |
2. | .0018676 .0018676 |
3. | .0041634 .0041634 |
4. | .0030219 .0030219 |
5. | .0041634 .0041634 |
|---------------------|
6. | .0244843 .0244843 |
7. | .0011538 .0011538 |
8. | .0239256 .0239256 |
9. | .0152577 .0152577 |
10. | .0018676 .0018676 |
|---------------------|
11. | .0068868 .0068868 |
12. | .0039202 .0039202 |
13. | .0142891 .0142891 |
14. | .0088657 .0088658 |
15. | .0075425 .0075425 |
|---------------------|
16. | .0088119 .0088119 |
17. | .0024226 .0024226 |
18. | .0232304 .0232304 |
19. | .0017392 .0017392 |
20. | .0017392 .0017392 |
+---------------------+
The listing above shows that our manually computed mu2 is equal to the value mu generated by the predict command.
Let’s return to computing the predictive margins with random effect for males and females while holding reading at 50. In the code below rxb0 and rxb1 are the predictive margins in log odds for males and females respectively. To get the predictive margins in the probability metric we exponentiate the values to get p0 and p1. Again, the zero term is for males and the one term is for females.
generate rxb0 = _b[1.female]*0 + (_b[read]+re1)*50 + (_b[_cons]+re2)
generate rxb1 = _b[1.female]*1 + (_b[read]+re1)*50 + (_b[_cons]+re2)
generate p0 = exp(rxb0)/(1+exp(rxb0))
generate p1 = exp(rxb1)/(1+exp(rxb1))
summarize p0 p1
Variable | Obs Mean Std. Dev. Min Max
-------------+--------------------------------------------------------
p0 | 200 .1755775 .2356444 .0125389 .829953
p1 | 200 .3156575 .303865 .0449446 .9476115
These probabilities are rather different from the values we computed earlier using just the fixed effects. This highlights the fact that estimating predicated values while averaging over the fixed effects (e.g. predict) is not the same as estimating predicted values assuming the random effect is zero (e.g. margins) If we look at p0 and p1 within cid equal 1 we see that all the values for each variable are the same.
list p0 p1 if cid==1
+---------------------+
| p0 p1 |
|---------------------|
1. | .0239256 .0832776 |
2. | .0239256 .0832776 |
3. | .0239256 .0832776 |
4. | .0239256 .0832776 |
5. | .0239256 .0832776 |
|---------------------|
6. | .0239256 .0832776 |
7. | .0239256 .0832776 |
8. | .0239256 .0832776 |
9. | .0239256 .0832776 |
10. | .0239256 .0832776 |
|---------------------|
11. | .0239256 .0832776 |
+---------------------+
This means that the mean values computed by summarize above are influenced by the number of observations within each cid. That is, larger cids influence the mean more than smaller cids. To get around this we can use just one value from each cid. This will produce a more balanced predicted values. Here is a way to do this.
bysort cid: gen seq=_n
summarize p0 p1 if seq==1
Variable | Obs Mean Std. Dev. Min Max
-------------+--------------------------------------------------------
p0 | 20 .1790172 .2416508 .0125389 .829953
p1 | 20 .3212767 .3125198 .0449446 .9476115
Now we are ready to tackle a more complex model and set of predictive margins. We will add socst to our model as a covariate. But first, we will drop the variables we created above.
drop mu-p1
xtmelogit honors read female socst || cid: read, var covar(unstr)
Refining starting values:
Iteration 0: log likelihood = -118.63766 (not concave)
Iteration 1: log likelihood = -105.09158 (not concave)
Iteration 2: log likelihood = -88.297585
Performing gradient-based optimization:
Iteration 0: log likelihood = -88.297585 (not concave)
Iteration 1: log likelihood = -83.769166 (not concave)
[output omitted]
Iteration 8: log likelihood = -77.974955
Iteration 9: log likelihood = -77.974949
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(3) = 9.12
Log likelihood = -77.974949 Prob > chi2 = 0.0277
------------------------------------------------------------------------------
honors | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
read | .0980242 .0658639 1.49 0.137 -.0310668 .2271151
female | 1.263858 .514295 2.46 0.014 .2558583 2.271858
socst | .0243635 .0317934 0.77 0.443 -.0379504 .0866774
_cons | -8.900997 3.842395 -2.32 0.021 -16.43195 -1.370041
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
-----------------------------+------------------------------------------------
cid: Unstructured |
var(read) | .0209259 .0228458 .0024626 .1778206
var(_cons) | 89.13344 89.85624 12.35767 642.902
cov(read,_cons) | -1.353323 1.426992 -4.150175 1.443529
------------------------------------------------------------------------------
LR test vs. logistic regression: chi2(3) = 12.38 Prob > chi2 = 0.0062
Note: LR test is conservative and provided only for reference.
predict re*, reffects
This time we want the predictive margins for males and females holding reading at values running from 30 to 70 in increments of five while leaving socst free to vary. We will do this with some temporary variables and three matrices: m predictive margins for males, f predictive margins for females and d the difference in predictive margins between males and females. Because of the local variables in the code below you will need to run the block of code all at one time, not one line at a time.
*** begin block ***
tempvar rxb0 rxb1 p0 p1
generate `rxb0' = . // initialize temp var
generate `rxb1' = . // initialize temp var
generate `p0' = . // initialize temp var
generate `p1' = . // initialize temp var
matrix m = J(9,1,0) // initialize matrix for males
matrix f = J(9,1,0) // initialize matrix for females
matrix d = J(9,1,0) // initialize matrix for difference
local i = 1 // initialize local macro
forvalues j=30(5)70 {
quietly replace `rxb0' = _b[female]*0 + (_b[read]+re1)*`j' + _b[socst]*socst + _b[_cons]+re2
quietly replace `rxb1' = _b[female]*1 + (_b[read]+re1)*`j' + _b[socst]*socst + _b[_cons]+re2
quietly replace `p0' = exp(`rxb0')/(1+exp(`rxb0'))
quietly sum `p0'
matrix m[`i',1] = r(mean)
quietly replace `p1' = exp(`rxb1')/(1+exp(`rxb1'))
quietly sum `p1'
matrix f[`i',1] = r(mean)
matrix d[`i',1] = f[`i',1] - m[`i',1]
local i = `i' + 1 // increment local macro
}
*** end block ***
We can look at the predictive margins for males by listing matrix m.
matrix list m
m[9,1]
c1
r1 .17859283
r2 .17596309
r3 .17393653
r4 .17354885
r5 .17695161
r6 .18797502
r7 .21236354
r8 .2572271
r9 .32805038
This first value above is the predictive margins for males with a reading score of 30 averaged across the observed values of socst. So, it does not matter that the number of observations differ within cid because the values of socst vary within cid.
Next, we will save the matrices as observations to our dataset.
matrix r = (303540455055606570)
svmat r
svmat m
svmat f
svmat d
list r1 m1 f1 d1 in 1/9
+-------------------------------------+
| r1 m1 f1 d1 |
|-------------------------------------|
1. | 30 .1785928 .2603539 .0817611 |
2. | 35 .1759631 .2647262 .0887631 |
3. | 40 .1739365 .2729119 .0989754 |
4. | 45 .1735488 .2877812 .1142324 |
5. | 50 .1769516 .3135495 .1365979 |
|-------------------------------------|
6. | 55 .187975 .3555776 .1676026 |
7. | 60 .2123635 .4193873 .2070238 |
8. | 65 .2572271 .5054632 .2482362 |
9. | 70 .3280504 .5996285 .2715782 |
+-------------------------------------+
Now we are ready to graph the predictive margins for males and females.
twoway (line m1 r1)(line f1 r1), ytitle(Predictive Margins) xtitle(Reading Score) ///
legend(order(1 "Males" 2 "Females")) name(predictive, replace) ylabel(0(.2).6)
And finally, we can plot the difference between males and females in their predictive margins.
twoway (line d1 r1), ytitle(Difference in Predictive Margins) ///
xtitle(Reading Score) name(difference, replace)
This general approach can be extended to more complex model with a concomitant increase in the complexity of our code.