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

[Reinhold Kliegl]

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
>>
>
>
>