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.
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.
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.
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.
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.
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