Sphericity and Mauchly's test in repeated measures

Sphericity is the assumption that, in a repeated-measures ANOVA, the differences between every pair of levels of the within-subject factor have the same variance. Mauchly’s test is the usual check for it, but the test is a weak gate: the safer habit is to apply the Greenhouse-Geisser correction every time, or to use a mixed model that states its covariance assumption openly.

The sections below explain what the assumption says, what a violation does to the F test, how Mauchly’s test and the epsilon corrections work, and how to read epsilon straight off a reported result. The worked example is a corrected repeated-measures ANOVA on the first four weighings of R’s ChickWeight data. It sits under the guide to repeated-measures ANOVA and linear mixed models, which compares the two methods on the full data.

What sphericity means

Repeated-measures ANOVA compares the means of one outcome measured on the same subjects at several levels of a factor, such as days, doses or conditions. Each subject serves as their own control, so the analysis works with within-subject differences. Sphericity asks that those differences behave alike: the variance of (day 2 minus day 0) should equal the variance of (day 6 minus day 0), of (day 6 minus day 4), and so on for every pair.

With two levels there is only one difference, so the assumption cannot fail. The question starts at three levels.

A stronger condition, compound symmetry, says every level has the same variance and every pair of levels has the same correlation. Compound symmetry is enough for sphericity but not required: Huynh and Feldt (1970) showed that the F ratio has an exact F distribution under the weaker sphericity condition. In practice the distinction matters less than the reason both fail. Measurements taken close together usually agree more than measurements taken far apart, and in anything that grows, subjects spread out over time. Both patterns make some differences far more variable than others. Separately from sphericity, the F test also assumes the within-subject differences are roughly normal.

What a violation does to the F test

When sphericity fails, the F statistic no longer follows the F distribution its degrees of freedom suggest, and the uncorrected test rejects a true null hypothesis more often than the nominal 5%. Box (1954) described the effect with a single number, epsilon, that measures how far the covariance of the repeated measurements departs from sphericity.

Epsilon equals 1 under sphericity. Its smallest possible value is 1 divided by the number of levels minus one: 1/3 for four levels, 1/11 for twelve. The corrections all work the same way. Both degrees of freedom of the F test are multiplied by an estimate of epsilon, and the p value is read from the F distribution with those smaller degrees of freedom. The F statistic itself does not change. Only the reference distribution does.

Mauchly’s test and why it is a weak gate

Mauchly (1940) gave a test of the null hypothesis that the covariance of the within-subject contrasts is spherical. The statistic, W, lies between 0 and 1, with 1 meaning perfect sphericity, and is converted to a chi-square test. A small p value is evidence of a violation.

The common workflow runs Mauchly’s test first, keeps the uncorrected F if the test is not significant, and corrects otherwise. Three problems undermine that gate:

• A non-significant result does not show sphericity holds. The test has little power in small samples, so a real violation can pass undetected.
• The test depends on normality. Blanca and colleagues (2023) summarise earlier simulations in which Mauchly’s test did not control its error rate with non-normal data, and was not sensitive to small departures from sphericity.
• The gate itself is a second test. The final p value then depends on a preliminary decision that the reported F does not account for.

Blanca and colleagues therefore proposed a rule based on the epsilon estimate alone, noting that it may remove the need for the preliminary test. They also found that the uncorrected F became more liberal as the violation grew, while both corrected tests kept the Type I error rate within a liberal robustness criterion. Their simulations used normal data.

Greenhouse-Geisser, Huynh-Feldt and the lower bound

Three estimates of epsilon are in common use.

• Greenhouse-Geisser (1959) estimates epsilon from the sample covariance matrix. The estimate tends to be too low when the true value is near 1, so the corrected test leans conservative.
• Huynh-Feldt (1976) adjusts the Greenhouse-Geisser estimate upwards to remove that bias. The adjusted value can exceed 1, in which case it is set to 1 for the correction. The corrected test leans liberal.
• The lower bound uses the smallest possible epsilon, 1 divided by the levels minus one. With normal data the result is valid whatever the covariance, and usually too conservative to be useful.

For choosing between the first two, the 2023 simulations recommended Greenhouse-Geisser when epsilon is below 0.60 and Huynh-Feldt at or above it, and Greenhouse-Geisser as the general choice because it is the more conservative. The same study advised against the uncorrected F even with epsilon as high as 0.90. Whatever the rule, fix it before seeing the data and name it in the methods.

Worked example: the first four weighings of ChickWeight

ChickWeight ships with R. It comes from an experiment on diet and early growth (Crowder and Hand, 1990): 50 chicks weighed at hatching, every second day to day 20, and again on day 21. The example takes the first four weighings, days 0, 2, 4 and 6, and asks whether mean weight differs across those days. Day is the within-subject factor and the chick is the subject.

The repeated-measures ANOVA was run in R through the Celeus engine, restricted to the four days by a recorded population filter. Classic repeated-measures ANOVA needs every subject at every level, so one chick, missing the day 4 and day 6 weighings, was dropped whole, leaving 196 weighings from 49 chicks (196 divided by 4). The engine applies the Greenhouse-Geisser correction to every repeated-measures ANOVA it runs. Its output and script appear below the article.

