3.1 The logic of a test

A test asks one question: if nothing were going on, how often would chance alone produce a pattern at least this strong? If the answer is “rarely”, we doubt that nothing is going on.

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

from statly import style

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

spector = sm.datasets.spector.load_pandas().data  # 32 students, a new teaching method or not
new_method = spector.PSI.to_numpy().astype(int)  # 1 = taught with the new method
improved = spector.GRADE.to_numpy().astype(int)  # 1 = the grade improved after the course
pre_test = spector.TUCE.to_numpy()  # score on a test taken before the course

3.1.1 A question and a null hypothesis

In a small experiment, 32 economics students took the same course. Fourteen were taught with a new method, eighteen the usual way. Afterwards, 8 of the 14 had improved their grade, against 3 of the 18. That is 57% against 17%. Does the method work, or were the 14 just lucky?

To answer, we imagine the dullest possible world: the method does nothing. This is the null hypothesis. In that world, each student would have improved or not whatever the method, and which group a student was in tells us nothing. The labels “new” and “usual” could be shuffled without changing anything that matters.

So we shuffle them, and see what differences chance alone produces:

def difference(values: np.ndarray, labels: np.ndarray) -> float:
    """Average with the new method minus average with the usual one (for 0/1: a share)."""
    return values[labels == 1].mean() - values[labels == 0].mean()


def plot_students(ax: Axes, labels: np.ndarray, title: str) -> None:
    """Two rows of students, improved ones filled, with each group's share."""
    for row, (group, name) in enumerate(((0, "usual"), (1, "new"))):
        outcome = np.sort(improved[labels == group])[::-1]  # improved students first
        x = np.arange(len(outcome))
        ax.scatter(x, np.full(len(x), row), s=90, facecolors="none", edgecolors=style.INK)
        ax.scatter(x[outcome == 1], np.full(outcome.sum(), row), s=90, color=style.CORAL)
        ax.text(19, row, f"{outcome.sum()} of {len(x)}", va="center", fontsize=10)
        ax.text(-1.2, row, name, va="center", ha="right", fontsize=10)
    ax.set(title=f"{title}: difference {difference(improved, labels):+.0%}", xlim=(-4, 22))
    ax.set(ylim=(-0.7, 1.7), xticks=[], yticks=[])
    ax.grid(False)
    for side in ("left", "bottom"):
        ax.spines[side].set_visible(False)


fig, axes = style.subplots(2, height=1.9, width=10)
plot_students(axes[0], new_method, "The experiment")
plot_students(axes[1], rng.permutation(new_method), "One shuffle")
fig.tight_layout()
plt.show()

Figure 3.1.1: Each circle is a student; filled circles improved. Left: the real experiment. Right: the same 32 students with the group labels shuffled at random; the difference shrinks to what chance gives.

One shuffle tells us little. Ten thousand shuffles show the whole range of what chance can do:

observed = difference(improved, new_method)
shuffles = np.array([difference(improved, rng.permutation(new_method)) for _ in range(10_000)])
as_large = shuffles >= observed - 1e-9
p_value = as_large.mean()


def plot_shuffles(ax: Axes) -> None:
    """How often each difference turned up, with the ones at least as large as the real one."""
    values, counts = np.unique(np.round(shuffles, 6), return_counts=True)
    colors = np.where(values >= observed - 1e-9, style.CORAL, style.BLUE)
    ax.bar(values, counts, width=0.07, color=colors, alpha=0.7)
    ax.axvline(observed, color=style.CORAL, lw=1.6)
    style.note(
        ax,
        f"{as_large.sum()} of 10,000 shuffles\nreach +{observed:.0%}:\np = {p_value:.3f}",
        x=0.03,
    )
    ax.xaxis.set_major_formatter(style.PERCENT)
    ax.set(xlabel="difference in share improved (new − usual)", ylabel="shuffles")


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

Figure 3.1.2: The difference in share improved for 10,000 random shuffles of the group labels: what chance alone produces. The real difference (+40 percentage points, red line) is out in the tail; shuffles at least that large (red bars) are rare.

3.1.2 The p-value

That share, the fraction of shuffles at least as extreme as the real result, is the p-value. Here it is about 2%: if the method did nothing, a difference this large would turn up in about 2 out of 100 experiments like this one. That is rare enough to doubt the null hypothesis.

Three details matter:

• “At least as extreme”, not “exactly this”. Any single difference is unlikely; the question is how far into the tail the result lies.
• One side or two. Here we asked whether the new method is better, so we counted one tail. If either direction would have been news, count both tails: a two-sided test, with roughly double the p-value. Decide which before looking at the data.
• The null is a model. “Shuffle the labels” is a fine model when students were put in groups at random. In an observational study, the groups may differ in other ways, and a small p-value only says the groups differ, not why.

Shuffling the labels is called a permutation test. It needs no distribution and no formula, only the data and the null hypothesis. Every other test in this part is a shortcut to the same answer.

3.1.3 What a p-value is not

The p-value is the most misunderstood number in science. It is not the probability that the null hypothesis is true, and it is not the size of the effect. It is a statement about the data under one assumption: if nothing were going on, data like these would be this rare.

