[R-sig-ME] priors with MCMCglmm binomial analysis
Pierre B. de Villemereuil
bonamy at horus.ens.fr
Thu Dec 8 09:24:04 CET 2011
I don't understand why you are using MCMCglmm since you don't have any
random effect... Do you plan to add some in the future ? For now, the
only difference between glm and MCMCglmm is that MCMCglmm fit
overdispersion by estimating 'units' (which is set to zero in glm models).
I'm not an expert in real binomial data (by that I mean binomial but not
binary), so I'm not sure about the following, but wouldn't it be
possible that only 3 levels of observations (even if it totalizes 341
observations) is too small to get rid of prior sensitivity ? Especially
since your variance is an overdispersion variance.
I say so because by the look of you credible interval, you seem to get
an inverse-Gamma posterior distribution for 'units' ! Indeed, if you use
an informative prior, such as :
With the following (your burn-in is far too conservative and thin can be
decreased a bit) :
Iterations = 1001:99991
Thinning interval = 10
Sample size = 9900
post.mean l-95% CI u-95% CI eff.samp
units 0.7527 0.1359 1.863 9812
Location effects: cbind(pos, (tot - pos)) ~ specie
post.mean l-95% CI u-95% CI eff.samp pMCMC
(Intercept) -1.786720 -3.549124 0.002936 9900 0.0491 *
speciecap 0.774466 -1.828374 3.425216 9900 0.5123
specieov 1.659190 -0.924082 4.063666 11072 0.1608
There is no high auto-correlation and convergence actually happened (you
should try, it is far quicker that way !).
Anyway, my advice would be : if you don't need MCMCglmm (e.g. random
effects), then don't use it. If you still need it for some reason, try
to fix the residual variance to some value. I think 0 isn't possible due
to mixing necessity, but you could check, otherwise, maybe 1 ?
Good luck !
Le 07/12/2011 17:14, m.fenati at libero.it a écrit :
> dear R users,
> I don't know how to set a MCMCglmm for binomial analysis. I'm trying to fit
> several model but bad results were always observed if compared with the glm
> My model is:
> tc<-matrix(c(124,184,33,18,86,9), ncol=2)
> tco<-data.frame("specie"=c("bov","ov","cap"),"pos"= tc[,2],"tot"=tc[,1])
> Large coefficient values and intervals occurred after fitting the model:
> Iterations = 300001:899901
> Thinning interval = 100
> Sample size = 6000
> DIC: 401.758
> R-structure: ~units
> post.mean l-95% CI u-95% CI eff.samp
> units 594395401 0.0002638 1.418e+09 5106
> Location effects: cbind(pos, (tot - pos)) ~ specie
> post.mean l-95% CI u-95% CI eff.samp pMCMC
> (Intercept) 59.24 -20274.73 22702.05 6000 0.697
> speciecap -62.28 -33022.96 28154.63 6000 0.815
> specieov -142.82 -28811.83 30635.06 6000 0.717
> When I fit a glm model, different and more reasonable results are observed:
> m.1<-glm(cbind(pos,(tot-pos))~specie, data=tco, family="binomial")
> Estimate Std. Error z value Pr(>|z|)
> (Intercept) -1.7731 0.2549 -6.955 3.52e-12 ***
> speciecap 0.7922 0.4667 1.698 0.0896 .
> specieov 1.6424 0.2947 5.574 2.49e-08 ***
> Where is the error? a mistake in the model specification, priors setting, ...
> I dont'know!
> I tried with other priors but the results never agree with the glm results.
> Could someone help me?
> Thanks in advance!
> R-sig-mixed-models at r-project.org mailing list
More information about the R-sig-mixed-models