An estimate without its uncertainty is half an answer. A confidence interval gives the other half: a range of values the data agree with.
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)
visits = (
sm.datasets.randhie.load_pandas().data.mdvis.to_numpy()
) # doctor visits per person per year (RAND study)
crime = sm.datasets.statecrime.load_pandas().data # crime and poverty in 50 US states and DCTake the 20,190 people of the RAND health study as a population. On average they saw a doctor 2.9 times a year. Now pretend we could only afford to ask some of them. A sample of 25 gives an average, another sample of 25 gives another. How much do those averages vary?
TRUE_MEAN, SD = visits.mean(), visits.std()
def plot_sample_means(ax: Axes, n: int) -> None:
"""Averages of many samples of size n, with ± one standard error marked."""
means = rng.choice(visits, (5_000, n)).mean(axis=1)
step = (
np.ceil(0.1 * n) / n
) # averages of whole numbers land on steps of 1/n: whole steps per bin
style.hist(ax, means, bins=np.arange(-0.5 / n, 7, step))
ax.axvline(TRUE_MEAN, color=style.INK, ls="--", lw=1.2)
standard_error = SD / np.sqrt(n)
top = ax.get_ylim()[1] * 0.9
ax.plot(
[TRUE_MEAN - standard_error, TRUE_MEAN + standard_error],
[top, top],
color=style.CORAL,
lw=2,
marker="|",
ms=10,
)
style.note(ax, f"standard error\n{standard_error:.2f}", x=0.62)
ax.set(title=f"Samples of {n}", xlabel="average visits per year", ylabel="samples")
fig, axes = style.subplots(3, height=3, width=10.5, sharex=True)
for ax, n in zip(axes, (25, 100, 400)):
plot_sample_means(ax, n)
fig.tight_layout()
plt.show()
The spread of an estimate from sample to sample is its standard error. For an average it is
The standard deviation describes how much people differ; the standard error describes how much averages differ. A sample of 400 pins the average down twice as well as a sample of 100, not four times. In practice we have one sample, so we estimate the standard deviation from it.
By the central limit theorem, the averages are close to normal, so 95% of them land within about two standard errors of the true value. Turn that around: the range average ± 2 standard errors contains the true value for 95% of samples. That range is a 95% confidence interval.
SAMPLES, N = 50, 100
def interval(sample: np.ndarray) -> tuple[float, float, float]:
"""The average and its 95% confidence interval (with the t distribution)."""
mean = sample.mean()
half = stats.t.ppf(0.975, len(sample) - 1) * sample.std(ddof=1) / np.sqrt(len(sample))
return mean, mean - half, mean + half
def plot_coverage(ax: Axes) -> None:
"""One horizontal interval per sample, red if it misses the true average."""
hits = 0
for i in range(SAMPLES):
mean, low, high = interval(rng.choice(visits, N))
hit = low <= TRUE_MEAN <= high
hits += hit
color = style.BLUE if hit else style.CORAL
ax.plot([low, high], [i, i], color=color, lw=1.6)
ax.plot(mean, i, "o", ms=3, color=color)
ax.axvline(TRUE_MEAN, color=style.INK, ls="--", lw=1.2)
style.note(ax, f"{hits} of {SAMPLES} intervals\ncontain the true average", x=0.72, y=0.14)
ax.set(xlabel="average visits per year", ylabel="sample", yticks=[])
fig, ax = plt.subplots(figsize=(8, 4.6))
plot_coverage(ax)
plt.show()
The 95% belongs to the method, not to one interval. Any single interval either contains the true value or it does not; we just do not know which. What we know is that the recipe catches the truth 19 times out of 20.
The exact multiplier is not 2 but comes from the t distribution of Chapter 2.1, because the standard deviation is estimated: 2.26 for a sample of 10, 2.04 for 30, 1.98 for 100. For large samples it is 1.96. In full, the interval for a mean is
where is the point of the t distribution with n − 1 degrees of freedom that has 97.5% below it (2.5% in each tail leaves 95% in the middle). Two things set the width: the confidence level and the sample size.
SIZES = np.arange(3, 201)
def plot_widths(ax: Axes) -> None:
"""t multiplier times 1/√n for three confidence levels."""
for level, color in ((0.90, style.GREEN), (0.95, style.BLUE), (0.99, style.CORAL)):
half = stats.t.ppf(0.5 + level / 2, SIZES - 1) / np.sqrt(SIZES)
ax.plot(SIZES, half, color=color, label=f"{level:.0%} confidence")
ax.set_ylim(0, 1.6)
ax.legend(fontsize=8)
ax.set(xlabel="sample size n", ylabel="half-width (standard deviations)")
fig, ax = plt.subplots(figsize=(8, 3.2))
plot_widths(ax)
plt.show()
The same recipe works for any estimate whose standard error we know: estimate ± about 2 standard errors.
That second interval is wide. The experiment says the method helps, but it cannot say whether it helps a little or a lot. That is what 32 students buy.
The simple interval for a share, (named after Wald), has a flaw: near 0% or 100%, and with small samples, it catches the truth far less than 95% of the time. With 0 successes it even has width zero. Wilson’s interval fixes this by working out which values of p would make the observed share unsurprising:
TRIES = 20
truths = np.linspace(0.01, 0.99, 197)
def coverage(method: str) -> np.ndarray:
"""Exact coverage: the probability, summed over all outcomes, that the interval holds the truth."""
k = np.arange(TRIES + 1)
p_hat = k / TRIES
if method == "wald":
half = 1.96 * np.sqrt(p_hat * (1 - p_hat) / TRIES)
low, high = p_hat - half, p_hat + half
else:
z2 = 1.96**2
centre = (p_hat + z2 / (2 * TRIES)) / (1 + z2 / TRIES)
half = (
1.96 * np.sqrt(p_hat * (1 - p_hat) / TRIES + z2 / (4 * TRIES**2)) / (1 + z2 / TRIES)
)
low, high = centre - half, centre + half
return np.array(
[stats.binom.pmf(k, TRIES, p)[(low <= p) & (p <= high)].sum() for p in truths]
)
def plot_coverage(ax: Axes) -> None:
"""Coverage against the true chance, for both intervals."""
ax.plot(truths, coverage("wald"), color=style.CORAL, lw=1.2, label="Wald (simple)")
ax.plot(truths, coverage("wilson"), color=style.BLUE, lw=1.2, label="Wilson")
ax.axhline(0.95, color=style.INK, ls="--", lw=1)
ax.xaxis.set_major_formatter(style.PERCENT)
ax.yaxis.set_major_formatter(style.PERCENT)
ax.legend(fontsize=8, loc="lower center")
ax.set(xlabel="true chance", ylabel="intervals that contain it", ylim=(0.5, 1))
fig, ax = plt.subplots(figsize=(8, 3.2))
plot_coverage(ax)
plt.show()
Use Wilson’s interval (or an exact one) for shares, especially with few cases or shares near 0 or 1. The saw-tooth comes from the counts being whole numbers.
For many estimates there is no formula for the standard error, or the formula assumes a normal shape the data do not have. The bootstrap from Chapter 1.1 fills the gap: resample the data with replacement, compute the estimate each time, and take the middle 95% of the results. Chapter 3.3 goes through the bootstrap in detail.
x, y = crime.poverty.to_numpy(), crime.murder.to_numpy()
r = np.corrcoef(x, y)[0, 1]
picks = rng.integers(0, len(x), (10_000, len(x)))
resampled = np.array([np.corrcoef(x[i], y[i])[0, 1] for i in picks])
low, high = np.quantile(resampled, [0.025, 0.975])
def plot_states(ax: Axes) -> None:
"""Poverty against murder rate, one point per state."""
ax.scatter(x, y, s=16, color=style.BLUE, alpha=0.7)
ax.annotate(
"DC",
(x[crime.index == "District of Columbia"][0], y.max()),
xytext=(6, -3),
textcoords="offset points",
fontsize=9,
)
style.note(ax, f"correlation {r:.2f}", x=0.03)
ax.set(title="States", xlabel="poverty (% of people)", ylabel="murders per 100,000")
def plot_resampled(ax: Axes) -> None:
"""Bootstrap correlations, with the middle 95% marked."""
style.hist(ax, resampled, bins=50)
for edge in (low, high):
ax.axvline(edge, color=style.CORAL, lw=1.6)
style.note(ax, f"95% interval\n{low:.2f} to {high:.2f}", x=0.03)
ax.set(title="Bootstrap", xlabel="correlation in a resample", ylabel="resamples")
fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_states(axes[0])
plot_resampled(axes[1])
fig.tight_layout()
plt.show()
The bootstrap distribution is lopsided, and a formula assuming a bell would miss that. It also shows something useful: how much one unusual point matters. Resamples without Washington DC give clearly higher correlations (about 0.61 on average, against 0.55 with it): DC’s murder rate is so extreme that it sits far above the line and weakens the straight-line fit.
A 95% interval and a test at the 5% level always agree: the interval leaves out 0 exactly when the test is significant. But the interval says more. It shows how large the effect could be, which a p-value hides.
STUDIES = [ # (name, estimate, low, high): change in blood pressure, mmHg
("A: large and clear", 12, 7, 17),
("B: clear but tiny", 0.8, 0.3, 1.3),
("C: unclear, could be large", 6, -2, 14),
("D: clearly about zero", 0.2, -0.6, 1.0),
]
MATTERS = 5 # mmHg: the smallest effect a doctor would care about
def plot_studies(ax: Axes) -> None:
"""Estimates with intervals against zero and against what matters in practice."""
ax.axvspan(MATTERS, 20, color=style.GREEN, alpha=0.08, lw=0)
ax.axvline(0, color=style.INK, lw=1)
ax.text(MATTERS + 0.3, -0.75, "large enough to matter", fontsize=9, color=style.GREEN)
for i, (name, estimate, low, high) in enumerate(STUDIES):
significant = low > 0 or high < 0
color = style.BLUE if significant else style.GREY
ax.plot([low, high], [i, i], color=color, lw=2.5)
ax.plot(estimate, i, "o", ms=7, color=color)
ax.set_yticks(range(len(STUDIES)), [s[0] for s in STUDIES])
ax.invert_yaxis()
ax.set(
xlabel="drop in blood pressure (mmHg)", xlim=(-4, 20), ylim=(len(STUDIES) - 0.4, -1.1)
)
fig, ax = plt.subplots(figsize=(8, 3.2))
plot_studies(ax)
plt.show()
“Significant” means only that the interval leaves out zero. Whether an effect matters is a separate question, about where the interval lies. Report both: the estimate with its interval first, the p-value second.
Next: the standard tests for comparing two groups.