Introduction

In this guide, we’ll delve into the application of repeated measures ANOVA using R.

Repeated measures ANOVA becomes especially powerful when the focus extends beyond the individual effects of independent variables (or factors) to encompass how they interact with each other. This method allows researchers to examine not only the main effects of each factor but also the potential interactions between them, providing a comprehensive understanding of their combined influence on the dependent variable.

Create alcohol example data

Let’s use the alcohol example from factorial ANOVA, but with a revised scenario - now every participant completes the driving course for all three levels of alcohol (sounds like a fun experiment!), separated by a week so we can be sure the participants are completely sober. Thus, alcohol is a within factor. We’ll also randomly assign each participant to drive slow or fast, so speed is a between factor.

cones <- c(0,1,0,2,1,2,3,0,3,1,4,6,3,2,3,0,0,3,2,0,4,5,2,3,2,7,9,6,4,6)
speed <- factor(c(rep("slow",15),rep("fast",15)))
id <- c(rep(c(1,2,3,4,5),3),rep(c(6,7,8,9,10),3))
alcohol <- factor(c(rep("drink1",5),rep("drink3",5),rep("drink5",5),
                    rep("drink1",5),rep("drink3",5),rep("drink5", 5)))
alcdata <- data.frame(id,cones,speed,alcohol)
knitr::kable(alcdata, format = "markdown")
id cones speed alcohol
1 0 slow drink1
2 1 slow drink1
3 0 slow drink1
4 2 slow drink1
5 1 slow drink1
1 2 slow drink3
2 3 slow drink3
3 0 slow drink3
4 3 slow drink3
5 1 slow drink3
1 4 slow drink5
2 6 slow drink5
3 3 slow drink5
4 2 slow drink5
5 3 slow drink5
6 0 fast drink1
7 0 fast drink1
8 3 fast drink1
9 2 fast drink1
10 0 fast drink1
6 4 fast drink3
7 5 fast drink3
8 2 fast drink3
9 3 fast drink3
10 2 fast drink3
6 7 fast drink5
7 9 fast drink5
8 6 fast drink5
9 4 fast drink5
10 6 fast drink5

It’s easy to rearrange the data.

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
alcdata <- alcdata %>% arrange(id,alcohol)
knitr::kable(alcdata, format = "markdown")
id cones speed alcohol
1 0 slow drink1
1 2 slow drink3
1 4 slow drink5
2 1 slow drink1
2 3 slow drink3
2 6 slow drink5
3 0 slow drink1
3 0 slow drink3
3 3 slow drink5
4 2 slow drink1
4 3 slow drink3
4 2 slow drink5
5 1 slow drink1
5 1 slow drink3
5 3 slow drink5
6 0 fast drink1
6 4 fast drink3
6 7 fast drink5
7 0 fast drink1
7 5 fast drink3
7 9 fast drink5
8 3 fast drink1
8 2 fast drink3
8 6 fast drink5
9 2 fast drink1
9 3 fast drink3
9 4 fast drink5
10 0 fast drink1
10 2 fast drink3
10 6 fast drink5

Converting to wide format should be really easy. It is, sometimes. The more complex the data, the harderit is. This is easy.

library(tidyr)
alcdata_wide <- alcdata %>%
  pivot_wider(names_from = alcohol, values_from = cones)
knitr::kable(alcdata_wide, format = "markdown")
id speed drink1 drink3 drink5
1 slow 0 2 4
2 slow 1 3 6
3 slow 0 0 3
4 slow 2 3 2
5 slow 1 1 3
6 fast 0 4 7
7 fast 0 5 9
8 fast 3 2 6
9 fast 2 3 4
10 fast 0 2 6

Same idea with converting back to long format.

alcdata_long <- alcdata_wide %>%
  pivot_longer(cols = starts_with("drink"), names_to = "alcohol", values_to = "cones", names_prefix = "drink")
