Repeated-measures ANOVA in R: a worked example

A repeated-measures ANOVA tests whether the mean of one outcome differs across several conditions or time points when the same subjects are measured at each. In R the work is in the setup: long-format data, a factor for the within-subject variable, and a correction for sphericity. This post runs the whole analysis on R’s ChickWeight data and then fits the equivalent mixed model so you can see both outputs from one script.

The example uses the first four weighing days of an early-growth experiment. It sits beside the broader comparison of repeated-measures ANOVA and mixed models, which explains when to choose each. Here the focus is the R code and how to read every line of the output.

When a repeated-measures ANOVA fits

The design has one outcome measured on a continuous scale, one within-subject factor with three or more levels, and every subject measured at every level. The levels might be time points, doses or conditions. Subjects are independent of one another, but the measurements from one subject are correlated, which is why a plain one-way ANOVA is wrong here: it would treat four weighings of one chick as four unrelated chicks. The independence assumption explains the cost.

With only two levels, a paired t-test does the same job. With a between-subject factor as well, such as treatment arm, you have a mixed design, and the engine’s repeated-measures ANOVA has no between-subject factor to offer. The two-way ANOVA guide covers designs where subjects fall into groups without repeated measurements.

The test asks one omnibus question: are the mean outcomes at all levels equal? A clear answer says at least one level differs, not which. Pairwise follow-ups come afterwards, as in post hoc tests after ANOVA.

Get the data into long format

R’s repeated-measures tools differ on layout, which is the most common source of errors. The functions below, aov() and afex::aov_ez(), want long format: one row per subject and condition, with a column for the subject, a column for the condition and a column for the outcome. ChickWeight already ships that way, with Chick, Time and weight.

data(ChickWeight)
cw <- as.data.frame(ChickWeight)
cw$Chick <- factor(cw$Chick, ordered = FALSE)
cw <- subset(cw, Time %in% c(0, 2, 4, 6))   # first four weighing days
cw$Day <- factor(cw$Time)                    # a factor, not a number

The last line is the one that matters most. If Time stays numeric, a model fits a straight-line slope and has one degree of freedom for time. A factor gives one level per day and the three degrees of freedom a four-level repeated-measures ANOVA should have. Multivariate tools such as anova() on an mlm object want wide format instead, and the long versus wide layout post shows how to reshape between them.

Check the assumptions before you run it

A repeated-measures ANOVA rests on four assumptions, and each has a practical check.

Independence between subjects. One chick’s weights must not depend on another’s. Chicks reared in the same pen or hatched in the same batch can break this, and the data alone rarely show it. Decide it from the design, and if subjects are clustered, add that level to a mixed model.

Normality of the outcome within each level. With four levels and a few dozen subjects, the F test is fairly forgiving. Plot the weights by day, or the differences between days, and look for heavy skew or a few extreme points rather than relying on a normality test. The engine flagged 9 weights beyond Tukey’s fences in this example and changed nothing; the right response to an outlier is to check whether it is a real chick or a typing error.

No extreme outliers. One mistyped weight inflates the variance of its day and distorts every comparison involving that day. Look at the boxplot for each level before fitting.

Sphericity. The variances of the differences between every pair of days should be equal. This is the assumption that fails most often, and the correction in the next sections handles it. Plot each chick’s trajectory first: lines that fan out, as they do here, suggest the assumption will fail.

A fifth condition is design, not distribution: every subject needs a value at every level. Count the missing cells before you start, because each one removes a whole subject from the classical analysis.

How the variation is split

The repeated-measures ANOVA divides the total variation in weight into three parts. One part lies between chicks, because some chicks are simply heavier than others. A second lies between days, which is the effect being tested. The third is what remains after both are removed, the residual, which measures how much a chick’s pattern across days departs from the average pattern.

The day effect is tested against that residual, not against the raw scatter of all 196 weighings. Because the between-chick variation has been removed from the denominator, the test has more power than a one-way ANOVA would, which is the gain from measuring everyone repeatedly. The degrees of freedom follow from the design. With 196 weighings the total is 195. Of these, 48 belong to chicks (49 minus 1), 3 to days (4 minus 1), and the remaining 144 to the residual (48 times 3). Those are the uncorrected values the correction later multiplies by epsilon.

Run the analysis in R

Base R fits the classical analysis with an error stratum for the subject. It needs a complete, balanced data set, so keep only the chicks weighed on all four days first:

n_days <- table(cw$Chick)
cc <- droplevels(subset(cw, Chick %in% names(n_days)[n_days == 4]))
fit <- aov(weight ~ Day + Error(Chick/Day), data = cc)
summary(fit)

The term Error(Chick/Day) tells R that days are nested within chicks, so the day effect is tested against the variation between days within a chick. The result is the uncorrected F test. It assumes sphericity, the condition that all pairwise differences between days have equal variance, and it has no argument to relax it. Left in, the incomplete chick makes R warn that the Error() model is singular.

