Repeated Measures ANOVA using R

While so-called “between-subjects” ANOVA is absolutely straightforward in R, performing repeated measures (within-subjects) ANOVA is not so obvious. I have come across at least three different ways of performing repeated measures ANOVA in R. Which method you use depends on a number of things including (a) how you wish to perform follow-up tests (e.g. Tukey tests); (b) whether you care about (and want to correct for) heterogeneity of variances of differences between groups (sphericity).

Note: to adhere to the sum-to-zero convention for effect weights, you should always do this before running anovas in R:

> options(contrasts=c("contr.sum","contr.poly"))

Put the code right at the top of your script, run it each time you start a new R session.

This matters sometimes (not always). If you don’t do it, your sum of squares calculations may not match what you get, for example, in SPSS, or in many common textbooks on ANOVA (e.g. Maxwell & Delaney). See [obtaining the same ANOVA results in R as in SPSS – the difficulties with Type II and Type III sums of squares] for a more detailed account of why this is.

Of course, if you have some principled reason why you do not want to adhere to the sum-to-zero convention for effects weights, then by all means, use some other convention. The important thing to remember is that R, by default, does NOT adhere to the sum-to-zero convention, so if you fail to use the above code snippet, your results may not look like what you get in SPSS or in common texts on ANOVA.

On to some example data and procedures:

Let’s assume we have the following dataset: one dependent variable and a single repeated-measures factor with three levels:

> dv <- c(1,3,4,2,2,3,2,5,6,3,4,4,3,5,6)
> subject <- factor(c("s1","s1","s1","s2","s2","s2","s3","s3","s3",
+ "s4","s4","s4","s5","s5","s5"))
> myfactor <- factor(c("f1","f2","f3","f1","f2","f3","f1","f2","f3",
+ "f1","f2","f3","f1","f2","f3"))
> mydata <- data.frame(dv, subject, myfactor)

1. Univariate approach using aov()

This method gives you the “usual” looking ANOVA table

> am1 <- aov(dv ~ myfactor + Error(subject/myfactor), data=mydata)
> summary(am1)
Error: subject
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  4   12.4     3.1
Error: subject:myfactor
          Df  Sum Sq Mean Sq F value   Pr(>F)
myfactor   2 14.9333  7.4667  13.576 0.002683 **
Residuals  8  4.4000  0.5500
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

You get no test for sphericity, and if you try to use mauchly.test() you get an error message:

> mauchly.test(am1)
Error in UseMethod("mauchly.test", object) :
  no applicable method for "mauchly.test"

You also cannot use TukeyHSD() to perform posthoc Tukey tests:

> TukeyHSD(am1)
Error in UseMethod("TukeyHSD") : no applicable method for "TukeyHSD"

The way to perform posthoc tests is to form a linear contrast that tests the pairwise comparison of interest, compute an F ratio associated with the contrast, and then correct for Type-I error as you wish — e.g. if you want to use Tukey tests, convert the F into a q (studentized range statistic) and then look up a probability using ptukey(). The error term (denominator) of the F ratio is read off of the ANOVA table output. Let’s say you wish to compare the means of levels 1 and 2 of factor1:

> SScomp <- (mean(dv[myfactor=="f1"]) -
+ mean(dv[myfactor=="f2"]))^2
> dfcomp <- 1
> n <- 5 #5 subjects per group
> SSerr <- 4.4 #read off the ANOVA table
> dferr <- 8 # read off the ANOVA table
> MSerr <- SSerr / dferr
> Fcomp <- (n * SScomp/2)/MSerr
> pcomp <- 1-pf(Fcomp, dfcomp, dferr)
> Fcomp
[1] 11.63636
> pcomp
[1] 0.00920665

now compute a Tukey probability:

> qobs <- sqrt(2*Fcomp)
> ngroups <- length(levels(myfactor))
> ptk <- 1 - ptukey(qobs, ngroups, (ngroups - 1) * (n - 1))
> ptk
[1] 0.02234716

2. Univariate approach using lme()

This method fits a linear mixed-effects model, and produces output that looks more like a typical linear model that you get using lm(). First we need to load in the nlme package.

> require(nlme)
> am2 <- lme(dv ~ myfactor, random = ~1|subject/myfactor,
+ data=mydata)
> summary(am2)
Linear mixed-effects model fit by REML
 Data: NULL
       AIC      BIC    logLik
  50.62575 53.53519 -19.31288
Random effects:
 Formula: ~1 | subject
        (Intercept)
StdDev:   0.9219536
 Formula: ~1 | myfactor %in% subject

        (Intercept)  Residual
StdDev:   0.6727578 0.3120855

Fixed effects: dv ~ myfactor
            Value Std.Error DF  t-value p-value
(Intercept)   2.2 0.5291500  8 4.157611  0.0032
myfactorf2    1.6 0.4690417  8 3.411210  0.0092
myfactorf3    2.4 0.4690417  8 5.116815  0.0009
 Correlation:
           (Intr) myfct2
