Introduction

In a previous article, we presented the one-way ANOVA, the statistical test used to compare a quantitative variable between three groups or more. One of its assumptions was stated very explicitly in that article: the observations must be independent, both within and between the groups. It was also mentioned that if observations between samples are dependent (for example, if three measurements have been collected on the same individuals, as it is often the case in medical studies when a metric is measured (i) before, (ii) during and (iii) after a treatment), the repeated measures ANOVA should be preferred. The present article is dedicated to that test.

The repeated measures ANOVA compares the means of a quantitative variable measured several times on the same subjects, that is, under \(k \geq 3\) related conditions or at \(k \geq 3\) different points in time. The aim is exactly the same as the aim of the one-way ANOVA (testing whether the means are equal across the \(k\) conditions), but the design is different: instead of \(k\) independent groups formed by different subjects, we have a single group of subjects who go through all the conditions. The factor whose levels are the conditions is then called a within-subjects factor, as opposed to the between-subjects factor of the one-way ANOVA.

The relationship between the two tests mirrors the relationship between the two versions of the Student’s t-test:

In other words, the repeated measures ANOVA is to the paired Student’s t-test what the one-way ANOVA is to the Student’s t-test for independent samples: the generalization of the same logic to three groups or more. This is also why the repeated measures ANOVA appears in the list of related methods presented in the article about the two-way ANOVA, next to the mixed ANOVA (used when a between-subjects factor and a within-subjects factor are present at the same time).

Taking the repeated structure of the data into account is not a detail, and it is beneficial for two reasons:

1. Analyzing repeated measurements as if they came from independent groups violates the independence assumption of the one-way ANOVA, so the results of that test could simply not be trusted.
2. Since each subject serves as its own control, the variability between subjects (the fact that some patients are, in general, more sensitive to pain than others, for instance) is isolated and removed from the error term. The test is therefore usually more powerful than a one-way ANOVA run on the same number of measurements.

The price to pay for this gain is an additional assumption, called sphericity, which does not exist in the independent-groups case and which is discussed in detail in this article.

In the remaining of the post, we present the data, the aim, the hypotheses and the assumptions of the test, and we finally show how to perform it in R, how to complement it with post-hoc tests and how to interpret the results.

Data

Datasets with a genuinely repeated structure are not so common among the datasets shipped with R, so we simulate our own data. This has the additional advantage that we know exactly how the data have been generated.

Suppose that a treatment against chronic pain is administered to 30 randomly selected patients, and that the intensity of the pain is measured on each patient on a scale from 0 (no pain at all) to 100 (unbearable pain) at three different moments:

1. before the treatment,
2. during the treatment, and
3. one month after the end of the treatment, in order to see whether the benefit of the treatment persists over time.

# number of patients
n <- 30

# each patient has its own baseline level of pain
patient_effect <- rnorm(n, mean = 0, sd = 8)

# pain score at the 3 moments
before <- 70 + patient_effect + rnorm(n, mean = 0, sd = 6)
during <- 63 + patient_effect + rnorm(n, mean = 0, sd = 6)
after <- 61.5 + patient_effect + rnorm(n, mean = 0, sd = 6)

# dataset in the long format
dat <- data.frame(
  patient = factor(rep(1:n, times = 3)),
  time = factor(rep(c("before", "during", "after"), each = n),
    levels = c("before", "during", "after")
  ),
  pain = round(c(before, during, after), 1)
)

str(dat)

## 'data.frame':	90 obs. of  3 variables:
##  $ patient: Factor w/ 30 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ time   : Factor w/ 3 levels "before","during",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ pain   : num  83.7 69.7 79.1 71.4 76.3 58.8 77.4 64.1 71.7 69.7 ...

head(dat)

##   patient   time pain
## 1       1 before 83.7
## 2       2 before 69.7
## 3       3 before 79.1
## 4       4 before 71.4
## 5       5 before 76.3
## 6       6 before 58.8

(Note that a seed has been set in the background with set.seed(42), so the simulated data and all the results presented below are reproducible.)

Two points are worth being highlighted in the code above:

• The term patient_effect is added to the three measurements of the same patient. This is what creates the dependency between the three scores of a given patient: a patient with a high baseline level of pain tends to have a high score at all three moments.
• The data are stored in the long format, with one row per patient and per measurement, and three variables: the identifier of the subject (patient), the within-subjects factor (time) and the quantitative variable of interest (pain). This is the format expected by the functions used in the rest of the article. If your data are in the wide format (one row per subject and one column per condition), they can be reshaped with the pivot_longer() function of the {tidyr} package.

As always, it is a good practice to start with some descriptive statistics, here the mean and the standard deviation of the pain score at each of the three moments:

# install.packages("dplyr")
library(dplyr)

dat %>%
  group_by(time) %>%
  summarise(
    n = n(),
    mean = mean(pain),
    sd = sd(pain)
  )