The afex package adds the correction and an effect size in one call:

afex::aov_ez(id = "Chick", dv = "weight", within = "Day", data = cw,
             anova_table = list(correction = "GG", es = "ges"))

Here correction = "GG" requests the Greenhouse-Geisser correction and es = "ges" requests generalized eta-squared, both stated explicitly so the call does not depend on package defaults. The function warns when it removes subjects with missing cells; read that warning rather than suppressing it. The sphericity post shows the base R route on a wide-format multivariate fit: mauchly.test() for Mauchly’s test and anova(..., X = ~1, test = "Spherical") for both epsilons.

What the engine computed

The same model was run through the Celeus engine, restricted to days 0, 2, 4 and 6 by a recorded population filter. The repeated-measures ANOVA method applies the Greenhouse-Geisser correction on every run. ChickWeight has 50 chicks, and one chick lacks the day 4 and day 6 weighings. The engine dropped it whole, leaving 196 weighings, which is 49 chicks times four days. The engine flagged that dropped chick and reported 9 weights beyond Tukey’s outlier fences, descriptive only, with nothing removed.

The engine’s own note on the drop says that one subject with a missing cell was dropped and that the linear mixed model is preferable for unbalanced data. A later section returns to that advice.

The weighing days differ: F(1.24, 59.42) = 571.7, p < 0.001. Nobody needed a test to expect chicks to grow in their first week, so the numbers are best used for what they show about the method. The full output and the script that produced it are listed beneath the article.

Read the corrected degrees of freedom

With four days and 49 chicks, an uncorrected test would use 3 and 144 degrees of freedom: four days minus one, and 48 chicks times 3. The reported values are 1.24 and 59.42. Dividing either corrected value by its uncorrected counterpart gives the Greenhouse-Geisser epsilon, a little over four tenths (derived here from the reported degrees of freedom: 1.24 / 3). For four levels the smallest possible epsilon is 1/3, so the estimate sits close to the floor and signals a large sphericity violation.

Growth data typically produce a violation like this one. Chicks start at similar weights and spread apart as they grow, so differences between distant days vary far more than differences between neighbouring ones. An uncorrected test would be miscalibrated, and the correction cut the numerator degrees of freedom by more than half. Here the F value is so large that the conclusion does not move. In a study with a modest effect, the same cut can move p across 0.05.

The engine reports no Mauchly test, because a preliminary test is a weak gate: it has little power in small samples and flags trivial departures in large ones. The correction is applied every time, and the size of the correction is the diagnostic. The sphericity post explains this reasoning in full.

Report the effect size

Generalized eta-squared was 0.8483 (95% CI 0.8066 to 0.8777), and partial eta-squared was 0.9225 (95% CI 0.901 to 0.9377). Generalized eta-squared is the one to report for a within-subject factor, because it keeps variation between chicks in the denominator and so stays comparable with between-subject designs (Bakeman, 2005). The eta-squared post explains the difference between the two measures.

Both intervals are computed from the uncorrected degrees of freedom, as the engine’s output notes, so they do not carry the sphericity correction. Say so when you report them. A p value alone says only that some difference was detected; the effect size says how much of the variation the days account for. Here the days account for 84.83 percent of the total variation including differences between chicks, a large effect by conventional benchmarks, which are rough guides and not thresholds of importance.

For a copy-ready sentence and a reporting template, see how to report ANOVA. A write-up of this result states the test, the correction, both degrees of freedom, the statistic, the exact or bounded p value, the effect size with its interval, and the number of subjects analysed.

Write up the result

A methods and results sentence for this analysis names the design, the correction and the effect size, and it states how many subjects were analysed. For example: “Mean weight differed across the first four weighing days, F(1.24, 59.42) = 571.7, p < 0.001 (Greenhouse-Geisser corrected), generalized eta-squared = 0.8483, 95% CI 0.8066 to 0.8777, among the 49 chicks weighed on all four days; one chick with missing weighings was excluded.” Adding the restriction to days 0 to 6 matters, because the result says nothing about later weighings.

If you used a different tool, say which and which version, and name any default you relied on. A reader should be able to rerun the call and get the same degrees of freedom. The script beneath the article does that.

Fit the same data as a mixed model

The repeated-measures ANOVA threw away a chick that had two usable weighings. A linear mixed model keeps it. The companion run fits weight on day as a four-level factor with a random intercept for each chick, the mixed-model counterpart of the ANOVA’s subject term:

library(lmerTest)
mm <- lmer(weight ~ Day + (1 | Chick), data = cw, REML = TRUE)
anova(mm, type = 3, ddf = "Kenward-Roger")

