Skip to content

MCMCglmm evaluates distributions differently for different orderings of (unordered) response variables

2 messages · Christoph Terwitte, Jarrod Hadfield

Original
Christoph Terwitte
# 
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:
'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 ?
subject     theme condition choice
5       CC bauernhof     test1      a
9       CC bauernhof      base      c
14      CC bauernhof     test2      a
19      CC    buffet     test1      b
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:
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 .
Jarrod Hadfield
# 
Hi,

You should expect different answers depending on which category you  
have as base-line using the model you have. The model has two  
responses (trait): trait 1 is log(Pr(c1))-log(Pr(c3)) and trait 2 is  
log(Pr(c2))-log(Pr(c3)). These are the log-odds ratio of being c1  
versus c3 and the log-odds ratio of being c2 versus c3, where c3 is  
the base-line category. These two responses are indexed by the  
categorical variable `trait' in MCMCglmm.

Currently, you have ~condition as a fixed effect, which means that the  
intercept and the effect of condition is the same for the two  
responses. Depending on the choices you probably want something along  
the lines of ~trait:condition-1 so that separate effects are fitted  
for each.

~theme assumes that the between-theme variances for the two log-odds  
ratios are the same and the correlation between them is one. More  
likely you want or us(trait):theme (they have different variances and  
the correlation is to be estimated).

Only under the fully parameterised model (everything interacted with  
trait in the fixed effects, and us(trait): for the random effects)  
should you expect the models to be equivalent and not depend on the  
choice of base-line category. Even then they will be  
reparameterisations of each other and you will have to do some   
post-analysis manipulation to get the same set of numbers.

Cheers,

Jarrod







Quoting Christoph Terwitte <christerwi at gmail.com> on Sun, 2 Jun 2013  
18:44:23 +0200: