3 level model: PQL1 and MCMC very differents at level 2

[likestatistic]

Dear All,

I am running a 3 level model using Runmlwin,
My response variable is binomial. the PQL1 worked fine but not the PQL2 (I have corrected the issue regarding the PQL2 in runmlwin ado file) , I tried directly the MCMC and increased the chain (100 000) based on the results of the mcmcsum (attached). Could you please explain why the results of the PQL1 and the MCMC are so different particularly at level 2 (0.3255509 vs. 0.004321) ? If it is due to the ESS, How could I increase it at level 2? The MCSE of posterior mean are very low (from 0.0006), is that means that my MCMC results are valid despite the big difference with the PQL1?

Many thanks for your help

Cheers

L

Code: Select all

xi: quietly runmlwin response cons, level3(area: cons) level2(fam: cons)  level1(indiv:) discrete(distribution(binomial) link(logit) denominator(cons)) forcesort nopause mlwinpath(C:\Program Files (x86)\MLwiN v2.32\i386\MLwiN.exe)

xi: quietly runmlwin response cons, level3(area: cons) level2(fam: cons)  level1(indiv:) discrete(distribution(binomial) link(logit) denominator(cons) pql1) initsprevious forcesort nopause mlwinpath(C:\Program Files (x86)\MLwiN v2.32\i386\MLwiN.exe)


MLwiN 2.32 multilevel model                     Number of obs      =    449196
Binomial logit response model
Estimation algorithm: IGLS, PQL1

-----------------------------------------------------------
                |   No. of       Observations per Group
 Level Variable |   Groups    Minimum    Average    Maximum
----------------+------------------------------------------
                area |     6504          3       69.1        525
                fam  |   382701          1        1.2          7
-----------------------------------------------------------

Run time (seconds)   =     252.03
Number of iterations =         19
------------------------------------------------------------------------------
   bookb425w |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        cons |    3.15876   .0101478   311.28   0.000     3.138871    3.178649
------------------------------------------------------------------------------

------------------------------------------------------------------------------
   Random-effects Parameters |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
Level 3: area                     |
                          var(cons) |   .2533781   .0111685      .2314883    .2752679
-----------------------------+------------------------------------------------
Level 2: fam                   |
                       var(cons) |      .3255509   .0442894      .2387452    .4123565
------------------------------------------------------------------------------

After increasing the chain

Code: Select all

xi: runmlwin response cons, level3(area: cons) level2(fam: cons) level1(indiv:) discrete(distribution(binomial) link(logit) denominator(cons)) mcmc(burnin(10000) chain(100000)) initsprevious forcesort nopause mlwinpath(C:\Program Files (x86)\MLwiN v2.32\i386\MLwiN.exe)



MLwiN 2.32 multilevel model                     Number of obs      =    449196
Binomial logit response model
Estimation algorithm: MCMC

-----------------------------------------------------------
                |   No. of       Observations per Group
 Level Variable |   Groups    Minimum    Average    Maximum
----------------+------------------------------------------
    newdatazone |     6504          3       69.1        525
                 fam |   382701          1        1.2          7
-----------------------------------------------------------

Burnin                     =      10000
Chain                      =     100000
Thinning                   =          1
Run time (seconds)         =      46877
Deviance (dbar)            =  152855.48
Deviance (thetabar)        =  150206.94
Effective no. of pars (pd) =    2648.54
Bayesian DIC               =  155504.03
------------------------------------------------------------------------------
   bookb425w |      Mean    Std. Dev.     ESS     P       [95% Cred. Interval]
-------------+----------------------------------------------------------------
        cons |   3.235088   .0112402     5912   0.000      3.21332    3.257249
------------------------------------------------------------------------------

------------------------------------------------------------------------------
   Random-effects Parameters |     Mean   Std. Dev.   ESS     [95% Cred. Int]
-----------------------------+------------------------------------------------
Level 3: area                     |
                          var(cons) |  .2893535  .0123243   3598   .2658652  .3141105
-----------------------------+------------------------------------------------
Level 2: fam                   |
                       var(cons) |   .004321  .0015074     50   .0018132   .007235
------------------------------------------------------------------------------

[GeorgeLeckie]

Hi Likestatistic,

Please note that this question is more about the interpretation of MLwiN output for different models (albeit presented in Stata via runmlwin), rather than on difficulties in actually calling MLwiN from Stata via runmlwin. We try to answer the latter, but have very limited resources to answer the former.

Nonetheless, here are a couple of thoughts...

Note, normally when you move from MQL1/MQL2/PQL1/PQL2 to MCMC, the estimated variance components increase (see MLwiN manual or pervious posts for discussion why). Here you see the opposite pattern for the level-2 variance, which is somewhat unusual.

You have a binary response and very small clusters. In particular your level-2 clusters on average only have 1.2 level-1 units. Computation is always going to be challenging in these circumstances whatever software and estimation method you use. Basically, it is hard to identify the level-2 variance component.

You could try fitting the model with melogit, or two-level versions of this model using xtlogit to try to work out what is going on.

Ultimately, you may have simplify the model to a two-level individuals-within-areas model to proceed with your analysis

Sorry not to be of more help.

Best wishes

George

[likestatistic]

Hi George,

Many thanks for your response. I will apply your recommendations.

Best whises

L