Analysis of variance

You have three or more groups and one numeric outcome, and you want to know whether the group means differ by more than chance would produce. Blood pressure under three drug regimens. Birth weight across four clinics. Haemoglobin by type of sickle cell disease. The test for this is the one-way analysis of variance, written ANOVA. The name reads backwards at first, because the outcome you care about is a mean but the tool is variance. That is the idea worth learning here: you decide whether several means differ by comparing two kinds of variation. This lesson assumes you already know the two-sample t test and what a P-value is.

Why not just run a t test on every pair

With three groups you could test A versus B, A versus C, and B versus C: three separate t tests. With more groups the count climbs fast. For k groups the number of pairwise tests is k(k−1)/2, which is 3 tests for 3 groups, 10 for 5 groups, and 45 for 10 groups.

Each test carries its own 5% chance of a false positive. Run many and those chances pile up. If the comparisons were independent and the null were true everywhere, the probability that at least one test fires falsely is:

For m = 3 tests at α = 0.05 that is 1 − 0.95³ = 0.14. For 10 tests it is 1 − 0.95¹⁰ = 0.40. So with five groups you have about a 40% chance of declaring at least one difference that is not real, before you have learned anything. This is the multiple comparisons problem. ANOVA steps around it by asking one global question first: do any of the means differ? One test, one P-value, and the false-positive rate stays at 5%.

Two kinds of variation

Picture the data plotted by group. Two things vary at once. The group means scatter around the grand mean, which is the between-group variation. And individuals scatter around their own group mean, which is the within-group variation, also called residual variation. If the groups are truly identical, the spread of the group means is only the noise you would expect from sampling, so between-group variation is about the same size as within-group variation. If the groups really differ, the means get pushed apart and between-group variation grows larger than within. ANOVA turns that comparison into a single number.

Partitioning the sum of squares

Start with the total variation: the sum of squared deviations of every observation from the grand mean, the same total sum of squares that sits behind the ordinary variance. ANOVA splits it cleanly into two parts.

The total degrees of freedom, n−1, split the same way: k−1 go to the between part and n−k to the within part, where n is the total number of observations and k the number of groups. Dividing each sum of squares by its own degrees of freedom gives a variance per degree of freedom, called the mean square.

The F-statistic and the F-test

The F-statistic is the ratio of the two mean squares.

Read F as a signal-to-noise ratio. The numerator carries the spread between groups, the denominator the spread within groups. If the null hypothesis of equal means holds, both mean squares estimate the same underlying variance, so F sits near 1. If the means differ, the numerator inflates and F climbs above 1. Under the null, F follows an F distribution, which, unlike most distributions you have met, needs a pair of degrees of freedom: k−1 in the numerator and n−k in the denominator. The F-test compares your F against that distribution and returns a one-sided P-value, because only a large F counts as evidence. A large F with a small P means at least one group mean differs from the rest.

Worked example: three antihypertensive drugs

A small trial measures the drop in systolic blood pressure (SBP, in mmHg) after four weeks on each of three drugs, with five patients per drug. The numbers are kept clean so you can follow every step by hand.

1. The grand mean is 150/15 = 10.
2. Within-group SS. Inside every group the deviations are −2, −1, 0, 1, 2. Squaring and summing gives 10 per group, so SS within = 30 on n−k = 12 degrees of freedom. MS within = 30/12 = 2.5.
3. Between-group SS. 5×(10−10)² + 5×(13−10)² + 5×(7−10)² = 0 + 45 + 45 = 90 on k−1 = 2 degrees of freedom. MS between = 90/2 = 45.
4. F. F = 45 / 2.5 = 18 on (2, 12) degrees of freedom. The 0.1% critical value for F(2, 12) is about 12.97, so P < 0.001.

The check sums: SS total = 90 + 30 = 120, which equals the sum of squared deviations of all 15 values from the grand mean. The means really do differ, and Drug B gives the largest SBP reduction.

sbp  <- c(8,9,10,11,12, 11,12,13,14,15, 5,6,7,8,9)
drug <- factor(rep(c("A","B","C"), each = 5))

summary(aov(sbp ~ drug))

fit <- aov(sbp ~ drug)
summary(fit)              # drug: SS=90 df=2 MS=45; resid: SS=30 df=12 MS=2.5; F=18, p<0.001
tapply(sbp, drug, mean)   # 10, 13, 7