A simulation makes it concrete. Run 10,000 small experiments, each comparing two groups of 20, once in a world where nothing is going on and once where there is a real effect:

EXPERIMENTS, GROUP_SIZE = 10_000, 20


def simulated_p_values(effect: float) -> np.ndarray:
    """p-values of t-tests on simulated experiments with a given true effect."""
    a = rng.normal(0, 1, (EXPERIMENTS, GROUP_SIZE))
    b = rng.normal(effect, 1, (EXPERIMENTS, GROUP_SIZE))
    return stats.ttest_ind(a, b, axis=1).pvalue


def plot_p_values(ax: Axes, effect: float, title: str) -> None:
    """A histogram of p-values in bins of 0.05, with the first bin in red."""
    p = simulated_p_values(effect)
    counts, edges = np.histogram(p, bins=np.linspace(0, 1, 21))
    colors = [style.CORAL] + [style.BLUE] * 19
    ax.bar(edges[:-1], counts, width=0.05, align="edge", color=colors, alpha=0.7)
    ax.set_ylim(0, counts.max() * 1.3)
    style.note(ax, f"below 0.05: {(p < 0.05).mean():.0%}", x=0.45)
    ax.set(title=title, xlabel="p-value", ylabel="experiments")


fig, axes = style.subplots(2, height=3, width=9)
plot_p_values(axes[0], 0.0, "No effect")
plot_p_values(axes[1], 0.5, "A real effect")
fig.tight_layout()
plt.show()

Figure 3.1.3: p-values from 10,000 simulated experiments with two groups of 20. Left: with no effect, p-values are spread evenly between 0 and 1, so 5% of them land below 0.05 by chance alone (red). Right: with a real effect of half a standard deviation, they pile up near 0, but many still land above 0.05.

With no effect, every p-value is equally likely: the distribution is uniform, as promised in Chapter 2.1. The 5% below 0.05 are false alarms. With a real effect, small p-values become common, but a third or more of these experiments still miss it: groups of 20 are small. A p-value above 0.05 does not show that there is no effect; it only fails to show that there is one.

The usual line is 0.05, the significance level. Results below it are called statistically significant. The line is a convention, not a law of nature: p = 0.049 and p = 0.051 are the same evidence. Report the p-value itself, and the size of the effect next to it.

3.1.4 Shortcuts: tests with a formula

Shuffling is easy for a computer, but for a century it was not. Classical tests use the distributions of Chapter 2.1 instead: the shuffles often pile up in a known shape, so the p-value can be read from a curve.

Before the course, the students took a test, so we can check whether the two groups started out alike. If the groups were formed at random, they should differ only by chance:

pre_difference = difference(pre_test, new_method)
pre_shuffles = np.array(
    [difference(pre_test, rng.permutation(new_method)) for _ in range(10_000)]
)


def plot_shortcut(ax: Axes) -> None:
    """Shuffled differences against the curve of the t-test, with both p-values."""
    style.hist(ax, pre_shuffles, bins=40, density=True)
    a, b = pre_test[new_method == 1], pre_test[new_method == 0]
    standard_error = np.sqrt(a.var(ddof=1) / len(a) + b.var(ddof=1) / len(b))
    x = np.linspace(-5, 5, 300)
    ax.plot(x, stats.norm.pdf(x, 0, standard_error), color=style.CORAL)
    ax.axvline(pre_difference, color=style.INK, ls="--", lw=1.5)
    shuffled_p = (abs(pre_shuffles) >= abs(pre_difference) - 1e-9).mean()
    t_test_p = stats.ttest_ind(a, b, equal_var=False).pvalue
    style.note(ax, f"shuffles: p = {shuffled_p:.2f}\nt-test: p = {t_test_p:.2f}", x=0.03)
    ax.set(xlabel="difference in average pre-test score (new − usual)", ylabel="how common")
    ax.set_yticks([])


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

Figure 3.1.4: The difference in average pre-test score between the groups, for 10,000 shuffles (bars), with the normal curve a t-test uses (red). The real difference (dashed) sits in the middle: the groups started out alike. The shuffles and the curve give almost the same p-value.

This time the real difference is ordinary, and both methods give p ≈ 0.55 (two-sided): no sign that the groups differed before the course. That supports reading the improvement as an effect of the method.

Every classical test follows the same recipe: a test statistic that measures the pattern (a difference, a ratio, a correlation), a null distribution that says what chance alone produces, and a p-value that says where the data fall in it. The next chapters fill in the recipe for common questions.

3.1.5 A test in symbols

Every test can be written the same way. The null hypothesis $H_0$ is the dull world, such as “the method makes no difference”. The alternative $H_1$ is what we suspect instead: “the new method is better” (one-sided) or “the methods differ” (two-sided). A test statistic $T$ measures the pattern, and the p-value is the probability, if $H_0$ were true, of a statistic at least as extreme as the one observed, $t_{\text{obs}}$:

$$p = P\big(T \ge t_{\text{obs}} \mid H_0\big) \quad \text{(one-sided)}, \qquad p = P\big(|T| \ge |t_{\text{obs}}| \mid H_0\big) \quad \text{(two-sided)}$$