myfactorf2 -0.443
myfactorf3 -0.443  0.500
Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max
-0.48688951 -0.31208553  0.05125158  0.28645975  0.44662036 

Number of Observations: 15
Number of Groups:
              subject myfactor %in% subject
                    5                    15

Note that this output is not like a traditional ANOVA table, it is more like what you get from the lm() function. You can get it to spit out the F test for the myfactor main effect using anova():

> anova(am2)
            numDF denDF  F-value p-value
(Intercept)     1     8 60.40869  0.0001
myfactor        2     8 13.57575  0.0027

You can see the F-value and p-value for the myfactor effect are the same as in the ANOVA table output using aov() above in 1.

Note there is no test of sphericity, and mauchly.test() gives an error message:

> mauchly.test(am2)
Error in UseMethod("mauchly.test", object) :
  no applicable method for "mauchly.test"

with TukeyHSD() you also get no joy:

> TukeyHSD(am2)
Error in UseMethod("TukeyHSD") : no applicable method for "TukeyHSD"

To perform pairwise tests you can follow the same procedure as above in 1, OR you can do the following using the multcomp package and the glht() function:

> require(multcomp)
> summary(glht(am2,linfct=mcp(myfactor="Tukey")))
     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts
Fit: lme.formula(fixed = dv ~ myfactor, data = mydata, random = ~1 |
    subject/myfactor)
Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)
f2 - f1 == 0    1.600      0.469   3.411  0.00185 **
f3 - f1 == 0    2.400      0.469   5.117  < 1e-04 ***
f3 - f2 == 0    0.800      0.469   1.706  0.20308
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Note that this method essentially uses the error term from the model fit as a whole … in other words it essentially assumes the sphericity assumption holds.

3. Multivariate approach using lm() / mlm()

The first step is to reorganize the data into a matrix form where rows are subjects, and columns are levels of the repeated measures factor:

> dvm <- with(mydata, cbind(dv[myfactor=="f1"],
+ dv[myfactor=="f2"], dv[myfactor=="f3"]))
> dvm
     [,1] [,2] [,3]
[1,]    1    3    4
[2,]    2    2    3
[3,]    2    5    6
[4,]    3    4    4
[5,]    3    5    6

The next step is to define a multivariate linear model:

> mlm1 <- lm(dvm ~ 1)
> mlm1
Call:
lm(formula = dvm ~ 1)
Coefficients:
             [,1]  [,2]  [,3]
(Intercept)  2.2   3.8   4.6

Note how the coefficients are equal to the means of each group.

The next step is to define a variable that defines the design of the study:

> rfactor <- factor(c("f1","f2","f3"))
> rfactor
[1] f1 f2 f3
Levels: f1 f2 f3

The next step is to load the car package, which contains the Anova() function (note the capital A in Anova())

> library(car)

Now we define a new anova model object mlm1.aov by a function call to Anova(), and then output the summary. Note in the summary we pass a flag, multivariate=FALSE, to avoid seeing the multivariate tests:

> mlm1.aov <- Anova(mlm1, idata = data.frame(rfactor),
+ idesign = ~rfactor, type="III")
> summary(mlm1.aov, multivariate=FALSE)
Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
                 SS num Df Error SS den Df      F   Pr(>F)
(Intercept) 187.267      1   12.400      4 60.409 0.001477 **
rfactor      14.933      2    4.400      8 13.576 0.002683 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Mauchly Tests for Sphericity

        Test statistic p-value
rfactor        0.26171 0.13388
Greenhouse-Geisser and Huynh-Feldt Corrections
 for Departure from Sphericity
         GG eps Pr(>F[GG])
rfactor 0.57528    0.01537 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
         HF eps Pr(>F[HF])
rfactor 0.65851    0.01086 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

So we get a traditional ANOVA table with the F test for the myfactor effect (called rfactor). Note that the F-value and p-value are the same as in the univariate approach above in 1. We also get, however, a test of sphericity (Mauchly test), and we get two corrected p-values, a Greenhouse-Geisser correction, and a Huynh-Feldt correction. Unfortunately TukeyHSD() still will not work here, so we have to do any followup tests by one of the methods outlined above in 1. or 2… OR, if we are concerned about the sphericity assumption, one may simply do a t.test on the difference scores between the two groups of interest, e.g.:

> contrastD <- dvm[,1]-dvm[,2]
> t.test(contrastD)

    One Sample t-test

data:  contrastD
t = -3.1379, df = 4, p-value = 0.03492
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -3.0157148 -0.1842852
sample estimates:
mean of x
     -1.6

Note that this is a test uncorrected for Type-I error. You could of course transform the t-statistic into an F (since dfnum=1, F = t^2), then convert the t into a q, and compute ptukey, as above.

Good luck