negative binomial multilevel model: different results in stata and mlwin

[pimverbunt]

Hello all,
I am modelling a multilevel negative binomial model. To check the consistency of my results, I ran the same model (same variables, same dataset) in both mlwin and stata. I was surprised to find out that the random intercepts and overdispersion factor drastically differs in stata (command: menbreg) and mlwin. For instance, in mlwin the random intercepts equals 0.75(0.187) and the overdispersion factor equals 1.269(0.009). In stata, the random intercepts and overdispersion factor equal 0.124(0.002) and 1.44, respectively.
Both programs do give very similar results when modelling a linear or logistic regression model.
What could be the cause of the discrepency in results?

[ChrisCharlton]

The closest MLwiN estimation method would be IGLS with PQL2 turned on. When I run a negative-binomial model using one of the MLwiN example data sets using the -menbreg- and -runmlwin- I get the following results:

Code: Select all

. use http://www.bristol.ac.uk/cmm/media/runmlwin/mmmec, clear
. generate lnexpected = ln(exp)

Code: Select all

. menbreg obs uvbi, offset(lnexpected) || nation: || region:

Fitting fixed-effects model:

Iteration 0:   log likelihood =  -1361.855  
Iteration 1:   log likelihood = -1230.0211  
Iteration 2:   log likelihood = -1211.0491  
Iteration 3:   log likelihood = -1202.5642  
Iteration 4:   log likelihood = -1202.5329  
Iteration 5:   log likelihood = -1202.5329  

Refining starting values:

Grid node 0:   log likelihood = -1209.6951

Fitting full model:

Iteration 0:   log likelihood = -1209.6951  (not concave)
Iteration 1:   log likelihood = -1195.0761  (not concave)
Iteration 2:   log likelihood = -1189.7235  (not concave)
Iteration 3:   log likelihood = -1167.5801  (not concave)
Iteration 4:   log likelihood = -1145.4326  (not concave)
Iteration 5:   log likelihood = -1138.4473  
Iteration 6:   log likelihood = -1088.3883  
Iteration 7:   log likelihood = -1086.7993  
Iteration 8:   log likelihood = -1086.4087  
Iteration 9:   log likelihood = -1086.3905  
Iteration 10:  log likelihood = -1086.3903  
Iteration 11:  log likelihood = -1086.3903  

Mixed-effects nbinomial regression              Number of obs     =        354
Overdispersion:            mean

-------------------------------------------------------------
                |     No. of       Observations per Group
 Group Variable |     Groups    Minimum    Average    Maximum
----------------+--------------------------------------------
         nation |          9          3       39.3         95
         region |         78          1        4.5         13
-------------------------------------------------------------

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(1)      =       8.73
Log likelihood = -1086.3903                     Prob > chi2       =     0.0031
-------------------------------------------------------------------------------
          obs |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------+----------------------------------------------------------------
         uvbi |  -.0335933   .0113725    -2.95   0.003     -.055883   -.0113035
        _cons |  -.0790615   .1295931    -0.61   0.542    -.3330592    .1749362
   lnexpected |          1  (offset)
--------------+----------------------------------------------------------------
     /lnalpha |  -4.182598   .3415021   -12.25   0.000     -4.85193   -3.513266
--------------+----------------------------------------------------------------
nation        |
    var(_cons)|   .1283613    .067897                      .0455186    .3619756
--------------+----------------------------------------------------------------
nation>region |
    var(_cons)|   .0401818   .0104856                      .0240938    .0670121
-------------------------------------------------------------------------------
LR test vs. nbinomial model: chi2(2) = 232.29             Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.

. // Convert overdispersion to same scale as reported by MLwiN
. display exp(_b[lnalpha:_cons])
.01525882

Code: Select all

. runmlwin obs cons uvbi, /// 
> level3(nation: cons) level2(region: cons) level1(county: ) ///
> discrete(distribution(nbinomial) pql2 offset(lnexpected)) nopause
 
MLwiN 2.36 multilevel model                     Number of obs      =       354
Negative binomial response model
Estimation algorithm: IGLS, PQL2

-----------------------------------------------------------
                |   No. of       Observations per Group
 Level Variable |   Groups    Minimum    Average    Maximum
----------------+------------------------------------------
         nation |        9          3       39.3         95
         region |       78          1        4.5         13
-----------------------------------------------------------

Run time (seconds)   =       1.66
Number of iterations =          6
------------------------------------------------------------------------------
         obs |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        cons |  -.0975399   .1293803    -0.75   0.451    -.3511205    .1560408
        uvbi |  -.0334264   .0112881    -2.96   0.003    -.0555508   -.0113021
------------------------------------------------------------------------------

------------------------------------------------------------------------------
   Random-effects Parameters |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
Level 3: nation              |
                   var(cons) |   .1272512   .0679558     -.0059397    .2604422
-----------------------------+------------------------------------------------
Level 2: region              |
                   var(cons) |   .0403602   .0103289      .0201159    .0606044
------------------------------------------------------------------------------
------------------------------------------------------------------------------
   Overdispersion Parameters |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
                    bcons2_1 |   .0152484   .0047656       .005908    .0245889
------------------------------------------------------------------------------

These results look fairly close. If you use these MLwiN estimation settings with your model do you see the same discrepancies that you reported before?

[pimverbunt]

Thanks Chris for your response. I have been selecting different integration methods in stata and it appears that Mode and curvature adaptive Gauss-Hermite quadrature (command: intmethod(mcaghermite)) and laplacian approximation (command: intmethod(laplace)) give results that are similar to IGLS/PQL2 in MlWin.