[R-sig-ME] multiple comparison for a LME with interaction
Kay Cecil Cichini
Kay.Cichini at uibk.ac.at
Thu Aug 4 19:05:23 CEST 2011
hi again,
here's an example - if you compare the coefficents between the models,
or between the models and the plot, as well as summing up the single
coeffecients in the custom model, you will see that this all makes
good sense!
if anyone of the "experts" would comment on this approach i'd be happy, too..
### code:
ChickWeight_new <- ChickWeight
ChickWeight_new$Time_factor <- as.factor(ifelse(ChickWeight_new$Time >
0, 1, 0))
str(ChickWeight_new)
with(ChickWeight_new, interaction.plot(Time_factor, Diet, weight, col = 1:4))
library(nlme)
library(multcomp)
mod1 <- lme(weight ~ Diet*Time_factor, random = ~ 1|Chick, data =
ChickWeight_new)
ChickWeight_new$Diet <- with(ChickWeight_new, relevel(Diet, ref = "2"))
mod2 <- lme(weight ~ Diet * Time_factor, random = ~ 1|Chick, data =
ChickWeight_new)
ChickWeight_new$Diet <- with(ChickWeight_new, relevel(Diet, ref = "3"))
mod3 <- lme(weight ~ Diet * Time_factor, random = ~ 1|Chick, data =
ChickWeight_new)
ChickWeight_new$Diet <- with(ChickWeight_new, relevel(Diet, ref = "4"))
mod4 <- lme(weight ~ Diet * Time_factor, random = ~ 1|Chick, data =
ChickWeight_new)
c1 <- rbind("Time effect within Diet1" = c(0,0,0,0,1,0,0,0),
"Time effect within Diet2" = c(0,0,0,0,1,1,0,0),
"Time effect within Diet3" = c(0,0,0,0,1,0,1,0),
"Time effect within Diet4" = c(0,0,0,0,1,0,0,1),
"Time Effect: Diet1 vs. Diet2" = c(0,0,0,0,0,1,0,0),
"Time Effect: Diet1 vs. Diet3" = c(0,0,0,0,0,0,1,0),
"Time Effect: Diet1 vs. Diet4" = c(0,0,0,0,0,0,0,1),
"Time Effect: Diet2 vs. Diet3" = c(0,0,0,0,0,-1,1,0),
"Time Effect: Diet2 vs. Diet4" = c(0,0,0,0,0,-1,0,1),
"Time Effect: Diet3 vs. Diet4" = c(0,0,0,0,0,0,-1,1))
summary(glht(mod1, c1))
summary(mod1)
summary(mod2)
summary(mod3)
summary(mod4)
###
Zitat von Kay Cecil Cichini <Kay.Cichini at uibk.ac.at>:
> hi,
>
> you will need to provide example data, otherwise it will be hard for
> anyone to offer help.
>
> in any case you could try if the contrasts are set up properly by
> releveling your factors (relevel-function) using the default
> treatment contrasts, then compare with your custom contrasts and see
> if the coefficients are the same between the two models where they
> should be...
>
> interactionplot() is also nice for checking if your coefficients make sense.
>
> best,
> kay
>
>
>
> Zitat von Poulin Nicolas <poulin at math.unistra.fr>:
>
>> Dear all,
>> I would like to know if you think that what I am doing is correct. I
>> have always the same self-questionning when working with (generalized)
>> linear mixed models with interactions.
>>
>> First let me introduce you my data.
>> A group of monkeys have been experimented within 3 different
>> conditions (variable Nb with values 1, 2 and 3).
>> As a social group of monkeys, the rank among the group of each monkey is
>> important in the way of explaining a reaction (variable Rk).
>> I want to see the effect of the condition, the rank and the interaction
>> of these two variables on the outcome of the experiment (variable Hz).
>>
>> There are 540 observations among 15 individuals (36 per individual)
>>
>> The selected model :
>> - the variable Nb have been transformed into a factor
>> - the square root of the variable Hz was used to be able to use a LME
>> instead of a GLMM (valided by a Q-Q plot)
>> - Hz seemed to be lineary dependent of the rank so I didn't used a
>> factor version of Rk (hence Rk is numeric)
>>
>> I only tried the random intercept based on the identity (ID) of each
>> individual.
>>
>> So the model I am using is :
>> model.sqrt<-lme(sqrtHz~Rk*Nb, random = ~ 1 | ID,data=data1)
>>
>> Here is the summary of the fixed part of the model :
>> Fixed effects: sqrtHz ~ Rk * Nb
>>
>> Value Std.Error DF
>> t-value p-value
>> (Intercept) 0.5031499 0.04741670 521 10.611238 0.0000
>> Rk -0.0456437 0.00641704 13 -7.112890 0.0000
>> Nb2 0.0551886 0.02609983 521 2.114520 0.0349
>> Nb3 0.0026095 0.02609983 521 0.099982 0.9204
>> Rk:Nb2 -0.0004107 0.00353217 521 -0.116277 0.9075
>> Rk:Nb3 0.0144365 0.00353217 521 4.087143 0.0001
>>
>> To obtain a p-value for the factor Nb and the intercation, I use :
>> > anova(model.sqrt)
>> numDF denDF F-value p-value
>> (Intercept) 1 521 113.85619 <.0001
>> Rk 1 13 45.33842 <.0001
>> Nb 2 521 18.49245 <.0001
>> Rk:Nb 2 521 11.46233 <.0001
>>
>>
>> Hence both Nb and the interaction have significant effects.
>>
>> My first question is do I have the right to do such a thing?
>>
>> I want to compare the different experiment and the interactions.
>>
>> First, I tried to use summary(glht(model.sqrt,linfct=mcp(Nb="Tukey")))
>> that give me the result :
>>
>> Estimate Std. Error z value Pr(>|z|)
>> 2 - 1 == 0 0.05519 0.02610 2.115 0.0868 .
>> 3 - 1 == 0 0.00261 0.02610 0.100 0.9945
>> 3 - 2 == 0 -0.05258 0.02610 -2.015 0.1088
>>
>> This output came with the warning message : covariate interactions found
>> -- default contrast might be inappropriate
>>
>> Moreover I can not apply this to have multiple comparisons among the
>> interactions.
>>
>> Browsing the mailing list, I tried to follow the method given by :
>> http://thebiobucket.blogspot.com/2011/06/glmm-with-custom-multiple-comparisons.html
>>
>> Hence I did :
>> > c <- rbind("Condition 1 - Condition 2" = c(0,0,1,0,0,0),
>> + "Condition 1 - Condition 3" = c(0,0,0,1,0,0),
>> + "Condition 2 - Condition 3" = c(0,0,1,-1,0,0),
>> + "Rank:Condition 1 - Rank:Condition 2" = c(0,0,0,0,1,0),
>> + "Rank:Condition 1 - Rank:Condition 2" = c(0,0,0,0,0,1),
>> + "Rank:Condition 1 - Rank:Condition 2" = c(0,0,0,0,1,-1)
>> + )
>>
>> and
>>
>> > summary(glht(model.sqrt,c))
>>
>> Simultaneous Tests for General Linear Hypotheses
>>
>> Fit: lme.formula(fixed = sqrtHz ~ Rk * Nb, data = data1, random = ~1 |
>> ID)
>>
>> Linear Hypotheses:
>>
>> Estimate Std. Error z value Pr(>|z|)
>> Condition 1 - Condition 2 == 0 0.0551886
>> 0.0260998 2.115 0.131
>> Condition 1 - Condition 3 == 0 0.0026095
>> 0.0260998 0.100 1.000
>> Condition 2 - Condition 3 == 0 0.0525791
>> 0.0260998 2.015 0.162
>> Rank:Condition 1 - Rank:Condition 2 == 0 -0.0004107 0.0035322
>> -0.116 1.000
>> Rank:Condition 1 - Rank:Condition 2 == 0 0.0144365 0.0035322
>> 4.087 <0.001 ***
>> Rank:Condition 1 - Rank:Condition 2 == 0 -0.0148472 0.0035322
>> -4.203 <0.001 ***
>>
>> So my real question is : is this correct?
>> If so, it's quite in contradiction with the results of
>> summary(model.sqrt) and anova(model.sqrt), for the multiple comparisons
>> among the different conditions. I know that the p-values of
>> summary(model.sqrt) can be seen as Dunnett comparison with Condition 1.
>>
>> Here I don't know what is the correction that is used but I think it
>> takes into consideration the fact that I have done 6 different tests.
>> These corrections seems to be too conservative as there is no
>> significant differences among the different conditions.
>> Is there a way to select a less conservative correction?
>>
>> If not I don't know how to explain the differences between the results
>> of both summary(model.sqrt) and anova(model.sqrt) and the last
>> displayed one.
>>
>> Thanks in advance for your answers, and appologies to have written such
>> a long message.
>> Best,
>> Nicolas
>> --
>> Nicolas Poulin
>> Ingénieur de Recherche
>> Centre de Statistique de Strasbourg (CeStatS)
>> http://www-math.u-strasbg.fr/CeStatS/
>> <http://www-math.u-strasbg.fr/CentreStat/>
>> IRMA, UMR 7501
>> Université de Strasbourg et CNRS
>> 7 rue René-Descartes
>> 67084 Strasbourg Cedex
>>
>> Tél. 33 (0)3 68 85 01 89
>>
>> [[alternative HTML version deleted]]
>>
>>
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
More information about the R-sig-mixed-models
mailing list