multilevel logistic model with MCMCglmm

Dear list members,

I fit a random intercept morel with lme4 first, then I tried to fit the same
model with MCMCglmm. I didn't specify the prior and ran 300000 iterations.
The model still didn't converge as the autocorrelation is very high (>0.4).
May I ask (1) do I have to specify my own prior? and (2) what is the
usually/average number of iterations in Bayesian multilevel model? 300000
seems quite large a number. Do I have to rerun the model if the previous one
doesn't converge?

Best with thanks.

data(Contraception,package="mlmRev")
library(lme4)
fm1 <-

glmer(use~urban+age+livch+(1|district),data=Contraception,family=binomial,nAGQ=7)

MC1 <-

MCMCglmm(use~urban+age+livch,random=~district,data=Contraception,family="categorical",nitt=300000,burnin=240000,thin=100)

summary(MC1)

 Iterations = 240001:299901
 Thinning interval  = 100
 Sample size  = 600

 DIC: 47.28245

 G-structure:  ~district

         post.mean l-95% CI u-95% CI eff.samp
district      3242      953     6490      132

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units     38025    13640    57500    85.53

 Location effects: use ~ urban + age + livch

            post.mean l-95% CI u-95% CI eff.samp  pMCMC
(Intercept)  -197.901 -262.707 -130.966    92.68 <0.002 **
urbanY         86.448   53.374  119.831   131.85 <0.002 **
age            -3.139   -5.439   -1.200   166.36 <0.002 **
livch1        128.027   76.737  177.583   126.78 <0.002 **
livch2        159.752  102.289  222.468   109.28 <0.002 **
livch3+       155.254   92.145  216.916   123.80 <0.002 **
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

autocorr(MC1$Sol)

<snip>

--
Wincent Ronggui HUANG
Sociology Department of Fudan University
PhD of City University of Hong Kong
http://asrr.r-forge.r-project.org/rghuang.html

	[[alternative HTML version deleted]]

Hi Professor Hadfield, I see.

In the multilevel logistic model, the level-1 error variance is fixed
pi^2/3, which I perhaps should follow.

Thank you very much

BTW, the MCMCglmm is a great package.

Regards,
Ronggui

On 3 April 2011 23:53, Jarrod Hadfield <j.hadfield at ed.ac.uk> wrote:

Hi Wincent,

The model you ran is meaningless, because with categorical data or ordinal
data the residual variance cannot be estimated from the data. The best
option is to fix the residual variance at some value (I use 1):
  
MC1<-MCMCglmm(use~urban+age+livch,random=~district,data=Contraception,family="categorical",prior=list(R=list(V=1,
fix=1), G=list(G1=list(V=1, nu=0))), slice=T)

I have left the prior on the district variance flat and improper, and I
have used slice sampling (slice=T) instead of Metropolis-Hastings updates
because the chain mixes faster:

summary(MC1)

 Iterations = 3001:12991
 Thinning interval  = 10
 Sample size  = 1000

 DIC: 2391.208


 G-structure:  ~district

        post.mean l-95% CI u-95% CI eff.samp
district    0.3333   0.1407   0.5777    459.9


 R-structure:  ~units

     post.mean l-95% CI u-95% CI eff.samp
units         1        1        1        0


 Location effects: use ~ urban + age + livch

           post.mean l-95% CI u-95% CI eff.samp  pMCMC
(Intercept)  -2.02718 -2.38427 -1.69340    559.7 <0.001 ***
urbanY        0.87874  0.61785  1.16036    640.9 <0.001 ***
age          -0.03233 -0.05013 -0.01523    531.4 <0.001 ***
livch1        1.32248  0.96341  1.66027    578.6 <0.001 ***
livch2        1.63986  1.24274  2.05330    542.6 <0.001 ***
livch3+       1.61311  1.18382  2.00119    512.2 <0.001 ***

Note that you can rescale the estimates to approximate what they would be
if the residual variance was zero (i.e. the glmer default):

c2<-((16*sqrt(3))/(15*pi))^2

posterior.mode(MC1$Sol/sqrt(1+c2))
(Intercept)      urbanY         age      livch1      livch2     livch3+
-1.69345255  0.72259159 -0.03015856  1.14554895  1.38585993  1.34449154

which are closer to the ML estimates.

Cheers,

Jarrod



 Quoting Wincent <ronggui.huang at gmail.com> on Sun, 3 Apr 2011 20:32:27
+0800:

 Dear list members,

I fit a random intercept morel with lme4 first, then I tried to fit the
same
model with MCMCglmm. I didn't specify the prior and ran 300000 iterations.
The model still didn't converge as the autocorrelation is very high
(>0.4).
May I ask (1) do I have to specify my own prior? and (2) what is the
usually/average number of iterations in Bayesian multilevel model? 300000
seems quite large a number. Do I have to rerun the model if the previous
one
doesn't converge?

Best with thanks.

 data(Contraception,package="mlmRev")

library(lme4)
fm1 <-

glmer(use~urban+age+livch+(1|district),data=Contraception,family=binomial,nAGQ=7)

 MC1 <-

MCMCglmm(use~urban+age+livch,random=~district,data=Contraception,family="categorical",nitt=300000,burnin=240000,thin=100)

summary(MC1)

 Iterations = 240001:299901
 Thinning interval  = 100
 Sample size  = 600

 DIC: 47.28245

 G-structure:  ~district

        post.mean l-95% CI u-95% CI eff.samp
district      3242      953     6490      132

 R-structure:  ~units

     post.mean l-95% CI u-95% CI eff.samp
units     38025    13640    57500    85.53

 Location effects: use ~ urban + age + livch

           post.mean l-95% CI u-95% CI eff.samp  pMCMC
(Intercept)  -197.901 -262.707 -130.966    92.68 <0.002 **
urbanY         86.448   53.374  119.831   131.85 <0.002 **
age            -3.139   -5.439   -1.200   166.36 <0.002 **
livch1        128.027   76.737  177.583   126.78 <0.002 **
livch2        159.752  102.289  222.468   109.28 <0.002 **
livch3+       155.254   92.145  216.916   123.80 <0.002 **
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

 autocorr(MC1$Sol)

<snip>

--
Wincent Ronggui HUANG
Sociology Department of Fudan University
PhD of City University of Hong Kong
http://asrr.r-forge.r-project.org/rghuang.html

       [[alternative HTML version deleted]]

--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.

--
Wincent Ronggui HUANG
Sociology Department of Fudan University
PhD of City University of Hong Kong
http://asrr.r-forge.r-project.org/rghuang.html