08 - ANOVA and Multiple Comparisons

You have four variants of a web page and you want to know whether the mean loading times are really different. The obvious route is to compare the variants in pairs with a t-test: that is six comparisons, each at the 5% level. The problem is that, even if the four means are identical, the probability that at least one of the six comparisons comes out “significant” is about 20%, not 5% (the simulation below gives 19.9%). The analysis of variance (ANOVA) solves the first half of the problem with a single test: it asks whether the means are all equal, with a guaranteed level of 5%.

The second half, working out which ones differ, requires a pairwise comparison corrected for multiple comparisons. This article, the eighth in the Probability and Statistics series, covers both, for the one-factor case, and shows you with simulations where the test breaks down. You need the t-test and the p-value from the previous article.

The idea: two estimates of the same variance

Imagine three fertilizers, with five plots each. The yields change for two reasons: because the fertilizers are different (variability between groups) and because even plots with the same fertilizer give different yields (variability within groups, the noise). ANOVA measures the two variabilities in the same unit and takes their ratio. Let x̄_i be the mean of group i with n_i observations, x̄ the mean of all N observations and k the number of groups.

The sum of squares between groups is SS_B = Σn_i(x̄_i − x̄)², the sum of squares within groups is SS_W = Σ(x_ij − x̄_i)², and the total sum of squares SS_T = Σ(x_ij − x̄)² is always equal to SS_B + SS_W: an algebraic identity, not an approximation.

Dividing each sum by its degrees of freedom (k − 1 for SS_B, N − k for SS_W) gives the mean squares MS_B and MS_W. MS_W estimates the variance σ² of the data whatever the truth about the means; MS_B estimates it only if the means are equal, and if they are not it tends to be larger. So F = MS_B/MS_W is about 1 under H₀ and grows with the differences between the means.

If the data in each group are Normal with the same variance and independent, the exact distribution of F under H₀ is a Fisher F with (k − 1, N − k) degrees of freedom; the procedure and its assumptions are described in the NIST statistics handbook, section 7.4.3.

The program works through a small exercise by hand and verifies it: three fertilizers with five plots, means 23, 28 and 24, chosen so that the arithmetic comes out round. The grand mean is 25, SS_B = 5·(4 + 9 + 1) = 70, SS_W = 10 + 10 + 10 = 30, so MS_B = 70/2 = 35, MS_W = 30/12 = 2.5 and F = 14. The program repeats the calculation, uses stats.f.sf for the p-value and compares with scipy.stats.f_oneway; it also computes three indices of effect size.

# Worked exercise: one-way ANOVA on three groups of five measurements, by hand and with scipy
import numpy as np
from scipy import stats

gruppi = [np.array([23, 25, 21, 24, 22.0]),    # mean 23
          np.array([27, 29, 26, 28, 30.0]),    # mean 28
          np.array([24, 22, 26, 25, 23.0])]    # mean 24
k, N = len(gruppi), sum(len(g) for g in gruppi)
tutti = np.concatenate(gruppi)
gran_media = tutti.mean()

SS_B = sum(len(g) * (g.mean() - gran_media) ** 2 for g in gruppi)     # between groups
SS_W = sum(((g - g.mean()) ** 2).sum() for g in gruppi)               # within groups
SS_T = ((tutti - gran_media) ** 2).sum()                              # total
MS_B, MS_W = SS_B / (k - 1), SS_W / (N - k)
F = MS_B / MS_W
p = stats.f.sf(F, k - 1, N - k)                                       # sf = 1 - cdf, without loss of precision

print(f"means: {[float(g.mean()) for g in gruppi]}, grand mean {float(gran_media)}")
print(f"SS_B = {SS_B:.1f}, SS_W = {SS_W:.1f}, SS_T = {SS_T:.1f}, SS_B + SS_W = {SS_B + SS_W:.1f}")
print(f"MS_B = {MS_B:.1f}, MS_W = {MS_W:.1f}, F = {F:.1f}, degrees of freedom ({k - 1}, {N - k}), p = {p:.5f}")
r = stats.f_oneway(*gruppi)
print(f"scipy.stats.f_oneway: F = {r.statistic:.1f}, p = {r.pvalue:.5f}")
eta2 = SS_B / SS_T
omega2 = (SS_B - (k - 1) * MS_W) / (SS_T + MS_W)
print(f"eta^2 = {eta2:.4f}, omega^2 = {omega2:.4f}, Cohen's f = {np.sqrt(eta2 / (1 - eta2)):.4f}")

