Skip to content

priors with MCMCglmm binomial analysis

4 messages · m.fenati at libero.it, Eryn McFarlane, Jakub Szymkowiak +1 more

Original
m.fenati at libero.it
# 
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:
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")
.........

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
Eryn McFarlane
# 
Hi Massimo,

I'm not sure, but I was under the impression that binomial analysis in MCMCglmm should be done with the family ="categorical". I don't know if this will make a huge difference, but hopefully it helps!

Good Luck,

Eryn

Eryn McFarlane
MSc Candidate
Department of Integrative Biology
University of Guelph
Guelph, ON, Canada N1G 2W1
On 2011-12-07, at 11:14 AM, m.fenati at libero.it wrote:

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Jakub Szymkowiak
    

    

      
      #
    

  


  

      
Hi,
did You try to run chain with more/less number of iteration?

Cheers,
Kuba


-----Oryginalna wiadomo??----- 
From: Eryn McFarlane
Sent: Wednesday, December 07, 2011 5:19 PM
To: m.fenati at libero.it
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] priors with MCMCglmm binomial analysis

Hi Massimo,

I'm not sure, but I was under the impression that binomial analysis in 
MCMCglmm should be done with the family ="categorical". I don't know if this 
will make a huge difference, but hopefully it helps!

Good Luck,

Eryn

Eryn McFarlane
MSc Candidate
Department of Integrative Biology
University of Guelph
Guelph, ON, Canada N1G 2W1

        
On 2011-12-07, at 11:14 AM, m.fenati at libero.it wrote:

        

            

      
      
      
    

        

      
      
    

        
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

  
  


  


      

    

    

      
      

        


  

        

      Pierre 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 :