Manual calculation of average marginal effects after melogit

[Dario Maimone Ansaldo Patti]

Dear All,

I am trying to calculate the average marginal effects (ame) after running a multilevel logit model. However, it takes ages, even using the noci option. The point is that I have around 300,000 observations nested in 90 groups. While searching for a solution to speed up the process, I discovered a nice post by Jeff Pitblado (StataCorp) , which explains the procedure adopted by margins to calculate the ame and, more importantly, the standard errors using delta method after a logit estimation. Then I thought to adjust the procedure for calculating the marginal effects after melogit. The reason to calculate ame manually is that it can allow me to save a lot og time. Suppose I run the following model (results are not important. In fact I am creating a fake new variable) and I am interested in the ame of age:

webuse bangladesh, clear

set seed 12345

gen x2 = runiform(1, 100)

melogit c_use x2 age || district:

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(2)      =       2.56
Log likelihood = -1265.7725                     Prob > chi2       =     0.2787
------------------------------------------------------------------------------
       c_use | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
          x2 |  -.0002711   .0016756    -0.16   0.871    -.0035552     .003013
         age |   .0085437   .0053778     1.59   0.112    -.0019966    .0190841
       _cons |  -.5264599   .1220406    -4.31   0.000     -.765655   -.2872648
-------------+----------------------------------------------------------------
district     |
   var(_cons)|   .2522657   .0804515                      .1350196    .4713241
------------------------------------------------------------------------------

margins, dydx(age)

Average marginal effects                                 Number of obs = 1,934
Model VCE: OIM

Expression: Marginal predicted mean, predict()
dy/dx wrt:  age

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         age |   .0018937   .0011888     1.59   0.111    -.0004363    .0042237
------------------------------------------------------------------------------

As in Jeff's post, I save the two Jacobian and the estimated variance matrixces corresponding to the reported predictive margins. Taking the squared values of the diagonal elements holds the deltha methods standard errors reported in margins table:

mat J=r(Jacobian)
mat rV=r(V)

di sqrt(rV[1,1])
.00118818

Then I tried to replicate the results above using the following code (adjusted from the one suggested by Jeff in his post):

predict p, mu
gen dpdxb = p*(1-p)
gen dpdx = dpdxb*_b[age]
sum dpdx
gen d2pdxb2 = p*(1-p)*(1-p) - p*p*(1-p)
predict rint, reffects
matrix vecaccum Jac = d2pdxb2 age urban rint
matrix Jac = Jac*_b[age]/1933
qui sum dpdxb
matrix Jac[1,1] = Jac[1,1] + r(mean)
matrix V = Jac*e(V)*Jac'

Summarizing dpdx, I should obtain the ame. However, as you can notice below they are different:

. sum dpdx

    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
        dpdx |      1,934    .0020262    .0002195   .0012703   .0022563

In addition, the standard errors, computed using the matrix V are different:

display sqrt(rV[1,1])
.02592107

I cannot understand what is going on. While I suspect a mistake in the computation of the s.e., I cannot understand why the magnitude of the ame are different. Any suggestion?

Thanks in advance,

Dario