# Output:
# means: [23.0, 28.0, 24.0], grand mean 25.0
# SS_B = 70.0, SS_W = 30.0, SS_T = 100.0, SS_B + SS_W = 100.0
# MS_B = 35.0, MS_W = 2.5, F = 14.0, degrees of freedom (2, 12), p = 0.00073
# scipy.stats.f_oneway: F = 14.0, p = 0.00073
# eta^2 = 0.7000, omega^2 = 0.6341, Cohen's f = 1.5275

The numbers coincide with those by hand: SS_B = 70, SS_W = 30, SS_T = 100 = 70 + 30, F = 14 with (2, 12) degrees of freedom and p = 0.00073, identical to f_oneway. A p-value this low says that differences between the means like these are very unlikely if the fertilizers were equivalent, but it does not say which one is different or by how much.

The indices answer the second part: η² = SS_B/SS_T = 0.70 is the share of the total variability explained by the group; ω² = (SS_B − (k − 1)MS_W)/(SS_T + MS_W) = 0.6341 is a less inflated version (η² tends to overestimate with little data); f = √(η²/(1 − η²)) = 1.5275 is the effect expressed as the standard deviation of the means relative to the standard deviation of the data.

Level and power of the F test

That F really follows the F law under H₀ is checked by simulating: 200,000 experiments with four groups of 30 Normal observations with the same mean. Under H₁, the statistic follows a noncentral F with noncentrality parameter λ = Σn_i(μ_i − μ̄)²/σ², and scipy.stats.ncf gives its probabilities: the SciPy documentation defines the distribution. The true means (10, 10.7, 10.7, 11.4 with σ = 2) are chosen to obtain a “medium” effect according to Cohen, f ≈ 0.25. The program compares exact and simulated power, and looks for how much data per group is needed to reach 80%.

# Under H0 the F statistic follows the F(k-1, N-k) law; under H1 a noncentral F: level and power
import numpy as np
from scipy import stats

def statistica_F(X):
    """X has shape (R, k, n): R repetitions, k groups, n observations per group."""
    R, k, n = X.shape
    medie = X.mean(axis=2)
    SS_B = n * ((medie - medie.mean(axis=1, keepdims=True)) ** 2).sum(axis=1)
    SS_W = ((X - medie[:, :, None]) ** 2).sum(axis=(1, 2))
    return (SS_B / (k - 1)) / (SS_W / (k * (n - 1)))

k, n, sigma, R = 4, 30, 2.0, 200_000
N, soglia = k * n, stats.f.ppf(0.95, k - 1, k * n - k)
rng = np.random.default_rng(42)

F0 = statistica_F(rng.normal(10.0, sigma, size=(R, k, n)))             # H0: four equal means
print(f"H0: threshold F(0.95; {k - 1}, {N - k}) = {soglia:.4f}; you reject {(F0 > soglia).mean():.4f} (expected 0.05)")
for q in (0.50, 0.90, 0.99):
    sotto = (F0 <= stats.f.ppf(q, k - 1, N - k)).mean()
    print(f"    share of F below the {q:.2f} quantile of the F law: {sotto:.4f} +/- {np.sqrt(q * (1 - q) / R):.4f}")

mu = np.array([10.0, 10.7, 10.7, 11.4])                                  # H1: four different means
lam = n * ((mu - mu.mean()) ** 2).sum() / sigma ** 2                     # noncentrality parameter
F1 = statistica_F(rng.normal(mu[None, :, None], sigma, size=(R, k, n)))
esatta = stats.ncf.sf(soglia, k - 1, N - k, lam)
sim = (F1 > soglia).mean()
print(f"H1: lambda = {lam:.3f}, Cohen's f = {np.sqrt(lam / N):.4f}")
print(f"exact power (noncentral F) = {esatta:.4f}; simulated = {sim:.4f} +/- {np.sqrt(sim * (1 - sim) / R):.4f}")

def potenza(n_gruppo, f):
    """Exact power of F with k groups of n_gruppo observations and Cohen's effect f."""
    Ntot = k * n_gruppo
    return stats.ncf.sf(stats.f.ppf(0.95, k - 1, Ntot - k), k - 1, Ntot - k, Ntot * f ** 2)

print(f"f = 0.25, {n} per group: exact power = {potenza(n, 0.25):.4f}")
n80 = next(m for m in range(3, 500) if potenza(m, 0.25) >= 0.80)
print(f"f = 0.25: first n per group with power >= 0.80 = {n80} (power {potenza(n80, 0.25):.4f})")