The result

Mean weight differed across the four days: F(1.24, 59.42) = 571.7, p < 0.001, with the Greenhouse-Geisser correction. Nobody needed a test to expect chicks to grow. The point of the example is the degrees of freedom.

Reading epsilon off the degrees of freedom

Without a correction the test would use 3 and 144 degrees of freedom: four days minus one, and (49 chicks minus one) times three. The reported values are 1.24 and 59.42. Dividing either corrected value by its uncorrected counterpart gives the Greenhouse-Geisser epsilon: 1.24 / 3 and 59.42 / 144 both come to about 0.41.

An epsilon of 0.41 signals a large violation. The floor for four levels is 1/3, about 0.33, and the estimate sits much nearer the floor than 1. By the 0.60 rule above, Greenhouse-Geisser is the correction to use. The same arithmetic works on any published result: a repeated-measures F with fractional degrees of freedom has been corrected, and their ratio to the uncorrected ones is the epsilon used.

Here the correction does not change the conclusion, because F is enormous. In a study with a modest effect, cutting the degrees of freedom by more than half can move p across 0.05, which is precisely when an uncorrected test misleads.

Effect size

Generalized eta-squared was 0.8483 (95% CI 0.8066 to 0.8777), and partial eta-squared 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 between-chick variation in the denominator and so stays comparable with between-subject designs (Bakeman, 2005). The post on eta squared explains the difference between the two. Both intervals are computed from the uncorrected degrees of freedom, as the engine’s output notes, so they do not carry the sphericity correction.

What the engine reports, and what it does not

The engine does not run Mauchly’s test and reports no epsilon at all, only the corrected degrees of freedom, so it has no Huynh-Feldt figure either. Its assumption line for this run says the F, df and p use the Greenhouse-Geisser correction, so a sphericity violation is already accounted for. The engine also flagged 9 weight values beyond Tukey’s fences, descriptive only, and the dropped chick.

To see Mauchly’s test and both epsilons for the same window, base R is enough (the repeated-measures ANOVA in R post runs the full analysis):

cw <- as.data.frame(subset(ChickWeight, Time %in% c(0, 2, 4, 6)))
wide <- reshape(cw[, c("Chick", "Time", "weight")], idvar = "Chick",
                timevar = "Time", direction = "wide")
wide <- na.omit(wide)  # complete chicks only, as the ANOVA requires
fit <- lm(cbind(weight.0, weight.2, weight.4, weight.6) ~ 1, data = wide)
mauchly.test(fit, X = ~1)                  # Mauchly's W and its p value
anova(fit, X = ~1, test = "Spherical")     # both epsilons and corrected p values

In the second table, the row labelled (Intercept) is the test for day, because the model works on the differences between days. Its G-G column is the Greenhouse-Geisser p value, and the Greenhouse-Geisser epsilon printed above the table should match the 0.41 read off the degrees of freedom.

Why a mixed model sidesteps the test

Repeated-measures ANOVA has one covariance assumption built in and offers a correction when it fails. A linear mixed model turns the question around: you state a covariance pattern as part of the model. Candidate patterns, chosen in advance, are compared on the data with an information criterion, instead of one pattern being screened by a preliminary test.

The choice still matters. A mixed model with only a random intercept per subject implies compound symmetry, which is a special case of sphericity, so its F test for time carries the same assumption as the uncorrected ANOVA. A random slope for time lets subjects grow at different rates, which lets the variance rise over time and the correlation between visits change with their spacing: the kind of pattern that can push epsilon as low as the 0.41 seen here. Celeus’s linear mixed model can add that random slope for a numeric time variable. Covariance structures for the residuals, such as first-order autoregressive or unstructured, are another route, available in R’s nlme and mmrm packages.

A mixed model also keeps the chick this example had to drop, since it uses every observed weighing. The guide to repeated-measures ANOVA and mixed models runs both on all twelve weighings and shows what each one keeps.

How to report it

A results sentence for the example might read:

Mean weight differed across days 0, 2, 4 and 6 (repeated-measures ANOVA with the Greenhouse-Geisser correction, epsilon = 0.41, F(1.24, 59.42) = 571.7, p < 0.001; generalized eta-squared = 0.85, 95% CI 0.81 to 0.88). One chick without a complete set of weighings was excluded.

In the methods, say that the Greenhouse-Geisser correction was applied by default rather than because Mauchly’s test was significant, and keep the fractional degrees of freedom: they show the reader which test was run. Follow-up comparisons between days are separate tests that need their own adjustment; the guide to multiple comparisons covers that step, and the guide to paired tests covers the two-level case where sphericity never arises.

Rerunning the example

The script below the article is the R script exactly as the engine generated it, with the seed and the population filter. On the same package versions it reproduces the corrected F test, its p value and both eta-squared values with their intervals. The Mauchly and Huynh-Feldt lines come only from the base R code above, not from the engine.

The example is about the method. It says nothing about how chicks should be fed or raised.

# 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: 20261031

## 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.41.7  |  catalog c40-dc08
# Research context: not provided for this run.
# Reproduces the recorded result deterministically.

set.seed(20261031)
.data_path <- "data/sphericity-mauchlys-test.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)