Before looking at the data we fix the significance level α, usually 5%. Values of T so extreme that they occur with probability α under $H_0$ form the rejection region; if the observed statistic lands there (equivalently, if p < α), we reject $H_0$.

OBSERVED = 1.8
z = np.linspace(-4, 4, 400)


def plot_region(ax: Axes, two_sided: bool) -> None:
    """The null curve, the rejection region(s), and the observed value with its p-value."""
    curve = stats.norm.pdf(z)
    ax.plot(z, curve, color=style.BLUE)
    cut = stats.norm.ppf(0.975 if two_sided else 0.95)
    tails = [z >= cut] + ([z <= -cut] if two_sided else [])
    for tail in tails:
        ax.fill_between(z[tail], curve[tail], color=style.CORAL, alpha=0.45, lw=0)
    ax.axvline(OBSERVED, color=style.INK, ls="--", lw=1.4)
    p = 2 * stats.norm.sf(OBSERVED) if two_sided else stats.norm.sf(OBSERVED)
    style.note(
        ax,
        f"reject beyond {'±' if two_sided else ''}{cut:.2f}\np = {p:.3f}",
        x=0.03,
        color=style.INK,
    )
    ax.set(
        title="Two-sided" if two_sided else "One-sided", xlabel="test statistic z", yticks=[]
    )


fig, axes = style.subplots(2, height=3, width=9.5, sharey=True)
plot_region(axes[0], two_sided=False)
plot_region(axes[1], two_sided=True)
axes[0].set_ylabel("how common if H₀ is true")
fig.tight_layout()
plt.show()

Figure 3.1.5: A test statistic that follows a standard normal curve under the null hypothesis, observed at z = 1.8 (dashed). Left: one-sided, the 5% rejection region is everything above 1.64, and p = 3.6%: significant. Right: two-sided, 2.5% in each tail, beyond ±1.96, and p = 7.2%: not significant. The same data, a different question.

This is why the choice of one or two sides must be made before seeing the data: picking the side afterwards halves the p-value for free.

3.1.6 The first formula: a test for a share

In Chapter 1.1 a coin gave 15 heads in 20 tosses. Is it fair? Here $H_0$ is “the chance of heads is 50%”. Under $H_0$, the number of heads follows the binomial distribution of Chapter 2.2, so the p-value can be computed exactly: add up the binomial probabilities of all outcomes at least as far from 10 as 15 is. That is the binomial test.

For larger samples, the binomial is close to a normal curve (the central limit theorem), and the share of heads $\hat p$ can be turned into a z-score:

$$z = \frac{\hat p - p_0}{\sqrt{p_0(1 - p_0)/n}}$$

the distance from the null value $p_0$ in standard errors.

FLIPS, HEADS = 20, 15
k = np.arange(FLIPS + 1)
probability = stats.binom.pmf(k, FLIPS, 0.5)
extreme = abs(k - FLIPS / 2) >= abs(HEADS - FLIPS / 2)
exact = stats.binomtest(HEADS, FLIPS, 0.5).pvalue
z_value = (HEADS / FLIPS - 0.5) / np.sqrt(0.25 / FLIPS)


def plot_binomial_test(ax: Axes) -> None:
    """Binomial probabilities, extreme outcomes in red, with the normal curve."""
    ax.bar(k, probability, color=np.where(extreme, style.CORAL, style.BLUE), alpha=0.7)
    grid = np.linspace(0, FLIPS, 300)
    ax.plot(grid, stats.norm.pdf(grid, FLIPS / 2, np.sqrt(FLIPS / 4)), color=style.INK, lw=1.4)
    style.note(
        ax,
        f"exact binomial test: p = {exact:.3f}\nz = {z_value:.2f}, normal approximation: p = {2 * stats.norm.sf(z_value):.3f}",
        x=0.6,
        color=style.INK,
    )
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.set(xlabel="heads in 20 tosses", ylabel="probability if fair")


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

Figure 3.1.6: Heads in 20 tosses of a fair coin (bars). The red bars are the outcomes at least as far from 10 as the observed 15, in both directions; together they hold the exact two-sided p-value. The curve is the normal approximation used by the z-test.

The exact test gives p = 0.041, the normal shortcut 0.025: with only 20 tosses the shortcut is rough, and the exact test is better. A common rule: use the normal approximation only when both $np_0$ and $n(1 - p_0)$ are at least 10.

3.1.7 In one picture

• The null hypothesis is the world where nothing is going on.
• The p-value is how often that world produces a result at least as extreme as ours.
• A permutation test builds that world by shuffling labels; classical tests read it from a curve.
• With no effect, p-values are uniform: 5% fall below 0.05 by chance.
• A test has a null and an alternative hypothesis, a statistic, a level α and a rejection region; decide one or two sides before looking.
• The binomial test is exact for a share; the z-test is its large-sample shortcut.
• A p-value is not the chance that the null is true, and not the size of the effect. Large p means “no evidence”, not “no effect”.

Next: how sure we can be of a number, with confidence intervals.