# Output:
# H0: threshold F(0.95; 3, 116) = 2.6828; you reject 0.0498 (expected 0.05)
#     share of F below the 0.50 quantile of the F law: 0.5003 +/- 0.0011
#     share of F below the 0.90 quantile of the F law: 0.9002 +/- 0.0007
#     share of F below the 0.99 quantile of the F law: 0.9901 +/- 0.0002
# H1: lambda = 7.350, Cohen's f = 0.2475
# exact power (noncentral F) = 0.5968; simulated = 0.5952 +/- 0.0011
# f = 0.25, 30 per group: exact power = 0.6065
# f = 0.25: first n per group with power >= 0.80 = 45 (power 0.8040)

Under H₀ the threshold is 2.6828, you reject 4.98% of the time, and the fractions of F below the 50%, 90% and 99% quantiles of the F law are 0.5003, 0.9002 and 0.9901, all within the printed standard error: the theoretical distribution holds.

Under H₁ the noncentrality parameter is 7.350 (f = 0.2475) and the exact power is 0.5968 against 0.5952 simulated, with a standard error of 0.0011. The result that weighs most for whoever designs the experiment is the last: with 30 observations per group and a medium effect the power is 0.6065, that is, you will fail to see the effect four times out of ten, and for 80% you need 45 observations per group.

Which means differ: Bonferroni and Tukey

The F test answers “at least one mean differs”. To find out which, you compare the pairs, controlling the familywise error (FWER): the probability of having at least one false positive among all the comparisons. With k = 4 there are k(k − 1)/2 = 6 pairs. Bonferroni uses the level 0.05/6 for each: simple and valid with any test, but conservative.

Tukey HSD instead uses the distribution of the studentized range (the maximum difference between means in units of the standard error) and compares each pair with a common threshold. The NIST handbook (section 7.4.7.1) states that with equal sizes the confidence coefficient of the whole set is exactly 1 − α, and that with different sizes the method is conservative.

The program implements the three thresholds on the difference between means, using for all of them the same common estimate MS_W (so the comparison is homogeneous), and checks the Tukey calculation by comparing the p-values with scipy.stats.tukey_hsd on one sample. Then it simulates 100,000 experiments with four groups of 25: first with all means equal (to measure the familywise error), then with the fourth group shifted by 1.5 (0.75 standard deviations) to measure the probability of discovering the pair (1, 4).

# Pairwise comparisons among k = 4 groups: familywise error and power of uncorrected t, Bonferroni, Tukey
import numpy as np
from scipy import stats

k, n, sigma, R = 4, 25, 2.0, 100_000
N, m = k * n, k * (k - 1) // 2                     # m = 6 comparisons
rng = np.random.default_rng(50)
i, j = np.triu_indices(k, 1)                       # the 6 pairs (i, j)

t_unc = stats.t.ppf(1 - 0.05 / 2, N - k)           # uncorrected t threshold (alpha = 0.05 per comparison)
t_bon = stats.t.ppf(1 - 0.05 / (2 * m), N - k)     # Bonferroni: alpha / 6 per comparison
q_tuk = stats.studentized_range.ppf(0.95, k, N - k)
# Thresholds on the difference of means; MS_W is the common variance estimate
def rifiuti(medie_vere):
    X = rng.normal(np.asarray(medie_vere)[None, :, None], sigma, size=(R, k, n))
    medie = X.mean(axis=2)
    MS_W = ((X - medie[:, :, None]) ** 2).sum(axis=(1, 2)) / (N - k)
    diff = np.abs(medie[:, i] - medie[:, j])                       # shape (R, 6)
    se = np.sqrt(MS_W * 2 / n)[:, None]
    return diff > t_unc * se, diff > t_bon * se, diff > q_tuk * se / np.sqrt(2)

# Check of the Tukey implementation against scipy on one sample
x = rng.normal(0, sigma, size=(k, n))
mio = np.abs(x.mean(axis=1)[i] - x.mean(axis=1)[j]) / np.sqrt(((x - x.mean(axis=1, keepdims=True)) ** 2).sum() / (N - k) / n)
p_mio = stats.studentized_range.sf(mio, k, N - k)
print("check Tukey p (mine, scipy):", np.round(p_mio, 5), np.round(stats.tukey_hsd(*x).pvalue[i, j], 5))

