www.cmm.bristol.ac.uk/forum

Posted: **Mon Feb 01, 2016 2:11 pm**

[edit: Fixed the title. Cleaned up the Stata code.]

Hello,

I've been playing around with simulated 2-level count data. MQL1 estimates of the intercept and variance component were quite different from the true parameter values (xtmepoisson estimates were much better), so I suspected that was due to the estimation technique and moved to PQL2 instead. However, i received an error message that seems to suggest that some internal computations went wrong. Am I doing something wrong or is this indeed a bug?

I also tried using the MLwiN user interface and got a similar error. MCMC didn't really converge to anything sensible within the first 50,000 iterations.

Code to generate the simulated data:

Code: Select all

.         clear
.         set obs 200     // Number of lvl2 units
.         set seed 4744
.         gen u = rnormal(0,sqrt(1.8))	// Random effect with variance 1.8
.         gen lvl2_id = _n
.         expand 500      // Number of lvl1 units per lvl2 unit - overall sample is 100,000

.         gen lvl1_id = _n
.         gen x1 = invnormal(runiform())
.         gen x2 = invnormal(runiform())
.         gen x3 = invnormal(runiform())
.         gen xb = 2 + 0.75*x1 + 1.2*x2 + (-0.3)*x3 + u
.         gen Y = rpoisson(exp(xb))

.         sort lvl2_id lvl1_id
.         gen cons = 1

xtmepoisson results (look correct):

Code: Select all

.         xtmepoisson Y x1 x2 x3 || lvl2_id:, var intpoints(23)

Mixed-effects Poisson regression                Number of obs     =    100,000
Group variable: lvl2_id                         Number of groups  =        200

                                                Obs per group:
                                                              min =        500
                                                              avg =      500.0
                                                              max =        500

Integration points =  23                        Wald chi2(3)      =   9.15e+06
Log likelihood = -233925.35                     Prob > chi2       =     0.0000

------------------------------------------------------------------------------
           Y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
          x1 |   .7504717   .0004783  1568.96   0.000     .7495342    .7514092
          x2 |   1.199206    .000471  2545.87   0.000     1.198282    1.200129
          x3 |  -.2997022   .0004642  -645.69   0.000    -.3006119   -.2987925
       _cons |   1.884631   .0945808    19.93   0.000     1.699256    2.070006
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
lvl2_id: Identity            |
                  var(_cons) |   1.788779   .1789234      1.470331    2.176199
------------------------------------------------------------------------------
LR test vs. Poisson model: chibar2(01) = 8.6e+06      Prob >= chibar2 = 0.0000

MLwiN with MQL1 (quite far off target):

Code: Select all

.         qui runmlwin Y cons x1 x2 x3, ///
>                 level2(lvl2_id: cons) level1(lvl1_id:) ///
>                 discrete(distribution(poisson) link(log)) ///
>                 nopause batch mlwinpath("C:\Program Files (x86)\MLwiN v2.35\x64\mlnscript.exe") maxiterations(100) 
.         assert e(converged)==1
.         runmlwin, nolog
 
MLwiN 2.34 multilevel model                     Number of obs      =    100000
Poisson response model
Estimation algorithm: IGLS, MQL1

-----------------------------------------------------------
                |   No. of       Observations per Group
 Level Variable |   Groups    Minimum    Average    Maximum
----------------+------------------------------------------
        lvl2_id |      200        500      500.0        500
-----------------------------------------------------------

Run time (seconds)   =       1.16
Number of iterations =          7
------------------------------------------------------------------------------
           Y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        cons |   2.824754   .1524951    18.52   0.000     2.525869    3.123639
          x1 |   .7504708   .0004662  1609.86   0.000     .7495571    .7513844
          x2 |   1.199205   .0004587  2614.08   0.000     1.198305    1.200104
          x3 |  -.2997025   .0004604  -651.02   0.000    -.3006048   -.2988002
------------------------------------------------------------------------------

------------------------------------------------------------------------------
   Random-effects Parameters |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
Level 2: lvl2_id             |
                   var(cons) |   4.650818     .46491      3.739612    5.562025
------------------------------------------------------------------------------

MLwiN with PQL2 (estimation error):

Code: Select all

.         runmlwin Y cons x1 x2 x3, ///
>                 level2(lvl2_id: cons) level1(lvl1_id:) ///
>                 discrete(distribution(poisson) link(log) pql2) initsprevious ///
>                 nopause batch mlwinpath("C:\Program Files (x86)\MLwiN v2.35\x64\mlnscript.exe") maxiterations(100)
 
Model fitted using initial values specified as parameter estimates from previous model

 --- Begin MLwiN error log --- 
MLN - Software for N-level analysis.   02/01/16 13:59:20

C:\Users\ng526\AppData\Local\Temp\ST_03000002.tmp

error while obeying batch file C:\Users\ng526\AppData\Local\Temp\ST_03000008.tmp
>  at line :
calc 'F~(H)' = expo('H')
2500 numeric errors in calculate command, first 20 affected entries listed. 
Affected entries replaced with system missing.

error while obeying batch file C:\Users\ng526\AppData\Local\Temp\ST_03000008.tmp
>  at line :
calc 'P' = expo('H')
2500 numeric errors in calculate command, first 20 affected entries listed. 
Affected entries replaced with system missing.

error while obeying batch file C:\Users\ng526\AppData\Local\Temp\ST_03000008.tmp
>  at line :
calc g9 = g9 * 'P'^0.5
3500 numeric errors in calculate command, first 20 affected entries listed. 
Affected entries replaced with system missing.

C:\Users\ng526\AppData\Local\Temp\ST_03000006.tmp

C:\Users\ng526\AppData\Local\Temp\ST_03000007.tmp
 --- End MLwiN error log --- 

runmlwin has encountered an error importing the model results from MLwiN. Check 
> that the model has run properly in MLwiN.

r(198);

end of do-file

Posted: **Sun Feb 21, 2016 7:58 pm**

A little push...

Posted: **Wed Apr 13, 2016 3:03 pm**

Hi Nils,

Your DGP and runmlwin syntax all look fine. The problem appears to be with some particularly high counts generated by your DGP. Essentially if you specify a lower true value for your level-2 variance everything is fine, but as you pump up the level-2 variance you get more and more extreme counts. Your most extreme count is some 14,000 for whatever your hypothetical response is.

Basically the way MLwiN algorithms for these models go about doing all the internal calculations is more sensitive to working with such extreme counts than xtmepoisson. MLwiN holds data to less precision than Stata. Stata also works on more stable metrics than MLwiN.

I'm afraid, if you want to use MLwiN you will have to consider a more restricted range of true values for your level-2 variance.

Best wishes

George

www.cmm.bristol.ac.uk/forum

2-level Poisson model - MQL1 works fine, PQL2 crashes

2-level Poisson model - MQL1 works fine, PQL2 crashes

Re: 2-level Poisson model - MQL1 works fine, PQL2 crashes

Re: 2-level Poisson model - MQL1 works fine, PQL2 crashes