ANOVA 📊
⚙️ Setup¶
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy import stats
# statsmodels = the R/SAS-style stats library (pre-installed on Colab): it computes the
# standard errors, test statistics, intervals, and post-hoc comparisons, so we write less by hand.
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.proportion import proportions_ztest, confint_proportions_2indep, proportion_confint
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from statsmodels.stats.weightstats import DescrStatsW, CompareMeans
from scipy.stats.contingency import association # Cramer's V in one call
IND="#4f46e5"; DEEP="#4338ca"; LIGHT="#818cf8"; INK="#1a2138"; GRID="#e6e9f2"; GREEN="#059669"; RED="#ef4444"; AMBER="#d97706"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
"axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
"axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
rng = np.random.default_rng(78)
from itertools import combinations
# simulate: 5 groups, all identical (H0 true). How often does ANY pairwise t-test fire?
def any_pairwise_hit(k=5,n=30):
gs=[rng.normal(0,1,n) for _ in range(k)]
return any(stats.ttest_ind(a,b).pvalue<0.05 for a,b in combinations(gs,2))
rate=np.mean([any_pairwise_hit() for _ in range(3000)])
print(f"5 identical groups: P(at least one pairwise t-test significant) = {rate:.2f}")
print("Far above 0.05 -> multiple comparisons inflate the family-wise error. ANOVA fixes this.")
5 identical groups: P(at least one pairwise t-test significant) = 0.27 Far above 0.05 -> multiple comparisons inflate the family-wise error. ANOVA fixes this.
Ten independent 5% tests would give a 1−0.95¹⁰ ≈ 40% chance of at least one false alarm. The ten pairwise tests here share groups, so they are correlated, and the simulation finds a lower 27%, still more than five times the 5% each test allows. ANOVA restores a single, honest 5% by asking one global question first: is there any difference at all?
# three groups with genuinely different means
g1=rng.normal(50,8,40); g2=rng.normal(55,8,40); g3=rng.normal(62,8,40)
F,p=stats.f_oneway(g1,g2,g3)
print(f"means: {g1.mean():.1f}, {g2.mean():.1f}, {g3.mean():.1f}")
print(f"F={F:.2f}, p={p:.2e} -> reject H0: not all means are equal")
# contrast: three IDENTICAL groups -> F near 1, p large
h=[rng.normal(50,8,40) for _ in range(3)]; F0,p0=stats.f_oneway(*h)
print(f"identical groups: F={F0:.2f}, p={p0:.2f} -> fail to reject (as expected)")
means: 49.1, 56.2, 61.7 F=24.88, p=9.95e-10 -> reject H0: not all means are equal identical groups: F=1.09, p=0.34 -> fail to reject (as expected)
A large F means the gaps between group means dwarf the within-group noise, evidence the groups really differ. Under the null, between and within variation are the same kind of thing, so F hovers around 1. The p-value comes from the F distribution with (k−1, N−k) degrees of freedom.
# pairwise_tukeyhsd was imported in the setup cell; one call replaces a whole loop
g1=rng.normal(50,8,50); g2=rng.normal(52,8,50); g3=rng.normal(60,8,50)
vals=np.concatenate([g1,g2,g3]); labels=["A"]*50+["B"]*50+["C"]*50
F,p=stats.f_oneway(g1,g2,g3); print(f"ANOVA F={F:.2f}, p={p:.2e}\n")
print(pairwise_tukeyhsd(vals, labels).summary())
ANOVA F=18.24, p=8.37e-08
Multiple Comparison of Means - Tukey HSD, FWER=0.05
====================================================
group1 group2 meandiff p-adj lower upper reject
----------------------------------------------------
A B -0.7997 0.8681 -4.5351 2.9357 False
A C 7.8242 0.0 4.0888 11.5595 True
B C 8.6238 0.0 4.8885 12.3592 True
----------------------------------------------------
One function call, no hand-coded loops. Each row compares a pair: reject=True means that pair differs (its confidence interval excludes 0). Tukey widens each comparison just enough to keep the overall error at 5%, the disciplined way to follow up a significant ANOVA.
A shipper compares three carriers on door-to-door delivery time (anova--carrier_delivery.xlsx), 80 shipments each. Are the average times genuinely different, or is it noise? We explore first, then run ANOVA and Tukey, both from statsmodels.
try: d = pd.read_excel("../../data/anova--carrier_delivery.xlsx", sheet_name="Shipments")
except FileNotFoundError: d = pd.read_excel(BASE+"anova--carrier_delivery.xlsx", sheet_name="Shipments")
# EXPLORE FIRST: size, missing values, and a per-carrier summary
print("shape:", d.shape, "| missing:", d.isna().sum().sum())
print(d.groupby("carrier").delivery_hours.agg(["size","mean","std","median"]).round(2))
shape: (240, 4) | missing: 0
size mean std median
carrier
FastFreight 80 30.39 7.05 31.30
GroundLink 80 35.87 5.92 35.94
RegionalCo 80 37.40 7.59 36.31
# CHECK ASSUMPTIONS before trusting the F-test (independence holds by the study design).
from scipy.stats import levene, shapiro, kruskal
from statsmodels.stats.oneway import anova_oneway
groups = [g.delivery_hours.values for _, g in d.groupby("carrier")]
resid = np.concatenate([g - g.mean() for g in groups]) # pooled within-group residuals
lev = levene(*groups) # equal variance across carriers (robust to non-normality)
sh = shapiro(resid) # normality of the residuals
print(f"Levene equal-variance p = {lev.pvalue:.3f} -> "
+ ("no clear evidence of unequal variances" if lev.pvalue>0.05 else "variances differ; use Welch ANOVA"))
print(f"Shapiro residual-normal p = {sh.pvalue:.3f} -> "
+ ("residuals look normal" if sh.pvalue>0.05 else "non-normal; lean on the large n or use Kruskal-Wallis"))
print(f"Kruskal-Wallis cross-check p = {kruskal(*groups).pvalue:.1e} (distribution-free, should agree with ANOVA)")
w = anova_oneway(d.delivery_hours, d.carrier, use_var="unequal") # Welch ANOVA: no equal-variance assumption
print(f"Welch ANOVA cross-check F = {w.statistic:.2f} on (2, {w.df_denom:.1f}) df, p = {w.pvalue:.1e}")
# a residual QQ plot makes the normality check visual
fig, ax = plt.subplots(figsize=(4.8,4)); stats.probplot(resid, dist="norm", plot=ax)
ax.set_title("Residual QQ plot (points on the line = normal)"); plt.tight_layout(); plt.show()
Levene equal-variance p = 0.095 -> no clear evidence of unequal variances Shapiro residual-normal p = 0.276 -> residuals look normal Kruskal-Wallis cross-check p = 3.1e-09 (distribution-free, should agree with ANOVA) Welch ANOVA cross-check F = 21.46 on (2, 156.2) df, p = 5.8e-09
Before reading the F-test we check its assumptions, because a significant p-value is only trustworthy if the test's conditions hold. Independence is guaranteed by design (each row is one separate shipment). Levene's test (p = 0.10) finds no clear evidence that the carriers' variances differ, and the residual QQ plot plus Shapiro (p = 0.28) show the residuals are roughly normal, so the classic ANOVA is reasonable here. With equal groups of 80 it tolerates modest differences in spread, and Welch's ANOVA, which drops the equal-variance assumption, agrees (F ≈ 21.5, p about 6e-9). When those checks FAIL you do not drop the question, you switch tools: unequal variances call for Welch's ANOVA (statsmodels' anova_oneway(..., use_var="unequal"), run above as a cross-check), and heavy non-normality calls for the rank-based Kruskal-Wallis test (scipy.stats.kruskal). As a safety net we ran Kruskal-Wallis anyway and it agrees with ANOVA (p about 3e-9): the carriers really do differ.
# ANOVA via a formula (like R: hours ~ carrier), then read the table from statsmodels
model = ols("delivery_hours ~ C(carrier)", data=d).fit()
aov = sm.stats.anova_lm(model, typ=2)
aov["mean_sq"] = aov["sum_sq"] / aov["df"] # MS = SS / df: the variance each row represents
aov = aov[["sum_sq", "df", "mean_sq", "F", "PR(>F)"]]
ss_b = aov.loc["C(carrier)", "sum_sq"]; ms_w = aov.loc["Residual", "mean_sq"]; ss_t = aov["sum_sq"].sum()
eta2 = ss_b / ss_t
omega2 = (ss_b - aov.loc["C(carrier)", "df"] * ms_w) / (ss_t + ms_w) # bias-corrected effect size
print(aov.round(3))
print(f"\neta^2 = {eta2:.3f} | omega^2 = {omega2:.3f} (carrier explains about {omega2*100:.0f}% of delivery-time variance, bias-corrected)")
print("\nTukey HSD (which carriers differ?):")
print(pairwise_tukeyhsd(d.delivery_hours, d.carrier).summary())
fig,ax=plt.subplots(figsize=(7,3.3))
ax.boxplot([g.delivery_hours.values for _,g in d.groupby("carrier")], tick_labels=list(d.groupby("carrier").groups))
ax.set_ylabel("delivery hours"); ax.set_title("Delivery time by carrier (one-way ANOVA)")
plt.tight_layout(); plt.show()
sum_sq df mean_sq F PR(>F)
C(carrier) 2173.764 2.0 1086.882 22.914 0.0
Residual 11241.803 237.0 47.434 NaN NaN
eta^2 = 0.162 | omega^2 = 0.154 (carrier explains about 15% of delivery-time variance, bias-corrected)
Tukey HSD (which carriers differ?):
Multiple Comparison of Means - Tukey HSD, FWER=0.05
===========================================================
group1 group2 meandiff p-adj lower upper reject
-----------------------------------------------------------
FastFreight GroundLink 5.4793 0.0 2.9109 8.0476 True
FastFreight RegionalCo 7.0106 0.0 4.4423 9.579 True
GroundLink RegionalCo 1.5314 0.3393 -1.037 4.0997 False
-----------------------------------------------------------
The F statistic is about 22.9 with a p-value near 10⁻⁹, so we reject the hypothesis that all three carriers deliver equally fast; carrier explains about 16% of the variation (η² ≈ 0.16). The Tukey table shows FastFreight (≈ 30.4 h) is significantly faster than both GroundLink (≈ 35.9 h) and RegionalCo (≈ 37.4 h), while GroundLink and RegionalCo are not clearly distinguishable. The omnibus test plus the post-hoc gives the full picture: there is a difference, and it is FastFreight pulling ahead. Read the table left to right: each sum of squares (SS) is a pile of variation, df its degrees of freedom, and MS = SS/df the variance it stands for; the F ratio is the between-carrier MS (about 1087) over the within-carrier MS (about 47), so the between-carrier mean square is about 23 times the within-carrier one, where equal carriers would give a ratio near 1. The bias-corrected omega-squared (0.15) sits just below eta-squared and is the number to report. In publication form: a one-way ANOVA showed delivery time differed by carrier, F(2, 237) = 22.91, p < .001, omega-squared = 0.15; Tukey HSD placed FastFreight (30.4 h) below both GroundLink (35.9 h) and RegionalCo (37.4 h), which did not differ.