print(f"thresholds: t = {t_unc:.3f}, Bonferroni t = {t_bon:.3f}, Tukey q/sqrt(2) = {q_tuk / np.sqrt(2):.3f}")
unc, bon, tuk = rifiuti([10, 10, 10, 10])                          # H0: all means equal
print("H0 true, at least one rejection among the 6 comparisons (familywise error):")
for nome, r in (("uncorrected", unc), ("Bonferroni", bon), ("Tukey", tuk)):
    print(f"  {nome:17s} {r.any(axis=1).mean():.4f}")
print(f"  reference: 1 - 0.95^6 = {1 - 0.95 ** m:.4f} would hold for 6 independent comparisons; Tukey is 0.05 by construction (balanced groups)")

unc, bon, tuk = rifiuti([10, 10, 10, 11.5])                        # H1: only group 4 is different (1.5 = 0.75 sigma)
coppia = [(a, b) for a, b in zip(i, j)].index((0, 3))
print("H1: group 4 has mean +1.5; fraction of samples in which the pair (1, 4) is significant:")
for nome, r in (("uncorrected", unc), ("Bonferroni", bon), ("Tukey", tuk)):
    print(f"  {nome:17s} {r[:, coppia].mean():.4f}")

# Output:
# check Tukey p (mine, scipy): [0.33901 0.15543 0.28332 0.97242 0.99951 0.98823] [0.33901 0.15543 0.28332 0.97242 0.99951 0.98823]
# thresholds: t = 1.985, Bonferroni t = 2.694, Tukey q/sqrt(2) = 2.615
# H0 true, at least one rejection among the 6 comparisons (familywise error):
#   uncorrected       0.1992
#   Bonferroni        0.0405
#   Tukey             0.0495
#   reference: 1 - 0.95^6 = 0.2649 would hold for 6 independent comparisons; Tukey is 0.05 by construction (balanced groups)
# H1: group 4 has mean +1.5; fraction of samples in which the pair (1, 4) is significant:
#   uncorrected       0.7460
#   Bonferroni        0.4853
#   Tukey             0.5174

With all means equal, the uncorrected pairwise t-test gives at least one false positive in 19.92% of the experiments. The reference 1 − 0.95⁶ = 0.2649 would hold for six independent comparisons, but the pairs share groups and the real rate is lower. Bonferroni is at 4.05%, Tukey at 4.95%, very close to the 5% it guarantees by construction with balanced groups.

The price of the correction shows under H₁: the pair with a true difference is discovered in 74.6% of cases without correction, in 48.5% with Bonferroni and in 51.7% with Tukey. The 74.6% is an illusory gain, bought with false positives at a rate you cannot declare; between the two honest methods, Tukey is more powerful than Bonferroni because it exploits the structure of the comparisons. If the comparisons are only against a control group, there is Dunnett’s test, available as scipy.stats.dunnett.

When the assumptions do not hold

The F test requires independence between observations, Normality within the groups and equal variances (homoscedasticity); the documentation of f_oneway lists the three and indicates Kruskal-Wallis and the Alexander-Govern test as alternatives, with a little power lost. Of the three, the most insidious is the one about variances, especially with groups of different sizes: as with the two-sample t-test of the previous article, the trouble arises when the smallest group is also the most dispersed. The classical correction is Welch’s ANOVA, which weights each group by the precision of its mean (w_i = n_i/s_i²) and corrects the degrees of freedom.

The program implements Welch’s (1951) formula and checks it indirectly: with two groups the statistic must coincide with the square of Welch’s t and give the same p-value, which the comparison with ttest_ind(equal_var=False) verifies. Then it simulates three groups with equal means in three arrangements of sizes and standard deviations, and measures how often they reject the classic F, Welch’s F and Kruskal-Wallis. With 10,000 repetitions per case the standard error is about 0.002.

# Different variances and sizes: classic F, Welch's F and Kruskal-Wallis, under H0 (equal means)
import numpy as np
from scipy import stats

def welch_F(X_list):
    """Welch's ANOVA (1951) for groups of different sizes: returns statistic and p-value."""
    k = len(X_list)
    n = np.array([len(x) for x in X_list], dtype=float)
    medie = np.array([x.mean() for x in X_list])
    varianze = np.array([x.var(ddof=1) for x in X_list])
    w = n / varianze                                   # weight = precision of each group's mean
    media_pesata = (w * medie).sum() / w.sum()
    lam = (((1 - w / w.sum()) ** 2) / (n - 1)).sum()
    F = ((w * (medie - media_pesata) ** 2).sum() / (k - 1)) / (1 + 2 * (k - 2) * lam / (k ** 2 - 1))
    df2 = (k ** 2 - 1) / (3 * lam)
    return F, stats.f.sf(F, k - 1, df2)

