[R-sig-ME] MCMCglmm and contrast

m.fenati at libero.it m.fenati at libero.it
Thu Sep 22 20:02:32 CEST 2011


Dear R users,
I would like to have a suggestion about MCMCglmm and contr.sdif.
I want to test the difference between two succesive level of a categorical 
variable (data_cat) in a multivariate model.
I fit the model as follows:

contrasts(cimu$data_cat)=contr.sdif

prior<-list(R=list(V=1,fix=1),G=list(G1=list(V=1,nu=0.002)),B=list(mu=c(rep
(0,11)),V=diag(11)*(3+pi^2/3)))
m.1<-MCMCglmm(edv_tot~age+sex+data_cat,slice=T,prior=prior,random=~ID_an,
data=cimu,nitt=900000,thin=100,burnin=300000,family="categorical",
verbose=FALSE)

summary(m.1)

Sample size  = 6000 

 DIC: 35.83431 

 G-structure:  ~ID_an

      post.mean  l-95% CI u-95% CI eff.samp
ID_an     3.214 0.0002216    13.79    520.3

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units         1        1        1        0

 Location effects: edv_tot ~ age + sex + data_cat 

            post.mean l-95% CI u-95% CI eff.samp  pMCMC   
(Intercept)   -0.9443  -3.3426   1.4717     4861 0.4360   
age1          -2.2578  -6.2168   1.4007     5584 0.2390   
age2           2.2226  -0.6441   5.3706     5628 0.1270   
age3           0.5132  -2.4640   3.7798     4989 0.7343   
age4           0.5891  -2.8786   4.0271     6000 0.7270   
age5          -0.3487  -3.1947   2.4558     5610 0.8213   
age6          -0.3023  -3.3455   2.7215     5834 0.8553   
sexM          -0.7528  -3.0999   1.6399     5588 0.5220   
data_cat2-1    2.3829  -0.5235   5.5751     4317 0.1220   
data_cat3-2   -3.3894  -6.3156  -0.1769     6000 0.0323 * 
data_cat4-3   -4.6328  -7.8367  -1.6581     4217 0.0020 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 


I aspect to obtain for "data_cat2-1" the same coefficient that I obtain if I 
use contr.treatment , but it is not the case (see the following example)

contrasts(cimu$data_cat)=contr.treatment

prior<-list(R=list(V=1,fix=1),G=list(G1=list(V=1,nu=0.002)),B=list(mu=c(rep
(0,11)),V=diag(11)*(3+pi^2/3)))
m.1<-MCMCglmm(edv_tot~age+sex+data_cat,slice=T,prior=prior,random=~ID_an,
data=cimu,nitt=900000,thin=100,burnin=300000,family="categorical",
verbose=FALSE)


 Iterations = 300001:899901
 Thinning interval  = 100
 Sample size  = 6000 

 DIC: 34.31511 

 G-structure:  ~ID_an

      post.mean  l-95% CI u-95% CI eff.samp
ID_an     2.134 0.0002254    10.21    905.5

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units         1        1        1        0

 Location effects: edv_tot ~ age + sex + data_cat 

            post.mean l-95% CI u-95% CI eff.samp    pMCMC    
(Intercept)   -0.6479  -3.2589   1.9114     6000 0.624333    
age1          -1.9388  -5.6190   1.9704     5216 0.314000    
age2           2.2851  -0.7129   5.2292     5347 0.131667    
age3           0.5602  -2.8026   3.4921     6000 0.710000    
age4           0.5284  -2.9907   3.9128     6000 0.756667    
age5           0.0143  -2.5979   2.7103     5635 0.969000    
age6          -0.3660  -3.6664   2.7862     6000 0.847333    
sexM          -0.6188  -2.8901   1.6339     6000 0.582667    
data_cat2      3.9755   0.6346   7.4141     6000 0.016000 *  
data_cat3     -0.1361  -2.5624   2.4428     5570 0.894667    
data_cat4     -4.8936  -7.8206  -1.9578     5026 0.000333 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 


Can someone help me?

Thank you in advance

Massimo




More information about the R-sig-mixed-models mailing list