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

Jarrod Hadfield j.hadfield at ed.ac.uk
Mon Jun 3 17:12:46 CEST 2013


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:

> 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 .
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>


-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.



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