[R-sig-ME] priors with MCMCglmm binomial analysis

[Pierre B. de Villemereuil]

Dear Massimo,

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 :
prior<-list(R=list(V=0.5,nu=6))

With the following (your burn-in is far too conservative and thin can be 
decreased a bit) :
m1<-MCMCglmm(cbind(pos,(tot-pos))~specie,prior=prior,data=tco,nitt=100000,thin=10,burnin=1000,family="multinomial2")

summary(m1)

  Iterations = 1001:99991
  Thinning interval  = 10
  Sample size  = 9900

  DIC: 401.7165

  R-structure:  ~units

       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 !

Cheers,
Pierre.


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
> results.
>
> 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])
>
> prior<-list(R=list(V=1,nu=0.002))
>
> m.1<-MCMCglmm(cbind(pos,(tot-pos))~specie,prior=prior,data=tco,nitt=900000,
> thin=100,burnin=300000,family="multinomial2",verbose=FALSE)
>
>
> Large coefficient values and intervals occurred after fitting the model:
>
>> summary(m.1)
>   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")
>
>> summary(m.1)
> .........
>
> Coefficients:
>              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!
>
> Regards
>
>
>
> Massimo
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models