The model used 198 weighings, the 196 from the 49 complete chicks plus the 2 from the incomplete one. The Kenward-Roger test for day gave F(3, 145.4) = 560.5, p < 0.001. Kenward and Roger (1997) developed the small-sample adjustment behind that denominator degrees of freedom. The test statistic is close to the ANOVA’s 571.7 and the conclusion is identical, which is the usual outcome when the data are nearly complete and the time points are fixed.

The model also estimates how chicks differ. The intraclass correlation was 0.3784 (95% CI 0.2292 to 0.5335): of the variation left once the day means are accounted for, 37.84 percent lies between chicks rather than within them. The standard deviation of the chick intercepts was 3.313 grams and the residual standard deviation was 4.246 grams. A repeated-measures ANOVA does not report these.

Compare the two outputs

The two analyses answer the same omnibus question with different machinery:

The mixed model with only a random intercept implies compound symmetry, a special case of sphericity, so it is not assumption-free. These data violate sphericity badly, so the mixed model’s F(3, 145.4) carries no correction and keeps all 3 numerator degrees of freedom, like the uncorrected ANOVA rather than the corrected one. Adding a random slope for time, or modelling the residual covariance directly, relaxes it. The comparison guide discusses those choices and the missing-data assumptions behind each method.

The mixed model’s checks were mixed. Shapiro-Wilk on the residuals detected non-normality, while the check on the chick intercepts did not. The engine’s warning says to interpret p values cautiously and consider transforming the outcome, though with 198 observations the central limit theorem softens the concern. A detected departure is information about the fit, not a reason to hide the result.

Follow up with pairwise comparisons

An omnibus test does not say which days differ. The mixed model’s engine output includes Tukey-adjusted pairwise contrasts from estimated marginal means. Day 0 averaged 8.16 grams less than day 2 (95% CI for the difference -10.37 to -5.953) and 33.15 grams less than day 6 (95% CI -35.38 to -30.93). All six adjusted contrasts were clear, each with an adjusted p value far below 0.001.

Six pairwise comparisons from four levels make multiplicity real, which is why the contrasts are adjusted. Run the omnibus test first and the comparisons second, and report that the follow-ups were adjusted. With a growth curve like this one, a trend question, such as the average gain per day, is often more informative than six pairwise differences.

Which R route to choose

Choose by what you need to report. Base aov() with an error stratum is enough for a quick look at a balanced design, and it ships with R, but it gives no correction. The afex call returns the corrected test, the effect size and the missing-cell warning in one table, which is why it suits most reports. Calling summary() on the afex result also prints Mauchly’s test and both epsilons. The multivariate route through mauchly.test() and anova() on an mlm object gives the same numbers without extra packages, at the price of reshaping to wide format first.

A mixed model through lme4 and lmerTest is the better choice once any subject is incomplete, visits are unevenly timed, or the sphericity pattern itself matters. The cost is a longer specification: you decide the random effects and the denominator degrees of freedom method, and you state both. Whatever you pick, fix the choice before looking at the results, and report it in the methods section.

Mistakes to avoid

• Leaving the within-subject variable numeric. A numeric Time gives one slope, not one level per day.
• Reading the uncorrected aov() table as final. Without a correction the p value is too small whenever sphericity fails.
• Dropping subjects without saying so. State how many subjects the complete-case rule removed.
• Treating a non-significant Mauchly test as proof of sphericity. It means a violation was not detected.
• Forgetting the between-subject factor. ChickWeight records four diets, and this analysis sets them aside. Whether diet changes growth needs a model with diet and day together, which neither run here is.

Rerun the example

The reproduce script beneath the article recomputes the repeated-measures ANOVA from the bundled data and returns the recorded numbers. The companion mixed-model run uses the same four days. The analysis here covers one outcome, four time points and no diet, and describes mean weight across days in this experiment. It makes no claim beyond growth in these chicks, and it is research statistics, not guidance for any individual.

# Results

