-margins- after -melogit-

[Jeph Herrin]

I am trying to replicate manually the results of -margins- after estimating a model using -melogit-. Here are my data:

. d y x group

              storage   display    value
variable name   type    format     label      variable label
--------------------------------------------------------------------------------------------
y               double  %12.0g                
x               double  %12.0g                
group           str6    %6s

Let me tag the groups, then run the model (which converges in 5 iterations):

. egen tag=tag(group)

. melogit y i.x || group:

Fitting fixed-effects model:

[...]

Mixed-effects logistic regression               Number of obs     =      2,087
Group variable:           group                 Number of groups  =        328

                                                Obs per group:
                                                              min =          1
                                                              avg =        6.4
                                                              max =         54

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(1)      =       4.50
Log likelihood = -335.56699                     Prob > chi2       =     0.0340
------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
             |
         1.x |   1.154201   .5442717     2.12   0.034     .0874484    2.220954
       _cons |  -4.972222   .4338917   -11.46   0.000    -5.822634    -4.12181
-------------+----------------------------------------------------------------
group        |
   var(_cons)|   5.033932   1.512757                      2.793259    9.072011
------------------------------------------------------------------------------
LR test vs. logistic model: chibar2(01) = 106.07      Prob >= chibar2 = 0.0000

Then calculate the marginal probabilities of x, conditioned on the random effects:

. margins x, predict(mu)

Adjusted predictions                            Number of obs     =      2,087
Model VCE    : OIM

Expression   : Marginal predicted mean, predict(mu)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           x |
          0  |   .0420558   .0080066     5.25   0.000     .0263632    .0577484
          1  |   .0902771   .0272681     3.31   0.001     .0368325    .1437217
------------------------------------------------------------------------------

So far so good. But then I wanted to demonstrate to someone how these were calculated, so I am trying to do it manually; first get the random effects for each group, then calculate the probabilities:

. 
. predict re, reffects
(calculating posterior means of random effects)
(using 7 quadrature points)

. gen prx0=invlogit(_b[_cons]+re)

. gen prx1=invlogit(_b[_cons]+re+_b[1.x])

. sum prx* if tag

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
        prx0 |        328    .0235069    .0703413   .0015256   .6677329
        prx1 |        328    .0537861    .1228471   .0048224   .8643794

. sum prx*

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
        prx0 |      2,087    .0274458    .0721269   .0015256   .6677329
        prx1 |      2,087    .0630935    .1330303   .0048224   .8643794

Why are these so different from the output of -margins-? Notably, if do the same thing using -margins x, predict(fixedonly)- and omit the random effects in my calculations, they agree.

thanks,
Jeph

[Alfonso Sánchez-Peñalver]

I believe you can find your answer in this FAQ at UCLA. It uses xtmelogit which is the old version of the melogit estimator you're using. It is a direct explanation of how to calculate those probabilities.

[Jeph Herrin]

Alfonso -

I checked my calculations against those on the UCLA page, and confirmed that my calculations are the same.

cheers,
Jeph

[Jeff Pitblado (StataCorp)]

Jeph's manual calculation is computing the average predicted probability using empirical Bayes' predicted values as plug-in values for the random effects.

margins is averaging the predicted probabilities that are marginal with respect to the random effects. To reproduce margins, Jeph can try the following:

rename x x_orig
gen x = 0
predict mprx0, mu marginal
replace x = 1
predict mprx1, mu marginal
drop x
rename x_orig x
sum mprx*

Here is an example using a dataset from the margins manual entry.

. webuse margex
(Artificial data for margins)

. gen x = sex

. keep yc x

. gen group = floor((_n-1)/30) + 1

. 
. melogit y i.x || group:

Fitting fixed-effects model:

Iteration 0:   log likelihood =  -1556.705  
Iteration 1:   log likelihood = -1553.7949  
Iteration 2:   log likelihood = -1553.7931  
Iteration 3:   log likelihood = -1553.7931  

Refining starting values:

Grid node 0:   log likelihood = -1595.1929

Fitting full model:

Iteration 0:   log likelihood = -1595.1929  (not concave)
Iteration 1:   log likelihood = -1553.4511  
Iteration 2:   log likelihood = -1553.1596  
Iteration 3:   log likelihood =  -1553.158  
Iteration 4:   log likelihood =  -1553.158  

Mixed-effects logistic regression               Number of obs     =      3,000
Group variable:           group                 Number of groups  =        100

                                                Obs per group:
                                                              min =         30
                                                              avg =       30.0
                                                              max =         30

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(1)      =     160.93
Log likelihood =  -1553.158                     Prob > chi2       =     0.0000
------------------------------------------------------------------------------
          yc |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
             |
         1.x |   1.336689   .1053693    12.69   0.000     1.130169    1.543209
       _cons |  -1.931789   .0853559   -22.63   0.000    -2.099084   -1.764495
-------------+----------------------------------------------------------------
group        |
   var(_cons)|   .0354421   .0353623                      .0050146    .2504973
------------------------------------------------------------------------------
LR test vs. logistic model: chibar2(01) = 1.27        Prob >= chibar2 = 0.1299

. 
. margins x, predict(mu)

Adjusted predictions                            Number of obs     =      3,000
Model VCE    : OIM

Expression   : Marginal predicted mean, predict(mu)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           x |
          0  |   .1280114   .0090849    14.09   0.000     .1102053    .1458175
          1  |   .3566212   .0138608    25.73   0.000     .3294546    .3837879
------------------------------------------------------------------------------

. 
. rename x xhold

. gen x = 0

. predict mprx0, mu marginal
(using 7 quadrature points)

. replace x = 1
(3,000 real changes made)

. predict mprx1, mu marginal
(using 7 quadrature points)

. drop x

. rename xhold x

. sum mprx*

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
       mprx0 |      3,000    .1280114           0   .1280114   .1280114
       mprx1 |      3,000    .3566212           0   .3566212   .3566212

Notice that since x is fixed and the predictions are marginal with respect to the random effects, the marginal predictions are constant.

[Jeph Herrin]

Jeff,

Thanks for the reply and the explanation. This is very helpful. I had actually tried using -predict mu, mu marginal- but must have done so incorrectly because it gave an error.

cheers,
Jeph

[Caroline Thirukumaran]

I am curious if the codes suggested by Jeff are applicable to Stata 14 only or could these be used in Stata 13 as well.

I am using Stata 13, and I am particularly interested in using

 .margins x, predict(mu)

following a -melogit- model that I am estimating. However, when I use the code, I get an error message "prediction is a function of possibly stochastic quantities other than e(b)". I have reviewed the UCLA webpage that Alfonso has referenced above, but the webpage does not provide direction on calculating the statistical significance of the marginal effect. I would greatly appreciate guidance on calculating the statistical significance for the marginal effect.

Many thanks,
Caroline