Repeated Measures ANOVA – How ANOVA is Calculated for a Repeated Measures Design: aov() vs lm() in R

anovalinear modelrrepeated measures

The title says it all, and I'm confused. The following runs a repeated measures aov() in R, and runs what I thought was an equivalent lm() call, but they return different error residuals (although the sums of squares are the same).

Clearly the residuals and fitted values from aov() are the ones used in the model, because their sums of squares add up to each of the model/residual sums of squares reported in summary(my.aov). So what are the actual linear models that are applied to a repeated measures design?

set.seed(1)
# make data frame,
# 5 participants, with 2 experimental factors, each with 2 levels
# factor1 is A, B
# factor2 is 1, 2
DF <- data.frame(participant=factor(1:5), A.1=rnorm(5, 50, 20), A.2=rnorm(5, 100, 20), B.1=rnorm(5, 20, 20), B.2=rnorm(5, 50, 20))

# get our experimental conditions
conditions <- names(DF)[ names(DF) != "participant" ]

# reshape it for aov
DFlong <- reshape(DF, direction="long", varying=conditions, v.names="value", idvar="participant", times=conditions, timevar="group")

# make the conditions separate variables called factor1 and factor2
DFlong$factor1 <- factor( rep(c("A", "B"), each=10) )
DFlong$factor2 <- factor( rep(c(1, 2), each=5) )

# call aov
my.aov <- aov(value ~ factor1*factor2 + Error(participant / (factor1*factor2)), DFlong)

# similar for an lm() call
fit <- lm(value ~ factor1*factor2 + participant, DFlong)

# what's aov telling us?
summary(my.aov)

# check SS residuals
sum(residuals(fit)^2)       # == 5945.668

# check they add up to the residuals from summary(my.aov)
2406.1 + 1744.1 + 1795.46   # == 5945.66

# all good so far, but how are the residuals in the aov calculated?
my.aov$"participant:factor1"$residuals

#clearly these are the ones used in the ANOVA:
sum(my.aov$"participant:factor1"$residuals ^ 2)

# this corresponds to the factor1 residuals here:
summary(my.aov)


# but they are different to the residuals reported from lm()
residuals(fit)
my.aov$"participant"$residuals
my.aov$"participant:factor1"$residuals
my.aov$"participant:factor1:factor2"$residuals

Best Answer

One way to think about it is to treat the situation as a 3-factorial between subjects ANOVA with IVs participant, factor1, factor2, and a cell size of 1. anova(lm(value ~ factor1*factor2*participant, DFlong)) calculates all the SS for all effects in this 3-way ANOVA (3 main effects, 3 first-order interactions, 1 second-order interaction). Since there's only 1 person in each cell, the full model has no errors, and the above call to anova() cannot compute F-tests. But the SS are the same as in the 2-factorial within design.

How does anova() actually compute the SS for an effect? Through sequential model comparisons (type I): It fits a restricted model without the effect in question, and an unrestricted model which includes that effect. The SS associated with this effect is the difference in error SS between both models.

# get all SS from the 3-way between subjects ANOVA
anova(lm(value ~ factor1*factor2*participant, DFlong))

dfL <- DFlong   # just a shorter name for your data frame
names(dfL) <- c("id", "group", "DV", "IV1", "IV2")   # shorter variable names

# sequential model comparisons (type I SS), restricted model is first, then unrestricted
# main effects first
anova(lm(DV ~ 1,      dfL), lm(DV ~ id,         dfL))  # SS for factor id
anova(lm(DV ~ id,     dfL), lm(DV ~ id+IV1,     dfL))  # SS for factor IV1
anova(lm(DV ~ id+IV1, dfL), lm(DV ~ id+IV1+IV2, dfL))  # SS for factor IV2

# now first order interactions
anova(lm(DV ~ id+IV1+IV2, dfL), lm(DV ~ id+IV1+IV2+id:IV1,  dfL))  # SS for id:IV1
anova(lm(DV ~ id+IV1+IV2, dfL), lm(DV ~ id+IV1+IV2+id:IV2,  dfL))  # SS for id:IV2
anova(lm(DV ~ id+IV1+IV2, dfL), lm(DV ~ id+IV1+IV2+IV1:IV2, dfL))  # SS for IV1:IV2

# finally the second-order interaction id:IV1:IV2
anova(lm(DV ~ id+IV1+IV2+id:IV1+id:IV2+IV1:IV2,            dfL),
      lm(DV ~ id+IV1+IV2+id:IV1+id:IV2+IV1:IV2+id:IV1:IV2, dfL))

Now let's check the effect SS associated with the interaction id:IV1 by subtracting the error SS of the unrestricted model from the error SS of the restricted model.

