Post hoc comparison in two way ANOVA with covariate using R

Question

Does any body know how to run a post hoc comparison in a 2X2 ANOVA with covariate in R. multcomp package seems very nice but I could not find clear answer or example to my question with this package. Thanks so much in advance

Answer 1

This is data from Maxwell & Delaney (2004), artificially extended to include a second between-subjects IV, yielding a 3x2 design. Using multcomp's glht() function is easier once you switch to the associated one-factorial design by combining your two IVs into one with interaction().

DV is depression scores pre-treatment and post-treatment. Treatment is one of SSRI, Placebo or Waiting List. Pre-treatment score is the covariate. I included a second IV.

P        <- 3     # number of groups in IV1
Q        <- 2     # number of groups in IV2
Njk      <- 5     # cell size
SSRIpre  <- c(18, 16, 16, 15, 14, 20, 14, 21, 25, 11)
SSRIpost <- c(12,  0, 10,  9,  0, 11,  2,  4, 15, 10)
PlacPre  <- c(18, 16, 15, 14, 20, 25, 11, 25, 11, 22)
PlacPost <- c(11,  4, 19, 15,  3, 14, 10, 16, 10, 20)
WLpre    <- c(15, 19, 10, 29, 24, 15,  9, 18, 22, 13)
WLpost   <- c(17, 25, 10, 22, 23, 10,  2, 10, 14,  7)
IV1      <- factor(rep(1:3,  each=Njk*Q), labels=c("SSRI", "Placebo", "WL"))
IV2      <- factor(rep(1:2, times=Njk*P), labels=c("A", "B"))

# combine both IVs into 1 to get the associated one-factorial design
IVi      <- interaction(IV1, IV2)
DVpre    <- c(SSRIpre,  PlacPre,  WLpre)
DVpost   <- c(SSRIpost, PlacPost, WLpost)

Now do the ANCOVA with IVi as between-subjects factor and DVpre as covariate. Using the associated one-factorial design is possible since it has the same Error-MS as the two-factorial design.

> aovAncova1 <- aov(DVpost ~ IVi     + DVpre, data=dfAncova)  # one-factorial design
> aovAncova2 <- aov(DVpost ~ IV1*IV2 + DVpre, data=dfAncova)  # two-factorial design
> summary(aovAncova1)[[1]][["Mean Sq"]][3]  # Error MS one-factorial design
[1] 32.84399

> summary(aovAncova2)[[1]][["Mean Sq"]][5]  # Error MS two-factorial design
[1] 32.84399

Next comes the matrix defining 3 cell comparisons with sum-to-zero coefficients. The coefficients follow the order of 2*3 levels of IVi.

> levels(IVi)
[1] "SSRI.A" "Placebo.A" "WL.A" "SSRI.B" "Placebo.B" "WL.B"

> cntrMat <- rbind("SSRI-Placebo"  = c(-1, 1, 0, -1,  1,  0),
+                  "SSRI-0.5(P+WL)"= c(-2, 1, 1, -2,  1,  1),
+                  "A-B"           = c( 1, 1, 1, -1, -1, -1))

> library(multcomp)
> summary(glht(aovAncova1, linfct=mcp(IVi=cntrMat), alternative="greater"),
+         test=adjusted("none"))

Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: User-defined Contrasts
Fit: aov(formula = DVpost ~ IVi + DVpre, data = dfAncova)

Linear Hypotheses:
                    Estimate Std. Error t value  Pr(>t)   
SSRI-Placebo <= 0      8.887      5.135   1.731 0.04846 * 
SSRI-WL <= 0          12.878      5.129   2.511 0.00976 **
A-B <= 0               1.024      6.492   0.158 0.4380
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
(Adjusted p values reported -- none method)