knitr::kable(alcdata_long, format = "markdown")
id speed alcohol cones
1 slow 1 0
1 slow 3 2
1 slow 5 4
2 slow 1 1
2 slow 3 3
2 slow 5 6
3 slow 1 0
3 slow 3 0
3 slow 5 3
4 slow 1 2
4 slow 3 3
4 slow 5 2
5 slow 1 1
5 slow 3 1
5 slow 5 3
6 fast 1 0
6 fast 3 4
6 fast 5 7
7 fast 1 0
7 fast 3 5
7 fast 5 9
8 fast 1 3
8 fast 3 2
8 fast 5 6
9 fast 1 2
9 fast 3 3
9 fast 5 4
10 fast 1 0
10 fast 3 2
10 fast 5 6

Let’s get some descriptives.

alcdata_long %>%
group_by(alcohol) %>%
  summarize(count = n(),
            mean = mean(cones, na.rm = TRUE),
            sd = sd(cones, na.rm = TRUE))
## # A tibble: 3 × 4
##   alcohol count  mean    sd
##   <chr>   <int> <dbl> <dbl>
## 1 1          10   0.9  1.10
## 2 3          10   2.5  1.43
## 3 5          10   5    2.16

There are ways to test homogeneity of covariance, but as far as I know they are always in the context of testing homogeneity of variance and covariance simultaneously. The most common is Box’s M, which has known problems. More common are direct tests of sphericity, which are later in this script. For now, we’ll just eyeball it.

cor(alcdata_wide[,3:5]) # get the correlation matrix
##             drink1      drink3     drink5
## drink1  1.00000000 -0.03521039 -0.2804228
## drink3 -0.03521039  1.00000000  0.6816212
## drink5 -0.28042282  0.68162120  1.0000000
cov(alcdata_wide[,3:5]) # get the covariance matrix
##             drink1      drink3     drink5
## drink1  1.21111111 -0.05555556 -0.6666667
## drink3 -0.05555556  2.05555556  2.1111111
## drink5 -0.66666667  2.11111111  4.6666667

My eyeballs tell me that there may be issues.

Repeated measures ANOVA with aov

To run repeated measures ANOVA in aov, you have to specify how the error will work. That involves adding a term to the model labeled Error (capital E) with an indicator variable of the individual/matched group (here, the indicator variable is id) in the numerator and the name of the within factor in the denominator (here, the within factor is alcohol). This syntax isn’t totally intuitive but at least it’s not complicated.

Let’s start easy - just one predictor, the within factor of alcohol.

output<-aov(cones ~ alcohol + Error(factor(id)/alcohol),data=alcdata)
summary(output)
## 
## Error: factor(id)
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  9  32.13    3.57               
## 
## Error: factor(id):alcohol
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## alcohol    2  85.40   42.70   19.57 3.05e-05 ***
## Residuals 18  39.27    2.18                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The interpretation works just like between-subjects ANOVA even though the output looks more complicated. There were significant differences in cones hit across alcohol conditions. Examination of means indicated that more alcohol was associated with more cones hit.

The ANOVA summary table becomes increasingly complicated and therefore irrelevant the more within factors there are. That’s because there is no longer a single error term. I recommend just reporting the F with associated df and p-value.

With a repeated measures ANOVA, post-hoc tests work the same, although there a lot of potential problems especially if sphericity is violated. There is some evidence that a Bonferroni post-hoc test is the best choice with repeated measures ANOVA, but the evidence is nowhere near solid.

pairwise.t.test(alcdata$cones,alcdata$alcohol,p.adjust.method = "bonferroni")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  alcdata$cones and alcdata$alcohol 
## 
##        drink1  drink3
## drink3 0.1096  -     
## drink5 1.7e-05 0.0058
## 
## P value adjustment method: bonferroni

The post-hoc tests indicate that 1 drink is significantly different from 5 drinks and 3 drinks is significantly different from 5 drinks. 1 drink and 3 drink are not significantly different from each other.

The multcomp function doesn’t work with repeated measures ANOVA. There is no easy way to do the Tukey post-hoc test that I know of.

Repeated measures ANOVA with ezANOVA

With repeated measures ANOVA, the ezANOVA function is a lot easier to work with. It is in the ez package.

