[R-sig-ME] MCMCglmm evaluates distributions differently for different orderings of (unordered) response variables

Christoph Terwitte christerwi at gmail.com
Sun Jun 2 18:44:23 CEST 2013


Dear list members,

I am evaluating subject choices between three possible responses in three conditions. 
Each of 10 subject made the choice in six variations in the three conditions. I therefore want to include subject-offsets in my model, as well as variation-offsets (variations are repeated between conditions). I also factor out condition-specific subject preferences, even though this has little effect on the model outcomes. 
I have had some success modeling the outcome with MCMCglmm, but, to my surprise, the outcome of the MCMCglmm model depends on the ordering of the outcome factor (unordered  in my data frame). below are the summaries of two (converged) MCMCglmm models that I expected to give me the same results: note that the calls are the same, as are the data frames.

Any insights why this happens, if there are more appropriate packages for my purposes, or how to better interpret my results would be greatly appreciated.
Thank you, list!

Christoph
ps: since I am new to the list, i do not know if I can attach code or data frame files… any advice there, too, would be appreciated!
pps: these are my data, script, and summaries:

data frame summary:
> str(christerwi.data)
'data.frame':	154 obs. of  4 variables:
 $ subject  : Factor w/ 10 levels "CC","ChDe","CiDe",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ theme    : Factor w/ 12 levels "sonnenaufgang",..: 12 12 12 11 11 11 10 10 10 9 ...
 $ condition: Factor w/ 3 levels "base","test1",..: 2 1 3 2 1 3 2 1 3 2 ...
 $ choice   : Factor w/ 3 levels "a","b","c": 1 3 1 2 2 1 2 1 1 2 …

> head(christerwi.data,4)
   subject     theme condition choice
5       CC bauernhof     test1      a
9       CC bauernhof      base      c
14      CC bauernhof     test2      a
19      CC    buffet     test1      b

> with(christerwi.data, table(condition, choice))
         choice
condition  a  b  c
    base  24 16 13
    test1 30 10  6
    test2 26 24  5

# for prior specification
k <- 3
I <- diag(k-1)
J <- matrix(rep(1, (k-1)^2), c(k-1, k-1))

# here is the model specification:
mcmcglmm.christerwi<-MCMCglmm(fixed=choice~condition, random = ~idh(condition):subject+subject+theme,data = christerwi.data,rcov = ~us(trait):units,
                              family = "categorical",nitt = 200000,burnin = 5000,thin=150,singular.ok=TRUE, verbose=FALSE,
                              prior=list(R = list(fix=1,V = (1/3)*(I+J)),
                                         G = list(G1 = list(V = diag(3)*0.2,nu = 1),
                                                  G2 = list(V = 0.2, nu = 1), 
                                                  G3 = list(V = 0.2, nu = 1))))

autocorrelation(mcmcglmm.christerwi$Sol) shows that the autocorrelation between successive draws with this combination of 'nitt', 'burnin', and 'thin' is below 0.1. :
, , (Intercept)

          (Intercept) conditiontest1 conditiontest2
Lag 0     1.000000000   -0.610678249    -0.65585856
Lag 100   0.030298125   -0.031276267    -0.01295414 (truncated)

# the same is true of all VCV draws.

# and now the summary indicates that the distribution differences are not significant between the various levels:
summary(mcmcglmm.christerwi)
…
 Location effects: choice ~ condition 

               post.mean l-95% CI u-95% CI eff.samp  pMCMC  
(Intercept)     -0.57857 -1.67346  0.43472     1527 0.2290  
conditiontest1  -1.03532 -2.33849  0.02690     1456 0.0634 .
conditiontest2  -0.05335 -1.08068  1.10419     1333 0.8841  

# NOW BUILD MODEL ON REORDERED OUTCOME VARIABLE:
> christerwi.data$choice<-relevel(christerwi.data$choice, ref='c')
> with(christerwi.data, table(condition, choice))
         choice
condition  c  a  b
    base  13 24 16
    test1  6 30 10
    test2  5 26 24

# now running the same model as above on this data:
mcmcglmm.christerwi.c<-MCMCglmm(fixed=choice~condition, random = ~idh(condition):subject+subject+theme,data = christerwi.data,rcov = ~us(trait):units,
                                family = "categorical",nitt = 300000,burnin = 10000,thin=200,singular.ok=TRUE, verbose=FALSE,
                                prior=list(R = list(fix=1,V = (1/3)*(I+J)),
                                           G = list(G1 = list(V = diag(3)*0.2,nu = 1),
                                                    G2 = list(V = 0.2, nu = 1), 
                                                    G3 = list(V = 0.2, nu = 1))))


# after again verifying the low autocorrelation between successive draws, I check the summary and find this:
summary(mcmcglmm.christerwi.c)
 Location effects: choice ~ condition 

               post.mean l-95% CI u-95% CI eff.samp pMCMC  
(Intercept)      1.03562 -0.56164  2.78721   1047.5 0.183  
conditiontest1   1.47854 -1.01625  3.61806    696.1 0.139  
conditiontest2   1.75178 -0.08553  3.60983   1322.5 0.051 .


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