3.8 Checking assumptions

Every test assumes something about the data. Some assumptions matter a lot, some hardly at all. Look at the data first; use tests of assumptions as a second opinion, not as gatekeepers.

import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
from matplotlib.axes import Axes
from scipy import stats
from sklearn import datasets

from statly import style

style.use()
rng = np.random.default_rng(1)

engel = (
    sm.datasets.engel.load_pandas().data
)  # income and food spending, Belgian households (1857)
income = engel.income.to_numpy()
wine = datasets.load_wine(as_frame=True).frame  # 178 Italian wines from three cultivars
CULTIVAR_COLORS = [style.BLUE, style.GREEN, style.CORAL]

3.8.1 Which assumptions matter?

The t-test, ANOVA and regression make three kinds of assumption: the units are independent, the groups have equal spreads (for Student’s t-test and classical ANOVA), and the values, or the noise, are normal. What happens when each is broken? Simulate a world where the null hypothesis is true and count the false alarms: a valid test at the 5% level should raise one in twenty.

STUDIES = 4_000


def draw(kind: str, n: int, sd: float = 1, clusters: int = 6) -> np.ndarray:
    """STUDIES simulated samples of n values: normal, skewed, or made of a few clusters."""
    if kind == "skewed":
        return rng.exponential(1, (STUDIES, n)) - 1  # exponential, centred on 0
    if kind == "clustered":  # each cluster shares a random shift
        shifts = np.repeat(rng.normal(0, 1, (STUDIES, clusters)), n // clusters, axis=1)
        return shifts + rng.normal(0, 0.5, (STUDIES, n))
    return rng.normal(0, sd, (STUDIES, n))


def false_alarms(a: tuple, b: tuple, equal_var: bool) -> float:
    """Share of t-tests with p < 0.05 on simulated groups with no true difference."""
    p = stats.ttest_ind(draw(*a), draw(*b), axis=1, equal_var=equal_var).pvalue
    return float((p < 0.05).mean())


CASES = [  # (name, first group, second group, assume equal spreads?)
    ("all fine (normal, 30 + 30)", ("normal", 30), ("normal", 30), True),
    ("skewed values, 30 + 30", ("skewed", 30), ("skewed", 30), True),
    ("skewed values, 8 + 8", ("skewed", 8), ("skewed", 8), True),
    ("unequal spreads, 10 + 40, Student", ("normal", 10, 3), ("normal", 40), True),
    ("unequal spreads, 10 + 40, Welch", ("normal", 10, 3), ("normal", 40), False),
    ("clustered (6 clusters of 5), Welch", ("clustered", 30), ("clustered", 30), False),
]
rates = [false_alarms(a, b, equal) for _, a, b, equal in CASES]


def plot_rates(ax: Axes) -> None:
    """One bar per scenario: its false-alarm rate."""
    names = [name for name, *_ in CASES]
    colors = [style.CORAL if r > 0.08 else style.BLUE for r in rates]
    ax.barh(names, rates, color=colors, alpha=0.7)
    for i, r in enumerate(rates):
        ax.text(r + 0.005, i, f"{r:.1%}", va="center", fontsize=9)
    ax.axvline(0.05, color=style.INK, ls="--", lw=1)
    ax.xaxis.set_major_formatter(style.PERCENT)
    ax.invert_yaxis()
    ax.set(xlabel="false alarms at the 5% level", xlim=(0, max(rates) * 1.2))


fig, ax = plt.subplots(figsize=(8, 3.2))
plot_rates(ax)
plt.show()

Figure 3.8.1: Share of false alarms of two-group tests at the 5% level, in 4,000 simulated studies where the groups truly do not differ. Skewed data barely matter with 30 per group. Unequal spreads with unequal group sizes break Student’s t-test but not Welch’s. Dependence, here 30 values from only 6 clusters per group, breaks every test badly.

The order of importance is clear: independence first, then equal spreads (only for tests that assume them, and mostly when group sizes differ), and normality last, at least for averages of more than a handful of values, thanks to the central limit theorem (Chapter 1.1). With very small samples and strong skew, normality does start to matter, and so do outliers.

3.8.2 Normality

The best check of a distribution’s shape is a picture: a histogram, and a QQ plot (Chapter 2.3). Tests exist too. The Shapiro-Wilk test asks whether the sorted values line up with the values a normal distribution would give, and is the most powerful common test of normality. But a test of normality answers the wrong question at both ends:

def plot_qq(ax: Axes, values: np.ndarray, title: str, color: str) -> None:
    """A QQ plot against a fitted normal, with the Shapiro-Wilk p-value."""
    values = np.sort(values)
    expected = stats.norm.ppf(
        (np.arange(1, len(values) + 1) - 0.5) / len(values), values.mean(), values.std()
    )
    ends = [expected.min(), expected.max()]
    ax.plot(ends, ends, color=style.INK, ls="--", lw=1)
    ax.plot(expected, values, "o", ms=4 if len(values) < 100 else 1.5, color=color, alpha=0.7)
    p = stats.shapiro(values).pvalue
    style.note(
        ax,
        f"Shapiro-Wilk p = {p:.2f}" if p > 0.001 else "Shapiro-Wilk p < 0.001",
        x=0.03,
        y=0.98,
        color=style.INK,
    )
    ax.set(title=title, xlabel="expected if normal", ylabel="data")


small_skewed = np.random.default_rng(58).lognormal(0, 0.9, 12)
large_nearly = np.random.default_rng(11).standard_t(20, 5_000)
fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_qq(axes[0], small_skewed, "12 skewed values", style.CORAL)
plot_qq(axes[1], large_nearly, "5,000 almost normal values", style.BLUE)
fig.tight_layout()
plt.show()

Figure 3.8.2: QQ plots against the normal, with the Shapiro-Wilk p-value. Left: 12 values from a strongly skewed (log-normal) population pass the test (p = 0.14); about one such sample in four does. With few data the test often cannot see the skew, though that is exactly when skew matters. Right: 5,000 values from a t distribution with 20 degrees of freedom, almost normal and harmless for any test, fail it: with lots of data it flags trivial differences.

How does the test decide? Its statistic W is, roughly, the squared correlation between the sorted data and the values a normal distribution would give: the straightness of the QQ plot. W = 1 is a perfectly straight line. Its null distribution comes from samples that really are normal, and its width depends strongly on n:

null_small = np.array([stats.shapiro(rng.normal(size=12)).statistic for _ in range(10_000)])
null_large = np.array([stats.shapiro(rng.normal(size=5_000)).statistic for _ in range(1_000)])
fig, axes = style.subplots(2, height=3.1, width=10.5)
style.null_panel(
    axes[0],
    stats.shapiro(small_skewed).statistic,
    samples=null_small,
    tail="left",
    xlabel="Shapiro-Wilk W",
)
axes[0].set_title("12 values")
style.null_panel(
    axes[1],
    stats.shapiro(large_nearly).statistic,
    samples=null_large,
    tail="left",
    xlabel="Shapiro-Wilk W",
)
axes[1].set_title("5,000 values")
fig.tight_layout()
plt.show()

Figure 3.8.3: The Shapiro-Wilk statistic W for 10,000 samples drawn from a normal distribution (grey), with the observed W of the two samples above (dashed). Small W means a crooked QQ plot, so the p-value is the left tail. Left: with 12 values, even truly normal samples often give a W as low as the skewed sample’s, so the test cannot tell them apart. Right: with 5,000 values, normal samples give W within a hair of 1, and the slightest bend, here the t distribution’s slightly heavy tails, lies outside that range.

So do not choose a test by running a normality test first. Look at the plots, think about where the data come from (counts? amounts that cannot go below zero? a few extreme cases?), and remember that the t-test and ANOVA care about the shape of the averages, which become normal fast. Two other normality tests you will meet: Kolmogorov-Smirnov with Lilliefors’ correction (Chapter 2.3), and Anderson-Darling, which looks harder at the tails.

3.8.3 Equal spreads

Student’s t-test and classical ANOVA assume every group has the same spread. Levene’s test checks this by running an ANOVA on each value’s distance from its group’s centre (the median, in the robust version by Brown and Forsythe): if one group’s distances are larger, its spread is larger. Bartlett’s test does the same job but is thrown off by non-normal data.

groups = [wine.color_intensity[wine.target == k].to_numpy() for k in range(3)]
levene = stats.levene(*groups)


def plot_spread(ax: Axes, values: list[np.ndarray], title: str, ylabel: str) -> None:
    """Values per cultivar as dots with group means."""
    for k, (group, color) in enumerate(zip(values, CULTIVAR_COLORS)):
        ax.scatter(k + rng.uniform(-0.15, 0.15, len(group)), group, s=9, color=color, alpha=0.6)
        ax.plot([k - 0.3, k + 0.3], [group.mean()] * 2, color=color, lw=2.5)
    ax.set_xticks(range(3), ["1", "2", "3"])
    ax.set(title=title, xlabel="cultivar", ylabel=ylabel)


fig, axes = style.subplots(3, height=3.1, width=11)
plot_spread(axes[0], groups, "Colour intensity", "colour intensity")
plot_spread(
    axes[1],
    [abs(g - np.median(g)) for g in groups],
    "Distance from the group median",
    "distance",
)
n_wines = sum(len(g) for g in groups)
style.null_panel(
    axes[2],
    levene.statistic,
    dist=stats.f(2, n_wines - 3),
    tail="right",
    xlabel="F of the distances",
)
axes[2].set_title("If the spreads were equal")
fig.tight_layout()
plt.show()

Figure 3.8.4: Colour intensity of the wines by cultivar: the spreads differ a lot, cultivar 3 varying three times as much as cultivar 2. Left: the values. Middle: each wine’s distance from its group’s median, which is what Levene’s test compares. Right: Levene’s test is an ANOVA on those distances, so its statistic is an F; the observed F lies far beyond the F distribution with 2 and 175 degrees of freedom that equal spreads would give.

When spreads differ, use the versions that do not assume equality, Welch’s t-test and Welch’s ANOVA, rather than testing first and then choosing. The spreads themselves can be the finding: a treatment that makes outcomes more variable matters even if the average does not move.

3.8.4 Independence

Independence is the assumption that matters most and is hardest to check from the data alone. It fails when the units come in clusters (pupils in classes, patients in hospitals, repeated measurements of one person) or in sequence (days, years). Values within a cluster are more alike than values from different clusters, so the data hold less information than their count suggests, and every standard error computed as if they were independent is too small.

The fix is in the design and the analysis, not in a test: analyse the clusters (for example their averages), use methods built for clustered data (mixed models, cluster-robust standard errors, Chapter 4.3), or for sequences, time series methods (Chapter 9.1). The Durbin-Watson test checks whether the residuals of a regression correlate with their neighbours in sequence.

3.8.5 Transformations

When values are skewed because effects multiply, a logarithm often makes them close to normal (Chapter 2.1), and turns comparisons of differences into comparisons of ratios. A square root is gentler and suits counts. The Box-Cox family includes both and everything in between:

$$x^{(\lambda)} = \begin{cases} \dfrac{x^\lambda - 1}{\lambda} & \lambda \ne 0 \\[4pt] \log x & \lambda = 0 \end{cases}$$

λ = 1 leaves the data as they are, ½ is roughly a square root, 0 the log, −1 roughly 1/x. The λ that makes the data most normal can be found by maximum likelihood (Chapter 2.3).

best_values, best_lambda = stats.boxcox(income)
VERSIONS = [
    ("as measured, λ = 1", income, style.GREY),
    ("square root, λ = 0.5", np.sqrt(income), style.GOLD),
    ("log, λ = 0", np.log(income), style.GREEN),
    (f"Box-Cox, λ = {best_lambda:.2f}", best_values, style.BLUE),
]

fig, axes = style.subplots(4, cols=4, height=2.6, width=11)
for ax, (title, values, color) in zip(axes, VERSIONS):
    style.hist(ax, values, bins=30, color=color)
    skew = stats.skew(values)
    style.note(ax, f"skewness {skew:.2f}", x=0.5, color=style.INK)
    ax.set(title=title, yticks=[])
axes[0].set_ylabel("households")
fig.tight_layout()
plt.show()

Figure 3.8.5: The 235 household incomes as they are (λ = 1), square-rooted (λ = 0.5), logged (λ = 0), and with the Box-Cox λ chosen by maximum likelihood, here about −0.4. The skew shrinks step by step; the right two are close to a bell.

The skewness in each panel measures lopsidedness: 0 for a symmetric shape, positive for a long right tail. Transform for a reason (multiplying effects, a variance that grows with the mean), and report the results back on the original scale where you can: a difference of logs is a ratio.

3.8.6 What to do when an assumption fails

Problem	What it does	What to do
dependence (clusters, sequences)	standard errors far too small	model the structure: cluster summaries, mixed models, time series
unequal spreads	Student’s t and classical ANOVA unreliable	Welch’s t-test, Welch’s ANOVA, robust standard errors
strong skew, small samples	t-test intervals and p-values off	transform (log), rank tests, bootstrap
outliers	means and standard deviations dragged	check them; median, trimmed mean, rank tests; report with and without
counts, shares, yes/no	normal methods give impossible values	binomial, Poisson and logistic models (Chapters 4.4 and 4.6)

A trimmed mean drops a fixed share of the smallest and largest values (say 10% at each end) before averaging: a compromise between the mean and the median, and the basis of robust versions of the t-test.

3.8.7 In one picture

• Assumptions in order of importance: independence, equal spreads, normality.
• Check normality with QQ plots; Shapiro-Wilk misses problems in small samples and flags trivial ones in large samples.
• Check spreads with plots and Levene’s test, but simply use Welch’s versions by default.
• Clusters and sequences shrink the real amount of information; account for them.
• Transformations (log, square root, Box-Cox) tame skew; ranks and the bootstrap avoid the assumption altogether.

Next: comparing groups on several outcomes at once.