library(ez)
output <- ezANOVA(alcdata_long,dv = .(cones),wid = .(id),within=.(alcohol))
## Warning: Converting "id" to factor for ANOVA.
## Warning: Converting "alcohol" to factor for ANOVA.
output
## $ANOVA
##    Effect DFn DFd        F            p p<.05       ges
## 2 alcohol   2  18 19.57385 3.051147e-05     * 0.5446429
## 
## $`Mauchly's Test for Sphericity`
##    Effect         W          p p<.05
## 2 alcohol 0.5603581 0.09859677      
## 
## $`Sphericity Corrections`
##    Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2 alcohol 0.6946172 0.0003460526         * 0.7812843 0.0001731119         *

The ezANOVA function automatically gives the most common test of sphericity, called Mauchly’s test. Here, p = .099 for Mauchly’s test, so no evidence that sphericity is violated. However, the sample size is small so I wouldn’t totally trust this result.

The ezANOVA function also automatically gives results with corrections for violations of sphericity. I don’t know of any clear evidence that one is better than the other, so pick one or the other to report. But not both- there is no point to reporting both. In general, I think it is better to automatically use one of the corrections, since repeated-measures ANOVA is sensitive to violations of sphericity.

Here, the result for the ANOVA with the Greenhouse-Geisser correction is p = .000346; with the Huynh-Feldt correction, p = .000173. In this case, and in most cases, the choice doesn’t matter: there is a significant effect of alcohol.

Add speed as a between factor

Adding speed as a between factor is easy. No change to the error term. Just add speed in as a predictor, connected to alcohol with a * if you want to consider an interaction.

output<-aov(cones ~ alcohol*speed + Error(factor(id)/alcohol),data=alcdata_long)
summary(output)
## 
## Error: factor(id)
##           Df Sum Sq Mean Sq F value Pr(>F)  
## speed      1  16.13   16.13   8.067 0.0218 *
## Residuals  8  16.00    2.00                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: factor(id):alcohol
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## alcohol        2  85.40   42.70  22.182 2.44e-05 ***
## alcohol:speed  2   8.47    4.23   2.199    0.143    
## Residuals     16  30.80    1.92                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There are different error terms for the between and within factors. Note that the interaction is considered within if it involves a within factor, so the results look a little different. The between factors will be in a different section than the within factors, but the interpretation of each result does not change.

Effect sizes work the same.

library(sjstats)
knitr::kable(anova_stats(output)) %>%
  kableExtra::kable_styling()
etasq partial.etasq omegasq partial.omegasq epsilonsq cohens.f group term sumsq df meansq statistic p.value power
speed 0.103 0.344 0.090 0.197 0.091 0.724 factor(id) speed 16.133 1 16.133 8.067 0.022 0.822
Residuals 0.102 0.342 0.004 0.010 0.004 0.721 factor(id) Residuals 16.000 8 2.000 NA NA 0.504
alcohol 0.545 0.735 0.514 0.585 0.520 1.665 factor(id):alcohol alcohol 85.400 2 42.700 22.182 0.000 1.000
alcohol:speed 0.054 0.216 0.029 0.074 0.029 0.524 factor(id):alcohol alcohol:speed 8.467 2 4.233 2.199 0.143 0.445
1 NA NA NA NA NA NA factor(id):alcohol Residuals 30.800 16 1.925 NA NA NA

The output is messy, so be careful that columns line up. The easiest thing I’ve found is to cut and paste the results into Word and adjust fonts - change the font to a fixed width font (Courier New is a good one) and reduce the font size until things stay on a single line.

Focus on the term column to figure out which row to use. There are effect sizes associated with the different error terms (residuals) but they aren’t what you want. Here, the \(\omega^2\) values for speed, alcohol, and the interaction are .090, .514, and .029, respectively.

The syntax for ezANOVA is easier.

output <- ezANOVA(alcdata_long,dv = .(cones),wid = .(id),within=.(alcohol),between=.(speed))
## Warning: Converting "id" to factor for ANOVA.
## Warning: Converting "alcohol" to factor for ANOVA.
output
## $ANOVA
##          Effect DFn DFd         F            p p<.05       ges
## 2         speed   1   8  8.066667 2.180567e-02     * 0.2563559
## 3       alcohol   2  16 22.181818 2.436446e-05     * 0.6459909
## 4 speed:alcohol   2  16  2.199134 1.432892e-01       0.1531966
## 
## $`Mauchly's Test for Sphericity`
##          Effect         W         p p<.05
## 3       alcohol 0.6584022 0.2315898      
## 4 speed:alcohol 0.6584022 0.2315898      
## 
## $`Sphericity Corrections`
##          Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF]
## 3       alcohol 0.7453799 0.0001973507         * 0.8769716 6.675607e-05
## 4 speed:alcohol 0.7453799 0.1602423534           0.8769716 1.513548e-01
##   p[HF]<.05
## 3         *
## 4