rng = np.random.default_rng(15)
# check: with two groups Welch's ANOVA coincides with Welch's t-test (F = t^2, same p-value)
a, b = rng.normal(2, 3, 10), rng.normal(0, 1, 40)    # different means, so F is not close to zero
F2, p2 = welch_F([a, b])
tw = stats.ttest_ind(a, b, equal_var=False)
print(f"two groups: Welch F {F2:.4f} versus t^2 {tw.statistic ** 2:.4f}; p {p2:.6f} versus {tw.pvalue:.6f}")

R = 10_000
casi = {"small group = larger variance": ((10, 20, 40), (3.0, 2.0, 1.0)),
        "small group = smaller variance": ((10, 20, 40), (1.0, 2.0, 3.0)),
        "equal sizes, different variances": ((25, 25, 25), (3.0, 2.0, 1.0))}
print("case (3 groups, equal means)             classic F  Welch F  Kruskal-Wallis")
for nome, (ns, sig) in casi.items():
    r_f = r_w = r_k = 0
    for _ in range(R):
        g = [rng.normal(0, s, m) for m, s in zip(ns, sig)]
        r_f += stats.f_oneway(*g).pvalue < 0.05
        r_w += welch_F(g)[1] < 0.05
        r_k += stats.kruskal(*g).pvalue < 0.05
    print(f"{nome:38s} {r_f / R:9.4f}  {r_w / R:7.4f}  {r_k / R:13.4f}")
print(f"standard error with R = {R}: about {np.sqrt(0.05 * 0.95 / R):.4f}")

# Output:
# two groups: Welch F 3.0844 versus t^2 3.0844; p 0.111590 versus 0.111590
# case (3 groups, equal means)             classic F  Welch F  Kruskal-Wallis
# small group = larger variance             0.1931   0.0485         0.1121
# small group = smaller variance            0.0120   0.0485         0.0209
# equal sizes, different variances          0.0618   0.0489         0.0592
# standard error with R = 10000: about 0.0022

The check with two groups is perfect: Welch’s F = t² = 3.0844 and p = 0.111590 in both. When the smallest group (10 observations) has the largest standard deviation (3, against 2 and 1), the classic F rejects a true H₀ 19.3% of the time instead of 5%, Welch 4.85%. In the opposite case the classic one rejects 1.2%, so it loses power without warning.

With equal sizes the error is small (6.2%), and Welch stays at 4.9%. Kruskal-Wallis is not the solution: it compares ranks, not means, and with different dispersions it rejects 11.2% of the time even with equal means. The rule that comes out is the one already seen for the t-test: Welch costs almost nothing with equal variances and protects with different variances.

This article covers only one-factor ANOVA. Two-factor analysis, with interactions, and mixed-effects models for repeated measures on the same subjects use the same decomposition, but require a setup of their own and are not here; for repeated measures keep in mind that the first assumption, independence, does not hold and the one-way F is not applicable.

Four mistakes that are really made

• Comparing all the pairs with uncorrected t-tests: with four groups and equal means, at least one false positive in 19.9% of cases.
• Reading a significant F as “all the means are different”: it only says that they are not all equal.
• Using the classic F with groups of different sizes and variances: the true level is 19.3% when the small group is the most dispersed.
• Designing with “30 per group” and “medium effect” without computing the power: it is 0.61, not 0.80.

The mistakes are all of calibration: the test works as declared only if you respect its assumptions and declare the number of comparisons. Every number in the list comes from a simulation in this article, and none requires more than a few tens of seconds to be redone; this is why it is worth redoing it on your own data, instead of trusting “30 per group is enough”.

You can now decompose the variability into SS_B and SS_W and compute F, p-value and effect by hand and with SciPy, compute the power with the noncentral F and the number of data points per group, compare the pairs with Tukey or Bonferroni while controlling the familywise error, and choose Welch when the variances are not equal. In the next article you move from continuous data to data in categories, with the chi-square test, and in the article on linear regression you will see that ANOVA is a special case of a linear model.

The same series also includes Statistical Inference: Estimators and Properties, Nonparametric Statistics: Wilcoxon and Mann-Whitney and Exploratory Data Analysis (EDA) with Python.