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),
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]]