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
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
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
------------------------------------------------------------------------------
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