sum(residuals(lm(DV ~ id+IV1+IV2,        dfL))^2) -
sum(residuals(lm(DV ~ id+IV1+IV2+id:IV1, dfL))^2)

Now that you have all the "raw" effect SS, you can build the within-subjects tests simply by choosing the correct error term to test an effect SS against. E.g., test the effect SS for factor1 against the interaction effect SS of participant:factor1.

For an excellent introduction to the model comparison approach, I recommend Maxwell & Delaney (2004). Designing Experiments and Analyzing Data.

Related Solutions

Solved – How to specify specific contrasts for repeated measures ANOVA using car

This method is generally considered "old-fashioned" so while it may be possible, the syntax is difficult and I suspect fewer people know how to manipulate the anova commands to get what you want. The more common method is using glht with a likelihood-based model from nlme or lme4. (I'm certainly welcome to be proved wrong by other answers though.)

That said, if I needed to do this, I wouldn't bother with the anova commands; I'd just fit the equivalent model using lm, pick out the right error term for this contrast, and compute the F test myself (or equivalently, t test since there's only 1 df). This requires everything to be balanced and have sphericity, but if you don't have that, you should probably be using a likelihood-based model anyway. You might be able to somewhat correct for non-sphericity using the Greenhouse-Geiser or Huynh-Feldt corrections which (I believe) use the same F statistic but modify the df of the error term.

If you really want to use car, you might find the heplot vignettes helpful; they describe how the matrices in the car package are defined.

Using caracal's method (for the contrasts 1&2 - 3 and 1&2 - 4&5), I get

      psiHat      tStat          F         pVal
1 -3.0208333 -7.2204644 52.1351067 2.202677e-09
2 -0.2083333 -0.6098777  0.3719508 5.445988e-01

This is how I'd get those same p-values:

Reshape the data into long format and run lm to get all the SS terms.

library(reshape2)
d <- OBrienKaiser
d$id <- factor(1:nrow(d))
dd <- melt(d, id.vars=c(18,1:2), measure.vars=3:17)
dd$hour <- factor(as.numeric(gsub("[a-z.]*","",dd$variable)))
dd$phase <- factor(gsub("[0-9.]*","", dd$variable), 
                   levels=c("pre","post","fup"))
m <- lm(value ~ treatment*hour*phase + treatment*hour*phase*id, data=dd)
anova(m)

Make an alternate contrast matrix for the hour term.

foo <- matrix(0, nrow=nrow(dd), ncol=4)
foo[dd$hour %in% c(1,2) ,1] <- 0.5
foo[dd$hour %in% c(3) ,1] <- -1
foo[dd$hour %in% c(1,2) ,2] <- 0.5
foo[dd$hour %in% c(4,5) ,2] <- -0.5
foo[dd$hour %in% 1 ,3] <- 1
foo[dd$hour %in% 2 ,3] <- 0
foo[dd$hour %in% 4 ,4] <- 1
foo[dd$hour %in% 5 ,4] <- 0

Check that my contrasts give the same SS as the default contrasts (and the same as from the full model).

anova(lm(value ~ hour, data=dd))
anova(lm(value ~ foo, data=dd))

Get the SS and df for just the two contrasts I want.

anova(lm(value ~ foo[,1], data=dd))
anova(lm(value ~ foo[,2], data=dd))

Get the p-values.

> F <- 73.003/(72.81/52)
> pf(F, 1, 52, lower=FALSE)
[1] 2.201148e-09
> F <- .5208/(72.81/52)
> pf(F, 1, 52, lower=FALSE)
[1] 0.5445999

Optionally adjust for sphericity.

pf(F, 1*.48867, 52*.48867, lower=FALSE)
pf(F, 1*.57413, 52*.57413, lower=FALSE)

Solved – Power Analysis for repeated measures ANOVA

To calculate the effect size for a 2-way repeated ANOVA on both factors, you can use two formulas:

$$\eta^2_{partial} =\frac{SS} {SS+SS_{Error}}$$ where $SS$ is the sum of squares.

The number you'll get, must then be multiplied with 100, so you'll have the percentage of explained influence of the factor on your dependent variable.

The other way to do this, is with the formulas:

$f^2 = F \cdot \frac{ df}{df_{error}} \rightarrow \eta^2_{partial} = \frac{f^2} {1+f^2}$

Likewise you'll have to multiple your result with 100, so you'll get the explained variance.

Both these formulas tell you about the effect size for your sample, but not for the population.

Best Answer

Related Solutions

Solved – How to specify specific contrasts for repeated measures ANOVA using car

Solved – Power Analysis for repeated measures ANOVA

Related Question