The interpretation is straightforward. There are significant effects of speed and alcohol. The interaction is not significant.

Just for demonstration purposes, I recoded the data to make speed a within factor as well.

alcdata_long$withinid <- ifelse(alcdata_long$id<6,alcdata_long$id,alcdata_long$id-5) 

Now, each within factor needs its own error term. This gets put into the denominator in the aov function.

output<-aov(cones ~ alcohol*speed + Error(factor(withinid)/(alcohol*speed)),data=alcdata_long)
summary(output)
## 
## Error: factor(withinid)
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  4  12.47   3.117               
## 
## Error: factor(withinid):alcohol
##           Df Sum Sq Mean Sq F value Pr(>F)   
## alcohol    2  85.40   42.70   12.68 0.0033 **
## Residuals  8  26.93    3.37                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: factor(withinid):speed
##           Df Sum Sq Mean Sq F value Pr(>F)  
## speed      1 16.133  16.133   18.26 0.0129 *
## Residuals  4  3.533   0.883                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: factor(withinid):alcohol:speed
##               Df Sum Sq Mean Sq F value  Pr(>F)   
## alcohol:speed  2  8.467   4.233   8.759 0.00966 **
## Residuals      8  3.867   0.483                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Since each factor has its own error term, there are lots of sections to the output, but interpretation is the same.

It’s even easier to do this in ezANOVA.

output <- ezANOVA(alcdata_long,
                  dv = .(cones),
                  wid = .(withinid),
                  within=.(alcohol,speed))
## Warning: Converting "withinid" to factor for ANOVA.
## Warning: Converting "alcohol" to factor for ANOVA.
output
## $ANOVA
##          Effect DFn DFd         F           p p<.05       ges
## 2       alcohol   2   8 12.683168 0.003304653     * 0.6459909
## 3         speed   1   4 18.264151 0.012912536     * 0.2563559
## 4 alcohol:speed   2   8  8.758621 0.009661066     * 0.1531966
## 
## $`Mauchly's Test for Sphericity`
##          Effect         W         p p<.05
## 2       alcohol 0.6613322 0.5378108      
## 4 alcohol:speed 0.2853746 0.1524484      
## 
## $`Sphericity Corrections`
##          Effect       GGe       p[GG] p[GG]<.05       HFe       p[HF] p[HF]<.05
## 2       alcohol 0.7470113 0.008819606         * 1.0914132 0.003304653         *
## 4 alcohol:speed 0.5832178 0.032411654         * 0.6762115 0.024622139         *

What happened to the types of sums of squares? ezANOVA gives us easy control over the SS. Let’s try Type III. Note that the syntax uses the number 3 for Type III.

output <- ezANOVA(alcdata_long,
                  dv = .(cones),
                  wid = .(withinid),
                  within=.(alcohol,speed),
                  type=3)
## Warning: Converting "withinid" to factor for ANOVA.
## Warning: Converting "alcohol" to factor for ANOVA.
output
## $ANOVA
##          Effect DFn DFd         F           p p<.05       ges
## 2       alcohol   2   8 12.683168 0.003304653     * 0.6459909
## 3         speed   1   4 18.264151 0.012912536     * 0.2563559
## 4 alcohol:speed   2   8  8.758621 0.009661066     * 0.1531966
## 
## $`Mauchly's Test for Sphericity`
##          Effect         W         p p<.05
## 2       alcohol 0.6613322 0.5378108      
## 4 alcohol:speed 0.2853746 0.1524484      
## 
## $`Sphericity Corrections`
##          Effect       GGe       p[GG] p[GG]<.05       HFe       p[HF] p[HF]<.05
## 2       alcohol 0.7470113 0.008819606         * 1.0914132 0.003304653         *
## 4 alcohol:speed 0.5832178 0.032411654         * 0.6762115 0.024622139         *