## # A tibble: 3 × 4
##   time       n  mean    sd
##   <fct>  <int> <dbl> <dbl>
## 1 before    30  69.8  10.7
## 2 during    30  64.7  10.5
## 3 after     30  61.9  11.5

The boxplots below give a first visual comparison of the three moments:

# install.packages("ggplot2")
library(ggplot2)

ggplot(dat) +
  aes(x = time, y = pain) +
  geom_boxplot() +
  labs(
    x = "Moment of the measurement",
    y = "Pain score"
  )

Boxplots are, however, not entirely satisfactory here because they completely hide the repeated structure of the data. A plot which is much better suited to repeated measurements is the so-called spaghetti plot, where the successive scores of each patient are joined by a line (the thick blue line represents the mean at each moment):

ggplot(dat) +
  aes(x = time, y = pain, group = patient) +
  geom_line(alpha = 0.3) +
  geom_point(alpha = 0.3) +
  stat_summary(aes(group = 1),
    fun = mean, geom = "line",
    linewidth = 1.2, color = "steelblue"
  ) +
  stat_summary(aes(group = 1),
    fun = mean, geom = "point",
    size = 3, color = "steelblue"
  ) +
  labs(
    x = "Moment of the measurement",
    y = "Pain score"
  )

In our sample, the mean pain score decreased from 69.8 before the treatment to 64.7 during the treatment, and to 61.9 one month after it. The spaghetti plot also shows that most (but not all) patients follow this downward trend, and that the general level of pain varies a lot from one patient to another, which is precisely the between-subjects variability that the repeated measures ANOVA is able to set aside.

The question is now whether these differences are large enough to be generalized to the population, or whether they could be explained by sampling fluctuations alone. This is where the repeated measures ANOVA comes into play.

Aim and hypotheses

As explained in the introduction, the repeated measures ANOVA is used to compare the means of a quantitative variable measured on the same subjects under \(k \geq 3\) related conditions or at \(k \geq 3\) points in time.

The null and alternative hypotheses of the test are:

• \(H_0\): the population means are equal in all \(k\) related conditions, that is, \(\mu_1 = \mu_2 = \dots = \mu_k\)
• \(H_1\): at least one condition is different from the others in terms of mean

Be careful that, exactly as for the one-way ANOVA or the Kruskal-Wallis test, the alternative hypothesis is not that all means are different from each other. The opposite of “all means are equal” (\(H_0\)) is “at least one mean is different from the others” (\(H_1\)). So if the null hypothesis is rejected, we only know that at least one moment differs from the others, and post-hoc tests (covered later in this article) must be performed in order to know which ones actually differ.

In the context of our example, the repeated measures ANOVA helps us to answer the following question: “Is the mean pain score the same before, during and after the treatment?”. With the notations of our dataset, the hypotheses become:

• \(H_0\): \(\mu_{before} = \mu_{during} = \mu_{after}\)
• \(H_1\): at least one of the three moments differs from the others in terms of mean pain score

Note that, as its name suggests, the test still works by comparing variances: the variability observed between the conditions is compared to the residual variability, but this time after having removed the variability due to the subjects themselves. This is the reason why the error term of a repeated measures ANOVA is smaller than the error term of a one-way ANOVA computed on the same data, and why the test is generally more powerful.

Assumptions

As for many statistical tests, some assumptions must be met for the results to be valid.

# residuals of the model, taking the patient effect into account
res_lm <- lm(pain ~ time + patient, data = dat)

par(mfrow = c(1, 2)) # combine plots

hist(residuals(res_lm),
  main = "Histogram of the residuals",
  xlab = "Residuals"
)

# install.packages("car")
library(car)

qqPlot(residuals(res_lm),
  id = FALSE # remove point identification
)

shapiro.test(residuals(res_lm))

## 
## 	Shapiro-Wilk normality test
## 
## data:  residuals(res_lm)
## W = 0.98388, p-value = 0.331

# install.packages("tidyr")
library(tidyr)

# from the long format to the wide format
dat_wide <- dat %>%
  pivot_wider(names_from = time, values_from = pain)

c(
  "during - before" = var(dat_wide$during - dat_wide$before),
  "after - before" = var(dat_wide$after - dat_wide$before),
  "after - during" = var(dat_wide$after - dat_wide$during)
)

## during - before  after - before  after - during 
##        77.21154        67.05289        69.24033

Repeated measures ANOVA in R

# install.packages("rstatix")
library(rstatix)

res_aov <- anova_test(
  data = dat,
  dv = pain,
  wid = patient,
  within = time
)

res_aov

## ANOVA Table (type III tests)
## 
## $ANOVA
##   Effect DFn DFd      F        p p<.05   ges
## 1   time   2  58 13.464 1.57e-05     * 0.085
## 
## $`Mauchly's Test for Sphericity`
##   Effect     W   p p<.05
## 1   time 0.992 0.9      
## 
## $`Sphericity Corrections`
##   Effect   GGe      DF[GG]    p[GG] p[GG]<.05   HFe      DF[HF]    p[HF]
## 1   time 0.993 1.99, 57.57 1.67e-05         * 1.065 2.13, 61.78 1.57e-05
##   p[HF]<.05
## 1         *

