A test can be wrong in two ways: it can see an effect that is not there, or miss one that is. Small studies miss often and exaggerate when they hit; many tests find something by chance.
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.axes import Axes
from scipy import stats
from statsmodels.stats.power import TTestIndPower
from statly import style
style.use()
rng = np.random.default_rng(1)Every test ends in a decision, significant or not, and the truth is either “no effect” or “an effect”. That makes four outcomes, two of them wrong:
| no effect in truth | an effect in truth | |
|---|---|---|
| significant | false alarm (type I error) | a hit |
| not significant | correct | a miss (type II error) |
The chance of a false alarm is the significance level, usually 5%: we choose it. The chance of a hit, when there is an effect, is the power. We do not choose it directly. It depends on how large the effect is, how noisy the data are, and how many units we measure.
SHIFT, LINE = 2.2, 1.96 # the effect, in standard errors; the line for 5%, two-sided
x = np.linspace(-4, 6.5, 500)
def plot_errors(ax: Axes) -> None:
"""The two curves, the line, and the areas of false alarms, misses and power."""
null, effect = stats.norm.pdf(x), stats.norm.pdf(x, SHIFT)
for tail in (x >= LINE, x <= -LINE):
ax.fill_between(x[tail], null[tail], color=style.CORAL, alpha=0.5, lw=0)
missed = x < LINE
ax.fill_between(x[missed], effect[missed], color=style.GREY, alpha=0.35, lw=0)
ax.fill_between(x[~missed], effect[~missed], color=style.BLUE, alpha=0.2, lw=0)
ax.plot(x, null, color=style.GREY, ls="--")
ax.plot(x, effect, color=style.BLUE)
ax.axvline(LINE, color=style.INK, lw=1.4)
power = stats.norm.sf(LINE - SHIFT) + stats.norm.cdf(-LINE - SHIFT)
ax.text(LINE + 0.1, 0.42, "significant →", fontsize=9.5)
ax.text(SHIFT + 0.6, 0.2, f"power\n{power:.0%}", fontsize=10, color=style.BLUE)
ax.text(
LINE - 0.5, 0.05, f"misses\n{1 - power:.0%}", fontsize=10, ha="right", color=style.INK
)
ax.text(3.1, 0.015, "false alarms", fontsize=9, color=style.CORAL)
ax.set(xlabel="test statistic (standard errors from zero)", ylabel="how common", yticks=[])
fig, ax = plt.subplots(figsize=(8, 3.3))
plot_errors(ax)
plt.show()
The two errors pull against each other. Move the line to the right (a stricter level, like 1%) and false alarms fall, but misses rise. The only way to cut both is to separate the curves: measure more units, or measure them with less noise.
More units make the standard error smaller, so the same effect sits more standard errors from zero, and the blue curve moves right. Try it:
Before a study, turn it around: how many units do we need to have a good chance of finding the effect we care about? The usual target is 80% power.
sizes = np.arange(5, 401)
solver = TTestIndPower()
def plot_power(ax: Axes) -> None:
"""Power against group size for three effect sizes, with the size needed for 80%."""
label_offsets = {0.2: (-92, 8), 0.5: (8, -34), 0.8: (8, -16)} # keep the labels apart
for d, color in ((0.2, style.CORAL), (0.5, style.GOLD), (0.8, style.BLUE)):
power = solver.power(effect_size=d, nobs1=sizes, alpha=0.05)
ax.plot(sizes, power, color=color, label=f"d = {d}")
needed = int(np.ceil(solver.solve_power(effect_size=d, power=0.8, alpha=0.05)))
ax.plot(needed, 0.8, "o", ms=6, color=color)
ax.annotate(
f"{needed} per group",
(needed, 0.8),
xytext=label_offsets[d],
textcoords="offset points",
fontsize=9,
color=color,
)
ax.axhline(0.8, color=style.INK, ls="--", lw=1)
ax.yaxis.set_major_formatter(style.PERCENT)
ax.legend(fontsize=8, loc="lower right")
ax.set(xlabel="units per group", ylabel="power", ylim=(0, 1.02))
fig, ax = plt.subplots(figsize=(8, 3.4))
plot_power(ax)
plt.show()
A large effect needs 26 per group, a medium one 64, a small one almost 400. For a two-group comparison with n per group and effect size d (Chapter 3.4), the test statistic sits on average standard errors from zero, so with the normal approximation
where Φ is the cumulative normal curve, for a two-sided 5% test and for 80% power. For d = 0.5: n ≈ 2 × 2.8² / 0.25 ≈ 63 per group. Because d is squared in the denominator, halving the effect means four times the data. Many published studies have far less power than 80%, which matters for the next point.
A study with low power finds the effect only when chance happens to make it look large. So the significant results of small studies are, on average, too big. This is the winner’s curse.
TRUE_EFFECT, GROUP = 0.3, 20
a = rng.normal(0, 1, (10_000, GROUP))
b = rng.normal(TRUE_EFFECT, 1, (10_000, GROUP))
estimates = b.mean(axis=1) - a.mean(axis=1)
significant = stats.ttest_ind(b, a, axis=1).pvalue < 0.05
def plot_curse(ax: Axes) -> None:
"""All estimates, with the significant ones in red and their average."""
edges = np.linspace(-1, 1.6, 53)
ax.hist(estimates[~significant], bins=edges, color=style.BLUE, alpha=0.5)
ax.hist(estimates[significant], bins=edges, color=style.CORAL, alpha=0.7)
ax.axvline(TRUE_EFFECT, color=style.INK, ls="--", lw=1.4)
average = estimates[significant].mean()
ax.axvline(average, color=style.CORAL, lw=1.6)
style.note(
ax, f"{significant.mean():.0%} significant;\ntheir average: {average:.2f}", x=0.7
)
ax.set(xlabel="estimated effect (standard deviations)", ylabel="studies")
fig, ax = plt.subplots(figsize=(8, 3.2))
plot_curse(ax)
plt.show()
If journals publish mainly significant results, the published effects are inflated, and later, larger studies find smaller effects. Part of what is called the replication crisis is just this. The cure is more power, and reporting every result, not only the significant ones.
At the 5% level, one test in twenty gives a false alarm when nothing is going on. Run twenty tests and the chance of at least one false alarm is 1 − 0.95²⁰ = 64%. Test enough things, and something will always be significant.
TESTS, REAL = 100, 10
real = np.arange(TESTS) < REAL
shifts = np.where(real, 3.0, 0.0) # real effects sit 3 standard errors from zero
z = rng.normal(shifts)
p = 2 * stats.norm.sf(abs(z))
def plot_raw(ax: Axes) -> None:
"""Each test's p-value on a log scale, with the 5% line."""
order = rng.permutation(TESTS)
colors = np.where(real[order], style.BLUE, style.GREY)
ax.scatter(np.arange(TESTS), p[order], s=16, c=colors)
ax.axhline(0.05, color=style.INK, ls="--", lw=1)
false_alarms = int((p[~real] < 0.05).sum())
style.note(ax, f"false alarms at 5%: {false_alarms}", x=0.03, y=0.08, va="bottom")
ax.set_yscale("log")
ax.set(title="Uncorrected", xlabel="test", ylabel="p-value (log scale)")
def plot_corrected(ax: Axes, alpha: float = 0.05) -> None:
"""Sorted p-values against the Bonferroni line and the Benjamini-Hochberg line."""
order = np.argsort(p)
rank = np.arange(1, TESTS + 1)
colors = np.where(real[order], style.BLUE, style.GREY)
ax.scatter(rank, p[order], s=16, c=colors)
ax.axhline(alpha / TESTS, color=style.CORAL, lw=1.4, label="Bonferroni")
ax.plot(rank, alpha / (TESTS - rank + 1), color=style.GOLD, lw=1.4, ls="--", label="Holm")
ax.plot(rank, alpha * rank / TESTS, color=style.GREEN, lw=1.4, label="Benjamini-Hochberg")
ax.set_yscale("log")
ax.set_xlim(0, 40)
ax.legend(fontsize=8, loc="lower right")
ax.set(title="Corrected", xlabel="rank of the p-value", ylabel="p-value (log scale)")
fig, axes = style.subplots(2, height=3.3, width=9.5)
plot_raw(axes[0])
plot_corrected(axes[1])
fig.tight_layout()
plt.show()
There are two things to control. The family-wise error rate (FWER) is the chance of even one false alarm among all the tests. The false discovery rate (FDR) is the expected share of false alarms among the tests called significant. With m tests and V false alarms out of R discoveries:
Three ways to correct:
Use Holm (or Bonferroni) when a single false alarm is costly, Benjamini-Hochberg when screening many candidates (genes, products, features) for the ones worth a closer look. For all pairs of group means there is Tukey’s method (Chapter 3.5); for any comparison dreamed up after seeing the data, Scheffé’s even wider intervals. When tests are strongly correlated, permutation methods that shuffle all tests at once (keeping their correlation) are less wasteful than Bonferroni.
Next: meta-analysis, combining the results of many studies.
