Contrasts in repeated mesures ANOVA in R

Question

I am a refugee from SPSS in the process of re-learning how to do everything in R. Mostly it's been fun, as R is great for a lot of things but I've run into a snag that I can't seem to find a solution for; any help appreciated.

I want to run repeated-measures ANOVA with contrasts. I know linear mixed models are generally the better way to go but I work in a community where that is only just starting to catch on and my students need to understand what they're reading when they encounter older techniques. SPSS makes contrasts on repeated measures easy, but R does weird things I don't understand. Here's a sample:

df.wide = data.frame(Subj = 1:8, Day1 = c(6,5,5,6,7,4,4,5), 
                     Day2 = c(5,5,6,5,3,2,4,7), Day3 = c(2,4,3,4,3,1,1,2))
attach(df.wide)

## contrasts as t-tests
Linear = -Day1 + 0*Day2 + Day3
Day1vs23 = -2*Day1 + Day2 + Day3
t.test(Linear)
t.test(Day1vs23)

## repeated-measures with long-format data
library(reshape2)
df.long = melt(df.wide, id.vars="Subj", variable.name="Trial", value.name="DV")
contrasts(df.long$Trial) = contr.poly(3)
df.long$Subj = factor(df.long$Subj)
Anv = aov(DV ~ Trial + Error(Subj/Trial), data=df.long)
summary(Anv, split=list(Trial=list("Linear"=1, "Quad"=2)))

The t-tests show t-values of 7.51 for Linear and 1.27 for quadratic, but this does NOT match the aov() output, which provides F-values of 25.28 and 2.51, and of course substantially different p-values. SPSS contrasts on the repeated measures yield t-values and F-values that match up.

It looks like R isn't partitioning the error term like SPSS does for contrasts, and both contrasts are being tested against a 14 d.f. error term, which is why it doesn't match the t-test. That seems wrong to me.

So, what am I doing wrong in R and how do I fix it?

Answer 1

Your hunch about what is going on is on the right track. While the "t-test approach" tests each contrast against its own separate error term, the RM-ANOVA tests each contrast against a pooled error term (hence DF=14), which assumes that the variances of the contrast scores are approximately equal (i.e., sphericity assumption). The results from these two methods are noticeably different here because, as we will see, the variances of the contrasts are very clearly unequal, so the analysis that doesn't rely on sphericity is a better choice here. If you are dead-set on pooling the error terms, for example in order to test the 2-DF Trial effect, you can use one of the corrections that we will see below.

(Before moving on to the rest of the explanation, I first want to quickly point out only the Linear contrast is comparable between your two analyses, as the Day1v23 contrast and Quad contrasts that you defined are clearly not equivalent.)

First let's replicate the t-test approach using a multivariate linear model (mlm):

### assume that all objects from your example code are loaded already
(conts <- as.matrix(df.wide[,-1]) %*% contrasts(df.long$Trial))
#              .L            .Q
# [1,] -2.8284271 -8.164966e-01
# [2,] -0.7071068 -4.082483e-01
# [3,] -1.4142136 -1.632993e+00
# [4,] -1.4142136  2.220446e-16
# [5,] -2.8284271  1.632993e+00
# [6,] -2.1213203  4.082483e-01
# [7,] -2.1213203 -1.224745e+00
# [8,] -2.1213203 -2.857738e+00
mod <- lm(conts ~ 1)
summary(mod)
# Response .L :
#   
#   Call:
#   lm(formula = .L ~ 1)
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -0.8839 -0.3535 -0.1768  0.5303  1.2374 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  -1.9445     0.2588  -7.514 0.000136 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
# 
# Residual standard error: 0.7319 on 7 degrees of freedom
# 
# 
# Response .Q :
#   
#   Call:
#   lm(formula = .Q ~ 1)
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -2.2454 -0.7144  0.0000  0.7144  2.2454 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)  -0.6124     0.4818  -1.271    0.244
# 
# Residual standard error: 1.363 on 7 degrees of freedom

Now let's see what it looks like when explicitly using the pooled error term approach, just to demonstrate its equivalence to your second analysis:

(C <- colSums(conts^2))
# .L .Q 
# 34 16 
### these are the variances of the reduced or "small" models

(A <- colSums(scale(conts, scale=F)^2))
#   .L    .Q 
# 3.75 13.00 
### these are the variances of the full or "large" models
### these are what is being pooled together to form the error term. observe:
sum(A)
# [1] 16.75
### this is the pooled Residual Sum Sq from the RM-ANOVA

C - A
# .L    .Q 
# 30.25  3.00 
### and these are the SSRs for the contrasts in the RM-ANOVA

### now let's run the model specifying that we are assuming sphericity
### note that this is the same model object ("mod") from the first block of code above
anova(mod, test="Spherical")
# Analysis of Variance Table
# 
# Greenhouse-Geisser epsilon: 0.7492
# Huynh-Feldt epsilon:        0.9078
# 
#             Df      F num Df den Df     Pr(>F)    G-G Pr     H-F Pr
# (Intercept)  1 13.896      2     14 0.00047349 0.0018555 0.00078087
# Residuals    7                                                     
### note the equivalent F-ratio
### see also the sphericity tests and corrected p-values