I have two questions concerning planned contrasts:
- I would like to know how factor-based contrasts (obtained through an interaction term) compare to model-paramter-based contrasts (obtained by specifying model parameters).
- I would like to know this for a simple case of comparing one condition with another, but ultimately I am interested in comparing one condition vs all other conditions.
Below are my attempts at understanding factor-based contrasts and model-parameter based contrasts:
#create some dummy data
data <- mtcars
#create interaction terms
data$interaction <- interaction(mtcars$am, mtcars$vs, sep="X")
data$interaction <- gsub("^0", "am0", data$interaction)
data$interaction <- gsub("^1", "am1", data$interaction)
data$interaction <- gsub("0$", "vs0", data$interaction)
data$interaction <- gsub("1$", "vs1", data$interaction)
data$interaction <- factor(data$interaction)
levels(data$interaction)
#[1] "am0Xvs0" "am0Xvs1" "am1Xvs0" "am1Xvs1"
From Eric Fuchs' blogpost I think I figured out how to obtain factor-based contrasts. Let us assume for now that I am interested in the comparison of am0Xvs0 vs. am0Xvs1. To this purpose, I create a contrast matrix where I assign equal weights with opposing signs to my two levels of interest, and 0 to the other two levels:
#specify contrasts:
c.f <- c(-1, 1, 0, 0)
mat.f <- cbind(c.f)
contrasts(data$interaction) <- mat.f
#fit model
fit.f <- aov(mpg~interaction, data)
#get coefficients for contrasts
summary(fit.f, split=list(interaction=list("am0Xvs0 vs. am0Xvs1"=1)))
Output:
Df Sum Sq Mean Sq F value Pr(>F)
interaction 3 788.6 262.86 21.81 1.73e-07 ***
interaction: am0Xvs0 vs. am0Xvs1 1 232.3 232.28 19.27 0.000147 ***
Residuals 28 337.5 12.05
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
For comparison, my attempt at model-parameter-based contrasts:
fit.m <- aov(mpg~am*vs, data)
with $E[mpg]=b_0 + b_1 \text{am} + b_2 \text{vs} + b3 (\text{am} \times \text{vs}).$ Based on Matt Blackwell's answer I think that a comparison of am0Xvs0 vs. am0Xvs1 means that $H_0: b_1 = 0$ (i.e., the regression weight associated with am). Therefore:
## construct contrast matrices
mat.m <- rbind("am0:vs0 - vs1" = c(0, 0, 1, 0))
library(car)
lht(fit.m, mat.m)
Output:
Linear hypothesis test
Hypothesis:
am = 0
Model 1: restricted model
Model 2: mpg ~ am * vs
Res.Df RSS Df Sum of Sq F Pr(>F)
1 29 425.84
2 28 337.48 1 88.36 7.3311 0.01142 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The values differ, so one of the two solutions is not correct (and my suspicion is that I did not correctly translate Matt's answer to this case). My question now is: which one of the two approaches is correct and how can the other be rewritten correctly?
Ultimately, I am interested in the comparison of 1 of the 4 levels of the interaction term vs. the other 3. Let's say I want to compare am0Xvs0 vs. the other 3 conditions. The factor-based version is simply an extension of what I have written above (if what I wrote above was correct):
#create contrast matrix
c.f2 <- c(1, -1/3, -1/3, -1/3)
mat.f2 <- cbind(c.f2)
contrasts(data$interaction) <- mat.f2
#fit model
fit.f2 <- aov(mpg~interaction, data)
#get coefficients for contrasts
summary(fit.f2, split=list(interaction=list("am0Xvs0 vs. rest"=1)))
Output:
Df Sum Sq Mean Sq F value Pr(>F)
interaction 3 788.6 262.9 21.81 1.73e-07 ***
interaction: am0Xvs0 vs. rest 1 487.8 487.8 40.48 6.95e-07 ***
Residuals 28 337.5 12.1
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Unfortunately I would not know how to start with the parameters for the model-parameter-based version, which is what I am ultimately interested in.