The results are the same; type III is default in ezANOVA. Let’s try Type I.

output <- ezANOVA(alcdata_long,
                  dv = .(cones),
                  wid = .(withinid),
                  within=.(alcohol,speed),
                  type=1)
## Warning: Converting "withinid" to factor for ANOVA.
## Warning: Converting "alcohol" to factor for ANOVA.
output
## $ANOVA
##          Effect DFn DFd         F           p p<.05       ges
## 1       alcohol   2   8 12.683168 0.003304653     * 0.7132517
## 2         speed   1   4 18.264151 0.012912536     * 0.3196830
## 3 alcohol:speed   2   8  8.758621 0.009661066     * 0.1978193

And Type II.

output <- ezANOVA(alcdata_long,
                  dv = .(cones),
                  wid = .(withinid),
                  within=.(alcohol,speed),
                  type=2)
## Warning: Converting "withinid" to factor for ANOVA.
## Warning: Converting "alcohol" to factor for ANOVA.
output
## $ANOVA
##          Effect DFn DFd         F           p p<.05       ges
## 2       alcohol   2   8 12.683168 0.003304653     * 0.6459909
## 3         speed   1   4 18.264151 0.012912536     * 0.2563559
## 4 alcohol:speed   2   8  8.758621 0.009661066     * 0.1531966
## 
## $`Mauchly's Test for Sphericity`
##          Effect         W         p p<.05
## 2       alcohol 0.6613322 0.5378108      
## 4 alcohol:speed 0.2853746 0.1524484      
## 
## $`Sphericity Corrections`
##          Effect       GGe       p[GG] p[GG]<.05       HFe       p[HF] p[HF]<.05
## 2       alcohol 0.7470113 0.008819606         * 1.0914132 0.003304653         *
## 4 alcohol:speed 0.5832178 0.032411654         * 0.6762115 0.024622139         *

The choice of SS depends on research design and research questions in the same way as with factorial ANOVA from the earlier lecture.

Cluster-robust standard errors

Independence of observations is a critical assumption. If we ignore the clustering by treating alcohol as if it is a between factor, we will tend to get underestimated SEs, leading to p-values that are too low, leading to a Type I error rate that is greater than the alpha (.05).

output<-lm(cones ~ alcohol,data=alcdata)
summary(output)
## 
## Call:
## lm(formula = cones ~ alcohol, data = alcdata)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##   -3.0   -0.9   -0.2    1.0    4.0 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     0.9000     0.5142   1.750   0.0915 .  
## alcoholdrink3   1.6000     0.7272   2.200   0.0365 *  
## alcoholdrink5   4.1000     0.7272   5.638 5.52e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.626 on 27 degrees of freedom
## Multiple R-squared:  0.5446, Adjusted R-squared:  0.5109 
## F-statistic: 16.15 on 2 and 27 DF,  p-value: 2.442e-05

Here is the syntax for a cluster robust SEs.

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(sandwich)
coeftest(output, vcov = vcovCL, cluster = id)
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)    0.90000    0.36067  2.4954 0.0189940 *  
## alcoholdrink3  1.60000    0.44240  3.6166 0.0012088 ** 
## alcoholdrink5  4.10000    0.90679  4.5215 0.0001102 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The SEs are different and these results are more valid than if the clustering were ignored.

Another option is the lm.cluster function in miceadds. It uses the same functions as the coeftest function, but with a more intuitive syntax.

library(miceadds)
## Loading required package: mice
## 
## Attaching package: 'mice'
## The following object is masked from 'package:stats':
## 
##     filter
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
## * miceadds 3.17-44 (2024-01-08 19:08:24)
output <- lm.cluster(cones ~ alcohol,data=alcdata, cluster = "id")
summary(output)
## R^2= 0.54464 
## 
##               Estimate Std. Error  t value     Pr(>|t|)
## (Intercept)        0.9  0.3606692 2.495361 1.258290e-02
## alcoholdrink3      1.6  0.6023274 2.656362 7.898867e-03
## alcoholdrink5      4.1  0.8800720 4.658710 3.181976e-06