[R-sig-ME] testing simple main effects by contrasts

[Erich Studerus]

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
>