import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols
df = pd.DataFrame({
    "sbp":  [8,9,10,11,12, 11,12,13,14,15, 5,6,7,8,9],
    "drug": ["A"]*5 + ["B"]*5 + ["C"]*5,
})

model = ols("sbp ~ C(drug)", data=df).fit()
print(sm.stats.anova_lm(model, typ=2))

model = ols("sbp ~ C(drug)", data=df).fit()
print(sm.stats.anova_lm(model, typ=2))   # C(drug): sum_sq=90 df=2 F=18 PR(>F)≈0.0002
print(df.groupby("drug")["sbp"].mean())  # A=10, B=13, C=7

Worked example: haemoglobin by sickle cell type

Real data from Anionwu and colleagues (1981) recorded steady-state haemoglobin (Hb, in g/dL) for patients with three types of sickle cell disease.

Partitioning the total sum of squares gives SS total = 137.85 on 40 df. Of that, 99.89 (72.5%) is between-group and 37.96 is within-group. The mean squares are MS between = 99.89/2 = 49.94 and MS within = 37.96/38 = 1.00. So F = 49.94 / 1.00 = 49.9 on (2, 38) degrees of freedom, with P < 0.001. There is strong evidence that mean haemoglobin differs by disease type: lowest in Hb SS, intermediate in Hb S/β-thalassaemia, and highest in Hb SC.

ss   <- c(7.2,7.7,8.0,8.1,8.3,8.4,8.4,8.5,8.6,8.7,9.1,9.1,9.1,9.8,10.1,10.3)
sbt  <- c(8.1,9.2,10.0,10.4,10.6,10.9,11.1,11.9,12.0,12.1)
sc   <- c(10.7,11.3,11.5,11.6,11.7,11.8,12.0,12.1,12.3,12.6,12.6,13.3,13.3,13.8,13.9)
hb   <- c(ss, sbt, sc)
type <- factor(rep(c("HbSS","HbSbeta","HbSC"), c(16,10,15)))

summary(aov(hb ~ type))

summary(aov(hb ~ type))    # type: SS=99.9 df=2 MS=49.9; resid: SS=38.0 df=38 MS=1.0; F=49.9, p<0.001
tapply(hb, type, mean)     # HbSS 8.71, HbSbeta 10.63, HbSC 12.30

import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols
ss  = [7.2,7.7,8.0,8.1,8.3,8.4,8.4,8.5,8.6,8.7,9.1,9.1,9.1,9.8,10.1,10.3]
sbt = [8.1,9.2,10.0,10.4,10.6,10.9,11.1,11.9,12.0,12.1]
sc  = [10.7,11.3,11.5,11.6,11.7,11.8,12.0,12.1,12.3,12.6,12.6,13.3,13.3,13.8,13.9]
df = pd.DataFrame({"hb": ss+sbt+sc,
                   "type": ["HbSS"]*16 + ["HbSbeta"]*10 + ["HbSC"]*15})

model = ols("hb ~ C(type)", data=df).fit()
print(sm.stats.anova_lm(model, typ=2))

model = ols("hb ~ C(type)", data=df).fit()
print(sm.stats.anova_lm(model, typ=2))   # C(type): sum_sq≈99.9 df=2 F≈49.9 PR(>F)<0.001
print(df.groupby("type")["hb"].mean())

A boxplot makes the between-group separation visible before any test. Draw it live.

ss   <- c(7.2,7.7,8.0,8.1,8.3,8.4,8.4,8.5,8.6,8.7,9.1,9.1,9.1,9.8,10.1,10.3)
sbt  <- c(8.1,9.2,10.0,10.4,10.6,10.9,11.1,11.9,12.0,12.1)
sc   <- c(10.7,11.3,11.5,11.6,11.7,11.8,12.0,12.1,12.3,12.6,12.6,13.3,13.3,13.8,13.9)
hb   <- c(ss, sbt, sc)
type <- factor(rep(c("HbSS","HbSbeta","HbSC"), c(16,10,15)),
               levels = c("HbSS","HbSbeta","HbSC"))

boxplot(hb ~ type)

boxplot(hb ~ type, col = "lightsteelblue",
        ylab = "Haemoglobin (g/dL)", xlab = "Sickle cell type",
        main = "Hb by disease type")
