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

[m.fenati at libero.it]

thanks at all,
I use MCMCglmm because in more cases of my analysis I work with perfect 
separation due to zero frequence in several cells of my contingengy tables. 
With MCMCglmm this problem seem to be partially avoided. then I have tried to 
compare glm and MCMCglmm with a dataset where perfect separation is not 
present.

I will try with other prioors psecification.

Cheers


Massimo

>----Messaggio originale----
>Da: bonamy at horus.ens.fr
>Data: 08/12/2011 9.24
>A: <r-sig-mixed-models at r-project.org>
>Ogg: Re: [R-sig-ME] priors with MCMCglmm binomial analysis
>
>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
>
>_______________________________________________
>R-sig-mixed-models at r-project.org mailing list
>https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>