get_anova_table(res_aov, correction = "GG")

## ANOVA Table (type III tests)
## 
##   Effect  DFn   DFd      F        p p<.05   ges
## 1   time 1.99 57.57 13.464 1.67e-05     * 0.085

summary(aov(pain ~ time + Error(patient / time),
  data = dat
))

## 
## Error: patient
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 29   8265     285               
## 
## Error: patient:time
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## time       2  958.2   479.1   13.46 1.57e-05 ***
## Residuals 58 2063.9    35.6                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Post-hoc tests

We have just shown that at least one moment differs from the others, but we still do not know which one(s). To answer this question, we need post-hoc tests (in Latin, “after this”, so after having obtained significant results for the ANOVA), also referred as multiple pairwise-comparison tests.

The logic is the same as the one presented for the one-way ANOVA: we compare the conditions two by two, and we adjust the \(p\)-values because performing several tests on the same data increases the risk of finding a significant difference by chance alone. With 3 moments, there are 3 pairs to compare, and the probability of observing at least one significant result purely by chance would already be \(1 - (1 - 0.05)^3 = 14.3\%\) without any adjustment.

The only difference with the one-way ANOVA is that the comparisons must take the pairing into account: instead of the Tukey HSD test (which compares independent groups), we perform paired Student’s t-tests on each pair of moments, together with an adjustment of the \(p\)-values for multiple comparisons. This is done with the pairwise_t_test() function of the {rstatix} package, with the arguments paired = TRUE and the Holm adjustment method:¹

dat %>%
  pairwise_t_test(pain ~ time,
    paired = TRUE,
    p.adjust.method = "holm"
  )

## # A tibble: 3 × 10
##   .y.   group1 group2    n1    n2 statistic    df         p   p.adj p.adj.signif
## * <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>     <dbl>   <dbl> <chr>       
## 1 pain  before during    30    30      3.19    29 0.00343   6.86e-3 **          
## 2 pain  before after     30    30      5.27    29 0.0000120 3.61e-5 ****        
## 3 pain  during after     30    30      1.82    29 0.0793    7.93e-2 ns

(Be careful that, with paired = TRUE, this function pairs the observations according to their order in the dataset, so the subjects must appear in the same order in each condition, which is the case here.)

It is the p.adj column (the \(p\)-values adjusted for multiple comparisons) which is of interest, and not the p column (the unadjusted \(p\)-values). These adjusted \(p\)-values must be compared to the desired significance level, here 5%.

Based on the output, we conclude that:

• the mean pain score differs significantly between before and during the treatment (adjusted \(p\)-value = 0.007),
• it differs significantly between before the treatment and one month after it (adjusted \(p\)-value < 0.001), and
• it does not differ significantly between during the treatment and one month after it (adjusted \(p\)-value = 0.079).

Combined with the descriptive statistics computed earlier, these post-hoc tests give a much more precise picture than the ANOVA alone: the treatment is associated with a significant decrease of the pain score (from 69.8 to 64.7 on average), and this benefit is still visible one month after the end of the treatment (61.9 on average, still significantly below the initial level). The additional decrease observed between the measurement during the treatment and the measurement one month later is, on the other hand, too small to be considered as significant.

Note that if the normality assumption had not been satisfied, the post-hoc tests would have been pairwise Wilcoxon signed-rank tests (pairwise_wilcox_test() with paired = TRUE) instead of paired t-tests, following a Friedman test instead of a repeated measures ANOVA.

Summary

In this article, we reviewed the aim, the hypotheses and the assumptions of the repeated measures ANOVA, the extension of the one-way ANOVA to related samples, used when the same subjects are measured on a quantitative variable under three conditions or more. Beyond the usual requirements of independence between subjects and normality of the residuals, this test relies on the sphericity assumption (the variances of the differences between all pairs of conditions must be equal), which is tested with Mauchly’s test and, if violated, dealt with by applying the Greenhouse-Geisser or the Huynh-Feldt correction to the degrees of freedom.

In practice, the test is easily performed in R with the anova_test() function of the {rstatix} package, which reports the ANOVA table, Mauchly’s test and both corrections at once. As for any ANOVA, a significant result only indicates that at least one condition differs from the others, so it must be followed by post-hoc tests, here pairwise paired t-tests with adjusted \(p\)-values. If normality or sphericity cannot be assumed, remember that the Friedman test is the nonparametric alternative, and that with independent samples the one-way ANOVA or the Kruskal-Wallis test should be preferred. See also this overview of the most common statistical tests if you hesitate between several methods.

Thanks for reading.

I hope this article helped you to understand the repeated measures ANOVA and how to perform it in R.

As always, if you have a question or a suggestion related to the topic covered in this article, please add it as a comment so other readers can benefit from the discussion.