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")
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)
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)
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.