# How to setup and interpret ANOVA contrasts with the car package in R?

Let's say I have a simple 2x2 factorial experiment that I want to do ANOVA on. Like this, for example:

d   <- data.frame(a=factor(sample(c('a1','a2'), 100, rep=T)),
b=factor(sample(c('b1','b2'), 100, rep=T)));
d$y <- as.numeric(d$a)*rnorm(100, mean=.75, sd=1) +
as.numeric(d$b)*rnorm(100, mean=1.2, sd=1) + as.numeric(d$a)*as.numeric(d$b)*rnorm(100, mean=.5, sd=1) + rnorm(100);  1. In the absence of a significant interaction, by default (i.e. contr.treatment) the output of Anova() is the overall significance of a over all levels of b and of b over all levels of a, is that right? 2. How should I specify a contrast that would allow me to test the significance of effect a with b being held constant at level b1, of effect a with b being held constant at level b2, and of the interaction a:b? - Thanks, @chl, you're right, the new title is better. – f1r3br4nd Jul 13 '11 at 16:21 add comment ## 1 Answer Your example leads to unequal cell sizes, which means that the different "types of sum of squares" matter, and the test for main effects is not as simple as you state it. Anova() uses type II sum of squares. See this question for a start. There are different ways to test the contrasts. Note that SS types don't matter as we are ultimately testing in the associated one-factorial design. I suggest using the following steps: # turn your 2x2 design into the corresponding 4x1 design using interaction() > d$ab <- interaction(d$a, d$b)       # creates new factor coding the 2*2 conditions
> levels(d$ab) # this is the order of the 4 conditions [1] "a1.b1" "a2.b1" "a1.b2" "a2.b2" > aovRes <- aov(y ~ ab, data=d) # oneway ANOVA using aov() with new factor # specify the contrasts you want to test as a matrix (see above for order of cells) > cntrMat <- rbind("contr 01"=c(1, -1, 0, 0), # coefficients for testing a within b1 + "contr 02"=c(0, 0, 1, -1), # coefficients for testing a within b2 + "contr 03"=c(1, -1, -1, 1)) # coefficients for interaction # test contrasts without adjusting alpha, two-sided hypotheses > library(multcomp) # for glht() > summary(glht(aovRes, linfct=mcp(ab=cntrMat), alternative="two.sided"), + test=adjusted("none")) Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: User-defined Contrasts Fit: aov(formula = y ~ ab, data = d) Linear Hypotheses: Estimate Std. Error t value Pr(>|t|) contr 01 == 0 -0.7704 0.7875 -0.978 0.330 contr 02 == 0 -1.0463 0.9067 -1.154 0.251 contr 03 == 0 0.2759 1.2009 0.230 0.819 (Adjusted p values reported -- none method)  Now manually check the result for the first contrast. > P <- 2 # number of levels factor a > Q <- 2 # number of levels factor b > Njk <- table(d$ab)                   # cell sizes
> Mjk     <- tapply(d$y, d$ab, mean)       # cell means
> dfSSE   <- sum(Njk) - P*Q                # degrees of freedom error SS
> SSE     <- sum((d$y - ave(d$y, d\$ab, FUN=mean))^2)    # error SS
> MSE     <- SSE / dfSSE                   # mean error SS
> (psiHat <- sum(cntrMat[1, ] * Mjk))      # contrast estimate
[1] -0.7703638

> lenSq <- sum(cntrMat[1, ]^2 / Njk)       # squared length of contrast
> (SE   <- sqrt(lenSq*MSE))                # standard error
[1] 0.7874602

> (tStat <- psiHat / SE)                   # t-statistic
[1] -0.9782893

> (pVal <- 2 * (1-pt(abs(tStat), dfSSE)))  # p-value
[1] 0.3303902

-
THANK YOU!!! You've just answered a question that two semesters of graduate level stats have not. I even considered using one-way anova before, but couldn't find any confirmation that this was a legitimate approach. –  f1r3br4nd Jul 13 '11 at 18:19
@f1r3br4nd It's legitimate since the error MS is equal in the associated one-way and the original two-way design. –  caracal Jul 13 '11 at 18:33
One last follow-up question, if I may: how does the two-way interaction generalize to interactions of larger numbers of variables? If I had an A*B*C term, would I build that up out of A:B = (A|B=1 - A|B=2), C:B = (C|B=1 - C|B=2), A:B:C = A:B - C:B, and so on? –  f1r3br4nd Jul 13 '11 at 19:53
@f1r3br4nd In a 2x2x2 design, there is only one unique A*B*C interaction contrast (like there is only one in a 2x2 case). The coefficients in an A*B*C interaction-contrast have to sum to zero over rows (A), columns (B), and planes (C) in the "design cube". If the order of cells in the associated one-way design is a1.b1.c1, a2.b1.c1, a1.b2.c1, a2.b2.c1, a1.b1.c2, a2.b1.c2, a1.b2.c2, a2.b2.c2, then the coefficients are c(1, -1, -1, 1, -1, 1, 1, -1). If you have more than two groups in your factors, all contrasts following the sum-to-zero rule are 3-way interaction contrasts. –  caracal Jul 14 '11 at 9:28