- Test: Repeated-measures ANOVA (Greenhouse-Geisser)
- n: 196
- statistic: 571.7
- p-value: < 0.001
- effect size: generalized eta-squared = 0.8483 (95% CI 0.8066 to 0.8777)
  (two-sided CI at the run's confidence level, from the uncorrected degrees of freedom)
- seed: 20261111

## Assumption checks
- sphericity (Greenhouse-Geisser correction applied): passed=NA (F, df and p use the Greenhouse-Geisser correction, so a sphericity violation is already accounted for)

## Term tests
| Term | Test | Statistic | df | Residual df | p-value | Partial eta-squared (95% CI) |
|---|---|---|---|---|---|---|
| Time | F | 571.7 | 1.238 | 59.42 | < 0.001 | 0.9225 (0.901 to 0.9377) |
- effect-size interval: two-sided CI at the run's confidence level, from the uncorrected degrees of freedom

## Power
- Power is not calculated for this test.

# Celeus - reproduce.R  (run from this package folder; bundled data in data/)
# Test: Repeated-measures ANOVA (Greenhouse-Geisser)  |  engine 0.53.2  |  catalog c40-dc08
# Research context: not provided for this run.
# Reproduces the recorded result deterministically.

set.seed(20261111)
.data_path <- "data/repeated-measures-anova-in-r.csv"
if (!file.exists(.data_path)) stop("Run reproduce.R from the package folder (the bundled data is in data/).")
data <- readr::read_csv(.data_path, show_col_types = FALSE)  # needs the readr package
clean_cells <- function (df) 
{
    ctrl_collapse <- "[\r\n\t\v\f]+"
    zero_width_chars <- intToUtf8(c(8203L, 65279L), multiple = TRUE)
    uni_space_chars <- intToUtf8(c(160L, 8239L, 8199L), multiple = TRUE)
    clean_text <- function(x) {
        out <- x
        for (ch in zero_width_chars) out <- gsub(ch, "", out, 
            fixed = TRUE)
        out <- gsub(ctrl_collapse, " ", out)
        for (ch in uni_space_chars) out <- gsub(ch, " ", out, 
            fixed = TRUE)
        trimws(gsub(" +", " ", out))
    }
    changed_cols <- character(0)
    cells_cleaned <- 0L
    retyped <- character(0)
    for (j in seq_along(df)) {
        col <- df[[j]]
        if (!is.character(col)) 
            next
        cleaned <- clean_text(col)
        n_changed <- sum(!is.na(col) & !is.na(cleaned) & col != 
            cleaned, na.rm = TRUE)
        if (n_changed > 0L) {
            changed_cols <- c(changed_cols, names(df)[j])
            cells_cleaned <- cells_cleaned + n_changed
        }
        df[[j]] <- cleaned
        non_blank <- cleaned[!is.na(cleaned) & nzchar(cleaned)]
        if (length(non_blank) && all(grepl("^[+-]?(0|[1-9][0-9]*)([.][0-9]+)?$", 
            non_blank))) {
            parsed <- suppressWarnings(as.numeric(cleaned))
            was_blank <- is.na(cleaned) | !nzchar(cleaned)
            lost <- (is.na(parsed) | is.infinite(parsed)) & !was_blank
            if (!any(lost)) {
                finite <- parsed[!is.na(parsed)]
                all_whole <- all(finite == floor(finite))
                in_range <- all(abs(finite) <= .Machine$integer.max)
                df[[j]] <- if (all_whole && in_range) 
                  as.integer(parsed)
                else parsed
                retyped <- c(retyped, names(df)[j])
            }
        }
    }
    names(df) <- clean_text(names(df))
    if (length(changed_cols)) {
        attr(df, "ingest_clean_columns") <- changed_cols
        attr(df, "ingest_clean_cells_cleaned") <- cells_cleaned
    }
    if (length(retyped)) 
        attr(df, "ingest_clean_retyped") <- retyped
    df
}
data <- clean_cells(data)

# Analysis population filter (recorded in the spec): Time in 0, 2, 4, 6
data <- data[which(data[["Time"]] %in% c("0", "2", "4", "6")), , drop = FALSE]
# Complete-case analysis: drop rows missing any analysis value.
.acols <- c("weight", "Time", "Chick")
data <- data[stats::complete.cases(data[, .acols, drop = FALSE]), , drop = FALSE]
result <- { .rm_d <- stats::setNames(data.frame(data[["Chick"]], data[["weight"]], data[["Time"]], check.names = FALSE), c(".v1", ".v2", ".v3")); .rm_d[[".v1"]] <- as.factor(.rm_d[[".v1"]]); .rm_d[[".v3"]] <- as.factor(.rm_d[[".v3"]]); afex::aov_ez(id = ".v1", dv = ".v2", within = ".v3", data = .rm_d, anova_table = list(correction = "GG"), fun_aggregate = mean)$anova_table }
print(result)

# Effect size (ges (generalized eta-squared)) and its confidence interval.
effect <- { .es_fit <- { .rm_d <- stats::setNames(data.frame(data[["Chick"]], data[["weight"]], data[["Time"]], check.names = FALSE), c(".v1", ".v2", ".v3")); .rm_d[[".v1"]] <- as.factor(.rm_d[[".v1"]]); .rm_d[[".v3"]] <- as.factor(.rm_d[[".v3"]]); afex::aov_ez(id = ".v1", dv = ".v2", within = ".v3", data = .rm_d, anova_table = list(correction = "GG"), fun_aggregate = mean) }; list(generalized = effectsize::eta_squared(.es_fit, generalized = TRUE, ci = 0.95, alternative = "two.sided"), partial = effectsize::eta_squared(.es_fit, partial = TRUE, ci = 0.95, alternative = "two.sided")) }
print(effect)