Calculating average marginal effects after melogit

[Dario Maimone Ansaldo Patti]

Dear All,

I am really getting crazy about the following point. Suppose I have the following model:

Prob[Y=1|X]=b0+b1*Age+zeta_i+error

I estimate it using a multilevel logit model:

webuse bangladesh, clear
melogit c_use age || district:

Then I calculate the marginal impact of age:

margins, dydx(age) predict(mu)

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

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

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         age |   .0018956   .0011887     1.59   0.111    -.0004342    .0042253
------------------------------------------------------------------------------

I want to replicate the above result manually:

predict p, mu 

su p


    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
           p |      1,934    .3903534    .1014132   .1646984   .6360608

I understood (possibly erroneously) that I could get the same value of p if I use:

mat B= e(b)
mat score double xb = B /* This should refer to the fixed part only*/
predict rint, reffects
gen prob=invlogit(xb+rint)

su prob

  Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
        prob |      1,934    .3903534    .1014132   .1646984   .6360608

This is in fact the case. Then, I would like to calculate the average marginal effects of age:

gen dpdx=prob*(1-prob)
gen me_age=dpdx*_b[age]

su me_age

 Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
      me_age |      1,934    .0019474    .0001967   .0011766   .0021382

Surprisingly, the result is slightly different from the one obtained using margins:

margins, dydx(age) predict(mu)

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

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

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         age |   .0018956   .0011887     1.59   0.111    -.0004342    .0042253
------------------------------------------------------------------------------

Why this is the case? Is there any mistake in my procedure?

Thanks in advance for your insights.

Best

Dario