Hi,
I have the following study design: Each of two groups of subjects (group A
and group B) was tested repeatedly on two experimental days. At each
experimental day, placebo or drug were administered in a randomized order
and subjects were tested under the the influence of the drug in 3
consecutive blocks. I have therefore one between subjects factor (Group A vs
Group B) and two within subjects factors (Drug and Block). Here's a
reproducible example:
x <- expand.grid(Block= paste('B', 1:3, sep=''),
Drug = c('Placebo', 'Drug'), Subj=1:24)
df <- data.frame(x, Group = rep(c('A','B'), each=72),
value=rnorm(144, sd=10))
m <- model.matrix(lm(value~Group*Drug+Block,df))
df$value <- df$value+m[,2]*7+ m[,3]*2 + m[,4]*5 + m[,5]*8 + m[,6]-2
#delete some observations to simulate an unbalanced design
df <- df[-as.numeric(sample(rownames(df),20)),]
mod <- lme(value ~ Group*Drug+Block, random = ~1|Subj, data = df)
anova(mod)
Now, since the interaction between Drug and Group is significant, the main
effects of Drug and Group can not be interpreted. I therefore define
contrasts to test the simple main effects, that is, the effect of the Drug
in each Group, averaged over all levels of the Block factor.
I use the contrast package to set up the contrasts:
library(contrast)
con1 <- contrast(mod, list(Block=c('B1','B2','B3'), Group = 'A', Drug =
'Pla'),
list(Block=c('B1','B2','B3'), Group = 'A', Drug = 'Drug'),
type='average')
con2 <- contrast(mod, list(Block=c('B1','B2','B3'), Group = 'B', Drug =
'Pla'),
list(Block=c('B1','B2','B3'), Group = 'B', Drug = 'Drug'),
type='average')
Now, I use the multcomp-package to correct the p-values for multiple
comparisons:
library(multcomp)
summary(glht(mod, linfct = rbind(con1$X,con2$X)))
This approach works fine with lme, but not with lmer. In the model above I
modeled the block factor as a fixed effect. However, according to my
understanding, it would be more parsimonious to model the block factor as a
random effect. I'm not directly interested in its effect and it's just an
additional source of variance due to habituation effects . Unfortunately,
the Block and Subj factors are crossed and therefore can not easily be
modeled by lme. After comparing different models, the best fitting model
with lmer looks like this:
lmer(value ~ Group*Drug+(Drug|Subj)+(1|Block))
How can I set up contrasts to test simple main effects for this lmer model?
I know, that lmer has a contrasts argument. However, as far as I can see,
one can only test contrasts for the levels of one factor at a time. Any
comments are highly appreciated.
Erich
testing simple main effects by contrasts
4 messages · Erich Studerus, Reinhold Kliegl
From what I understand you want to test the effect of Drug within each
of the two Groups. Here is a simple way to do this:
# Convert 2 Drug x 2 Group design to 1 x 4 Cond design
df$Cond <- factor(paste(df$Group, df$Drug, sep="_"))
# ANOVA: one main, two nested fixed effects
cmat <- matrix(c( -1, -1, +1, +1, # Main effect of Group
-1, +1, 0, 0, # Effect of Drug | Group==A
0, 0, -1, +1), 4, 3) # Effect of Drug | Group==B
colnames(cmat) <- c(".group", ".drug_A", ".drug_B")
contrasts(df$Cond) <- cmat/2
print(mod <- lmer(value ~ Cond + (1 |Subj), data=df), cor=FALSE)
print(mod <- lmer(value ~ Cond + Block + (1 |Subj), data=df),
cor=FALSE) # Block as fixed effect
print(mod <- lmer(value ~ Cond + (1 |Subj) + (1|Block), data=df),
cor=FALSE) # Block as random effect
The fixed-effect coefficients reported correspond to
(Intercept) = mean of four conditions
Group = difference between the two groups
.drug_A = difference between Drug levels in Group A
.drug_B = difference between Drug levels in Group B
and associated standard errors and "t-values". Estimates line up with
table means for balanced designs.
You mention 2 Days of testing but this factor did not appear in the example.
The models you tested with lme and lmer where structurally different
with respect to Block. Given that you have only 3 Blocks, I would
include Block as a fixed effect, not as a random effect, but opinions
differ on this. Anyway, this has not much of a bearing on the
estimates of the Group difference and the estimates of Drug effects
within Groups.
Reinhold Kliegl
On Sat, Jun 13, 2009 at 4:57 PM, Erich
Studerus<erich.studerus at bli.uzh.ch> wrote:
Hi,
I have the following study design: Each of two groups of subjects (group A
and group B) was tested repeatedly on two experimental days. At each
experimental day, placebo or drug were administered in a randomized order
and subjects were tested under the the influence of the drug in 3
consecutive blocks. I have therefore one between subjects factor (Group A vs
Group B) and two within subjects factors (Drug and Block). Here's a
reproducible example:
x <- expand.grid(Block= paste('B', 1:3, sep=''),
?Drug = c('Placebo', 'Drug'), Subj=1:24)
df <- data.frame(x, Group = rep(c('A','B'), each=72),
?value=rnorm(144, sd=10))
m <- model.matrix(lm(value~Group*Drug+Block,df))
df$value <- df$value+m[,2]*7+ m[,3]*2 + m[,4]*5 + m[,5]*8 + m[,6]-2
#delete some observations to simulate an unbalanced design
df <- df[-as.numeric(sample(rownames(df),20)),]
mod <- lme(value ~ Group*Drug+Block, random = ~1|Subj, data = df)
anova(mod)
Now, since the interaction between Drug and Group is significant, the main
effects of Drug and Group can not be interpreted. I therefore define
contrasts to test the simple main effects, that is, the effect of the Drug
in each Group, averaged over all levels of the Block factor.
I use the contrast package to set up the contrasts:
library(contrast)
con1 <- contrast(mod, list(Block=c('B1','B2','B3'), Group = 'A', Drug =
'Pla'),
? ? list(Block=c('B1','B2','B3'), Group = 'A', Drug = 'Drug'),
type='average')
con2 <- contrast(mod, list(Block=c('B1','B2','B3'), Group = 'B', Drug =
'Pla'),
? ? list(Block=c('B1','B2','B3'), Group = 'B', Drug = 'Drug'),
type='average')
Now, I use the multcomp-package to correct the p-values for multiple
comparisons:
library(multcomp)
summary(glht(mod, linfct = rbind(con1$X,con2$X)))
This approach works fine with lme, but not with lmer. In the model above I
modeled the block factor as a fixed effect. However, according to my
understanding, it would be more parsimonious to model the block factor as a
random effect. I'm not directly interested in its effect and it's just an
additional source of variance due to habituation effects . Unfortunately,
the Block and Subj factors are crossed and therefore can not easily be
modeled by lme. After comparing different models, the best fitting model
with lmer looks like this:
lmer(value ~ Group*Drug+(Drug|Subj)+(1|Block))
How can I set up contrasts to test simple main effects for this lmer model?
I know, that lmer has a contrasts argument. However, as far as I can see,
one can only test contrasts for the levels of one factor at a time. Any
comments are highly appreciated.
Erich
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
1 day later
Thank you so much. I highly appreciate your help.
However, when I extend my example to the case of a design with 3 drug
levels, the results that I get with the contrast-package and the results
that I get with your approach seem to differ. I don't know why this is the
case. Could you please have a look at the following example and tell me what
I'm doing wrong.
#prepare an example data.frame
set.seed(5)
x <- expand.grid(Block= paste('B', 1:3, sep=''),
Drug = c('Pla', 'Drug1','Drug2'), Subj=1:24)
df <- data.frame(x, Group = rep(c('A','B'), each=108),
value=rnorm(216, sd=10))
m <- model.matrix(lm(value~Group*Drug+Block,df))
df$value <- df$value+m[,2]*7+ m[,3]*2 + m[,4]*5+
m[,5]*8 + m[,6]*-2+m[,7]*3+m[,8]*-5
df <- df[-as.numeric(sample(rownames(df),20)),]
mod <- lme(value~Group*Drug+Block, random=~1|Subj,data=df)
library(contrast)
#Effect of Drug1 | Group A
contrast(mod,
list(Block=c('B1','B2','B3'), Group = 'A', Drug = 'Pla'),
list(Block=c('B1','B2','B3'), Group = 'A', Drug = 'Drug1'),
type='average')
#Effect of Drug2 | Group A
contrast(mod,
list(Block=c('B1','B2','B3'), Group = 'A', Drug = 'Pla'),
list(Block=c('B1','B2','B3'), Group = 'A', Drug = 'Drug2'),
type='average')
#Effect of Drug1 | Group B
contrast(mod,
list(Block=c('B1','B2','B3'), Group = 'B', Drug = 'Pla'),
list(Block=c('B1','B2','B3'), Group = 'B', Drug = 'Drug1'),
type='average')
#Effect of Drug2 | Group B
contrast(mod,
list(Block=c('B1','B2','B3'), Group = 'B', Drug = 'Pla'),
list(Block=c('B1','B2','B3'), Group = 'B', Drug = 'Drug2'),
type='average')
df$Cond <- with(df, interaction(Drug, Group))
#Main effect of Group
group <- c(-1,-1,-1,1,1,1)
#Effect of Drug1 | Group A
Drug1.A <- c(-1,1,0,0,0,0)
#Effect of Drug2 | Group A
Drug2.A <- c(-1,0,1,0,0,0)
#Effect of Drug1 | Group B
Drug1.B <- c(0,0,0,-1,1,0)
#Effect of Drug2 | Group B
Drug2.B <- c(0,0,0,-1,0,1)
contrasts(df$Cond) <- cbind(group, Drug1.A, Drug2.A,
Drug1.B, Drug2.B)/2
summary(lme(value ~ Cond + Block, random=~1|Subj,data=df))
-----Urspr?ngliche Nachricht-----
Von: Reinhold Kliegl [mailto:reinhold.kliegl at gmail.com]
Gesendet: Samstag, 13. Juni 2009 20:20
An: Erich Studerus
Cc: r-sig-mixed-models at r-project.org
Betreff: Re: [R-sig-ME] testing simple main effects by contrasts
From what I understand you want to test the effect of Drug within each
of the two Groups. Here is a simple way to do this:
# Convert 2 Drug x 2 Group design to 1 x 4 Cond design
df$Cond <- factor(paste(df$Group, df$Drug, sep="_"))
# ANOVA: one main, two nested fixed effects
cmat <- matrix(c( -1, -1, +1, +1, # Main effect of Group
-1, +1, 0, 0, # Effect of Drug |
Group==A
0, 0, -1, +1), 4, 3) # Effect of Drug |
Group==B
colnames(cmat) <- c(".group", ".drug_A", ".drug_B")
contrasts(df$Cond) <- cmat/2
print(mod <- lmer(value ~ Cond + (1 |Subj), data=df), cor=FALSE)
print(mod <- lmer(value ~ Cond + Block + (1 |Subj), data=df),
cor=FALSE) # Block as fixed effect
print(mod <- lmer(value ~ Cond + (1 |Subj) + (1|Block), data=df),
cor=FALSE) # Block as random effect
The fixed-effect coefficients reported correspond to
(Intercept) = mean of four conditions
Group = difference between the two groups
.drug_A = difference between Drug levels in Group A
.drug_B = difference between Drug levels in Group B
and associated standard errors and "t-values". Estimates line up with
table means for balanced designs.
You mention 2 Days of testing but this factor did not appear in the example.
The models you tested with lme and lmer where structurally different
with respect to Block. Given that you have only 3 Blocks, I would
include Block as a fixed effect, not as a random effect, but opinions
differ on this. Anyway, this has not much of a bearing on the
estimates of the Group difference and the estimates of Drug effects
within Groups.
Reinhold Kliegl
On Sat, Jun 13, 2009 at 4:57 PM, Erich
Studerus<erich.studerus at bli.uzh.ch> wrote:
Hi, I have the following study design: Each of two groups of subjects (group A and group B) was tested repeatedly on two experimental days. At each experimental day, placebo or drug were administered in a randomized order and subjects were tested under the the influence of the drug in 3 consecutive blocks. I have therefore one between subjects factor (Group A
vs
Group B) and two within subjects factors (Drug and Block). Here's a
reproducible example:
x <- expand.grid(Block= paste('B', 1:3, sep=''),
?Drug = c('Placebo', 'Drug'), Subj=1:24)
df <- data.frame(x, Group = rep(c('A','B'), each=72),
?value=rnorm(144, sd=10))
m <- model.matrix(lm(value~Group*Drug+Block,df))
df$value <- df$value+m[,2]*7+ m[,3]*2 + m[,4]*5 + m[,5]*8 + m[,6]-2
#delete some observations to simulate an unbalanced design
df <- df[-as.numeric(sample(rownames(df),20)),]
mod <- lme(value ~ Group*Drug+Block, random = ~1|Subj, data = df)
anova(mod)
Now, since the interaction between Drug and Group is significant, the main
effects of Drug and Group can not be interpreted. I therefore define
contrasts to test the simple main effects, that is, the effect of the Drug
in each Group, averaged over all levels of the Block factor.
I use the contrast package to set up the contrasts:
library(contrast)
con1 <- contrast(mod, list(Block=c('B1','B2','B3'), Group = 'A', Drug =
'Pla'),
? ? list(Block=c('B1','B2','B3'), Group = 'A', Drug = 'Drug'),
type='average')
con2 <- contrast(mod, list(Block=c('B1','B2','B3'), Group = 'B', Drug =
'Pla'),
? ? list(Block=c('B1','B2','B3'), Group = 'B', Drug = 'Drug'),
type='average')
Now, I use the multcomp-package to correct the p-values for multiple
comparisons:
library(multcomp)
summary(glht(mod, linfct = rbind(con1$X,con2$X)))
This approach works fine with lme, but not with lmer. In the model above I
modeled the block factor as a fixed effect. However, according to my
understanding, it would be more parsimonious to model the block factor as
a
random effect. I'm not directly interested in its effect and it's just an additional source of variance due to habituation effects . Unfortunately, the Block and Subj factors are crossed and therefore can not easily be modeled by lme. After comparing different models, the best fitting model with lmer looks like this: lmer(value ~ Group*Drug+(Drug|Subj)+(1|Block)) How can I set up contrasts to test simple main effects for this lmer
model?
I know, that lmer has a contrasts argument. However, as far as I can see, one can only test contrasts for the levels of one factor at a time. Any comments are highly appreciated. Erich
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
You use a "contr.sum" specification for your Drug contrasts, leaving
out the Placebe level. That is, your coefficients represent the
difference of Drug1 and Drug2 from the mean of all three levels of
the Drug factor (within Group). This is probably not what you want.
You probably want to use Placebo as a reference condition for Drug1
and Drug2. If so, then you need a "contr.treatment" specification:
##
# Extension to three drug levels with treatment contrasts within groups
cmat.t <- matrix(c( -1, -1, -1, +1, +1, +1, # Main effect of Group
0, +1, 0, 0, 0, 0, # Effect of Drug.1
vs. Placebo | Group==A
0, 0, +1, 0, 0, 0, # Effect of Drug.2
vs. Placebo | Group==A
0, 0, 0, 0, +1, 0, # Effect of Drug.1
vs. Placebo | Group==B
0, 0, 0, 0, 0, +1), 6, 5) # Effect of Drug.1
vs. Placebo | Group==B
colnames(cmat.t) <- c(".group", ".drug1_A", ".drug2_A",".drug1_B", ".drug2_B")
contrasts(df$Cond) <- cmat.r
print(mod <- lmer(value ~ Cond + Block + (1 |Subj), data=df), cor=FALSE)
## Alternatively, perhaps you want to compare drug1 against placebo
and drug2 against drug1
# Extension to three drug levels with simple difference contrasts within groups
cmat.r <- matrix(c( -1, -1, -1, +1, +1, +1, # Main effect of Group
-2/3, +1/3, +1/3, 0, 0, 0, # Drug.1
vs. Placebo | Group==A
-1/3, -1/3, +2/3, 0, 0, 0, # Drug.2
vs. Drug1 | Group==A
0, 0, 0, -2/3, +1/3, +1/3, # Drug.1
vs. Placebo | Group==B
0, 0, 0, -1/3, -1/3, +2/3), 6, 5) # Drug.2
vs. Drug1 | Group==B
colnames(cmat.r) <- c(".group", ".drug1-plac_A",
".drug2-drug1_A",".drug1-plac_B", ".drug2-plac_B")
contrasts(df$Cond) <- cmat.r
print(mod <- lmer(value ~ Cond + Block + (1 |Subj), data=df), cor=FALSE)
############
I suspect the contrast package uses the treatment contrast, but I am
not familiar with it. Anyway, it does not matter what you use, but it
is important that coefficients estimate what you expect them to
estimate.
Reinhold Kliegl
On Mon, Jun 15, 2009 at 3:04 PM, Erich
Studerus<erich.studerus at bli.uzh.ch> wrote:
Thank you so much. I highly appreciate your help.
However, when I extend my example to the case of a design with 3 drug
levels, the results that I get with the contrast-package and the results
that I get with your approach seem to differ. I don't know why this is the
case. Could you please have a look at the following example and tell me what
I'm doing wrong.
#prepare an example data.frame
set.seed(5)
x <- expand.grid(Block= paste('B', 1:3, sep=''),
Drug = c('Pla', 'Drug1','Drug2'), Subj=1:24)
df <- data.frame(x, Group = rep(c('A','B'), each=108),
value=rnorm(216, sd=10))
m <- model.matrix(lm(value~Group*Drug+Block,df))
df$value <- df$value+m[,2]*7+ m[,3]*2 + m[,4]*5+
?m[,5]*8 + m[,6]*-2+m[,7]*3+m[,8]*-5
df <- df[-as.numeric(sample(rownames(df),20)),]
mod <- lme(value~Group*Drug+Block, random=~1|Subj,data=df)
library(contrast)
#Effect of Drug1 | Group A
contrast(mod,
list(Block=c('B1','B2','B3'), Group = 'A', Drug = 'Pla'),
list(Block=c('B1','B2','B3'), Group = 'A', Drug = 'Drug1'),
type='average')
#Effect of Drug2 | Group A
contrast(mod,
list(Block=c('B1','B2','B3'), Group = 'A', Drug = 'Pla'),
list(Block=c('B1','B2','B3'), Group = 'A', Drug = 'Drug2'),
type='average')
#Effect of Drug1 | Group B
contrast(mod,
list(Block=c('B1','B2','B3'), Group = 'B', Drug = 'Pla'),
list(Block=c('B1','B2','B3'), Group = 'B', Drug = 'Drug1'),
type='average')
#Effect of Drug2 | Group B
contrast(mod,
list(Block=c('B1','B2','B3'), Group = 'B', Drug = 'Pla'),
list(Block=c('B1','B2','B3'), Group = 'B', Drug = 'Drug2'),
type='average')
df$Cond <- with(df, interaction(Drug, Group))
#Main effect of Group
group <- c(-1,-1,-1,1,1,1)
#Effect of Drug1 | Group A
Drug1.A <- c(-1,1,0,0,0,0)
#Effect of Drug2 | Group A
Drug2.A <- c(-1,0,1,0,0,0)
#Effect of Drug1 | Group B
Drug1.B <- c(0,0,0,-1,1,0)
#Effect of Drug2 | Group B
Drug2.B <- c(0,0,0,-1,0,1)
contrasts(df$Cond) <- cbind(group, Drug1.A, Drug2.A,
Drug1.B, Drug2.B)/2
summary(lme(value ~ Cond + Block, random=~1|Subj,data=df))
-----Urspr?ngliche Nachricht-----
Von: Reinhold Kliegl [mailto:reinhold.kliegl at gmail.com]
Gesendet: Samstag, 13. Juni 2009 20:20
An: Erich Studerus
Cc: r-sig-mixed-models at r-project.org
Betreff: Re: [R-sig-ME] testing simple main effects by contrasts
From what I understand you want to test the effect of Drug within each
of the two Groups. Here is a simple way to do this:
# Convert 2 Drug x 2 Group design to 1 x 4 Cond design
df$Cond <- factor(paste(df$Group, df$Drug, sep="_"))
# ANOVA: one main, two nested fixed effects
cmat <- matrix(c( -1, -1, +1, +1, ? ? ? ? ? ? # Main effect of Group
? ? ? ? ? ? ? ? ? ? ? ? ? ? -1, +1, ?0, ?0, ? ? ? ? ? ? # Effect of Drug |
Group==A
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 0, ?0, -1, +1), ?4, ?3) ? ?# Effect of Drug |
Group==B
colnames(cmat) <- c(".group", ".drug_A", ".drug_B")
contrasts(df$Cond) <- cmat/2
print(mod <- lmer(value ~ Cond + (1 |Subj), data=df), cor=FALSE)
print(mod <- lmer(value ~ Cond + Block + (1 |Subj), data=df),
cor=FALSE) ?# Block as fixed effect
print(mod <- lmer(value ~ Cond + (1 |Subj) + (1|Block), data=df),
cor=FALSE) # Block as random effect
The fixed-effect coefficients reported correspond to
(Intercept) ?= mean of four conditions
Group ?= difference between the two groups
.drug_A = difference between Drug levels in Group A
.drug_B = difference between Drug levels in Group B
and associated standard errors and "t-values". ?Estimates line up with
table means for balanced designs.
You mention 2 Days of testing but this factor did not appear in the example.
The models you tested with lme and lmer where structurally different
with respect to Block. Given that you have only 3 ?Blocks, I would
include Block as a fixed effect, not as a random effect, but opinions
differ on this. Anyway, this has not much of a bearing on the
estimates of the Group difference and the estimates of Drug effects
within Groups.
Reinhold Kliegl
On Sat, Jun 13, 2009 at 4:57 PM, Erich
Studerus<erich.studerus at bli.uzh.ch> wrote:
Hi, I have the following study design: Each of two groups of subjects (group A and group B) was tested repeatedly on two experimental days. At each experimental day, placebo or drug were administered in a randomized order and subjects were tested under the the influence of the drug in 3 consecutive blocks. I have therefore one between subjects factor (Group A
vs
Group B) and two within subjects factors (Drug and Block). Here's a
reproducible example:
x <- expand.grid(Block= paste('B', 1:3, sep=''),
?Drug = c('Placebo', 'Drug'), Subj=1:24)
df <- data.frame(x, Group = rep(c('A','B'), each=72),
?value=rnorm(144, sd=10))
m <- model.matrix(lm(value~Group*Drug+Block,df))
df$value <- df$value+m[,2]*7+ m[,3]*2 + m[,4]*5 + m[,5]*8 + m[,6]-2
#delete some observations to simulate an unbalanced design
df <- df[-as.numeric(sample(rownames(df),20)),]
mod <- lme(value ~ Group*Drug+Block, random = ~1|Subj, data = df)
anova(mod)
Now, since the interaction between Drug and Group is significant, the main
effects of Drug and Group can not be interpreted. I therefore define
contrasts to test the simple main effects, that is, the effect of the Drug
in each Group, averaged over all levels of the Block factor.
I use the contrast package to set up the contrasts:
library(contrast)
con1 <- contrast(mod, list(Block=c('B1','B2','B3'), Group = 'A', Drug =
'Pla'),
? ? list(Block=c('B1','B2','B3'), Group = 'A', Drug = 'Drug'),
type='average')
con2 <- contrast(mod, list(Block=c('B1','B2','B3'), Group = 'B', Drug =
'Pla'),
? ? list(Block=c('B1','B2','B3'), Group = 'B', Drug = 'Drug'),
type='average')
Now, I use the multcomp-package to correct the p-values for multiple
comparisons:
library(multcomp)
summary(glht(mod, linfct = rbind(con1$X,con2$X)))
This approach works fine with lme, but not with lmer. In the model above I
modeled the block factor as a fixed effect. However, according to my
understanding, it would be more parsimonious to model the block factor as
a
random effect. I'm not directly interested in its effect and it's just an additional source of variance due to habituation effects . Unfortunately, the Block and Subj factors are crossed and therefore can not easily be modeled by lme. After comparing different models, the best fitting model with lmer looks like this: lmer(value ~ Group*Drug+(Drug|Subj)+(1|Block)) How can I set up contrasts to test simple main effects for this lmer
model?
I know, that lmer has a contrasts argument. However, as far as I can see, one can only test contrasts for the levels of one factor at a time. Any comments are highly appreciated. Erich
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models