abline(h = mean(hb), lty = 2)   # grand mean reference line

import numpy as np
import matplotlib.pyplot as plt
ss  = [7.2,7.7,8.0,8.1,8.3,8.4,8.4,8.5,8.6,8.7,9.1,9.1,9.1,9.8,10.1,10.3]
sbt = [8.1,9.2,10.0,10.4,10.6,10.9,11.1,11.9,12.0,12.1]
sc  = [10.7,11.3,11.5,11.6,11.7,11.8,12.0,12.1,12.3,12.6,12.6,13.3,13.3,13.8,13.9]

plt.boxplot([ss, sbt, sc], labels=["HbSS","HbSbeta","HbSC"])
plt.show()

grand = np.mean(ss + sbt + sc)
plt.boxplot([ss, sbt, sc], labels=["HbSS","HbSbeta","HbSC"])
plt.axhline(grand, ls="--", color="gray")   # grand mean reference
plt.ylabel("Haemoglobin (g/dL)"); plt.xlabel("Sickle cell type")
plt.title("Hb by disease type"); plt.show()

Assumptions of one-way ANOVA

The F-test rests on three conditions.

• Normality. The outcome is roughly normally distributed within each group. ANOVA is fairly robust here, so moderate departures can be ignored.
• Equal variances. The spread is the same in every group, also called homogeneity of variance. This one matters: badly unequal variances distort the F-test and its P-value.
• Independence. Observations are independent. Repeated measurements on the same patient break this and need a different model.

One detail ties ANOVA back to what you already know. With exactly two groups, one-way ANOVA and the two-sample t test give the same P-value, and F equals t squared. ANOVA is the generalisation of the t test to more than two groups.

After a significant F: post-hoc comparisons

A significant F is a starting gun, not a finish line. It says at least one mean differs, but not which ones. To locate the differences you run post-hoc comparisons, which test the pairs while still controlling the family-wise error rate you set out to protect. Tukey's honestly significant difference (HSD) compares every pair and widens the intervals just enough to hold the overall error at 5%. Bonferroni does the same job more bluntly by dividing α by the number of comparisons. Run these only after the overall F is significant, never as a fishing trip.

sbp  <- c(8,9,10,11,12, 11,12,13,14,15, 5,6,7,8,9)
drug <- factor(rep(c("A","B","C"), each = 5))

TukeyHSD(aov(sbp ~ drug))

TukeyHSD(aov(sbp ~ drug))   # B-A diff=3, C-A diff=-3, C-B diff=-6; adjusted p-values for each pair

import pandas as pd
from statsmodels.stats.multicomp import pairwise_tukeyhsd
df = pd.DataFrame({
    "sbp":  [8,9,10,11,12, 11,12,13,14,15, 5,6,7,8,9],
    "drug": ["A"]*5 + ["B"]*5 + ["C"]*5,
})

print(pairwise_tukeyhsd(df["sbp"], df["drug"]))

print(pairwise_tukeyhsd(df["sbp"], df["drug"]))   # meandiff and reject flag for A-B, A-C, B-C

Common mistakes

• Running many t tests instead of one ANOVA. Ten pairwise tests across five groups push the false-positive rate near 40%. Do the global F-test first.
• Stopping at the F-test. A significant F says some mean differs, not which. Without post-hoc comparisons you cannot name the groups that differ.
• Ignoring unequal variances. ANOVA tolerates non-normality but not badly unequal spreads. Check the group variances; transform or use Welch ANOVA when they differ.
• Treating an ordered exposure as unordered. If the groups have a natural order, like dose or stage, a one-way ANOVA across labels ignores the ordering. Consider a linear trend instead.
• Forgetting independence. Repeated measures on the same patient violate the independence assumption and inflate significance. Use a repeated-measures or mixed model.

Tips

• Plot a boxplot by group before testing. It shows the between-group separation and flags unequal spreads in one picture.
• Report F with both degrees of freedom, for example F(2, 38) = 49.9, P < 0.001. The pair of df is part of the result.
• Check that SS between plus SS within equals SS total. If it does not, the partition or the grouping is wrong.
• For two groups, confirm that F equals t squared. It is a quick sanity check that the ANOVA is set up correctly.
• Decide the post-hoc method before seeing the data, and run it only when the global F is significant.