Randomized algorithms are convenient to test statistically: run the algorithm many times, compute a p-value, and assert that the null hypothesis is not rejected. The problem arises when there are several such tests in the same suite.
The multiple comparisons problem
Each individual test has a false rejection rate of α under the null hypothesis. Across m independent tests, the probability that at least one spuriously rejects is
\[\text{FWER} = 1 - (1 - \alpha)^m\]
With \(\alpha = 0.05\) and \(m = 10\) tests this is about 40%. A CI system running such a suite would see a spurious rejection — and therefore a spurious test failure — on almost every other build.
The fix is to adjust the per-test threshold so that the family-wise error rate stays at \(\alpha\). Plain Bonferroni does this by testing each hypothesis at \(\alpha/m\): simple, but wasteful. The interesting corrections are all step-down procedures, and they are easier to understand as variations on a single scheme than as separate recipes.
Step-down procedures
Each raw p-value \(p_k\) is mapped to an adjusted value \(\tilde{p}_k\), and null hypothesis \(k\) is rejected when \(\tilde{p}_k \le \alpha\). The p-values are adjusted to be monotonically increasing, so the first time we fail to reject a null hypothesis, we’ll fail to reject the rest of the null hypotheses as well. To be specific, say we have \(K\) tests with p-values sorted ascending as \(p_1 \le p_2 \le \cdots \le p_K\). We set \(\tilde{p}_k = \max\left(\tilde{p}_{k-1},P( c(S_k) < p_k) \right)\) where \(S_k = \{k, k+1, \dotsc, K\}\), the random variable \(c(S) = \min_{j \in S} P_j\) , and \(P\) is distributed according to the joint null distribution.
Let \(I_0\) be the set of true null hypotheses. Say \(L\) is the first true null hypothesis in our ordering. That is, \(p_L\) is an observation of \(c(I_0)\). By monotonicity, the event that we reject a true null hypothesis is the same as the event that we reject hypothesis \(L\). We reject at \(L\) when \(P( c(S_L) \leq p_L) \leq \alpha\). Because \(L\) is the first true null, \(I_0 \subseteq S_L\). So \(c(S_L) \leq c(I_0)\) and \(P(c(S_L) \leq p_L ) \geq P(c(I_0) \leq p_L)\). Let \(F\) be the CDF of \(c(I_0)\). This means that \(F(p_L) \leq \alpha\) or \(p_L \leq F^{-1}(\alpha)\) when we reject. As \(p_L\) has CDF \(F\), this occurs with probability \(\alpha\).
Where the procedures differ is how they estimate \(P(c(S_k) \leq p_k)\). For Westfall-Young, we can simulate samples of \(c(S_k)\) and count the number of samples below \(p_k\). This means the number of samples must exceed \(1/\alpha\) for any rejection to be possible.
For Holm Bonferroni, we use a union bound: the probability the smallest \(p\) value is below a threshold is at most the probability that any of the p-values are below the threshold: \(P(c(S_k) \leq p_k) \leq |S_k|p_k\).
For the Sidak correction, we assume independence: \(P(c(S_k) \leq p_k) = 1 - (1 - p_k)^{|S_k|}\).
So the choice is really a choice about how much I claim to know about the joint null. A union bound assumes nothing and pays for it in power. Šidák assumes independence. Westfall-Young assumes I can simulate the joint null, which is true whenever I can say concretely what “nothing is going on” would look like for my data.
Why this is awkward in pytest
Any of these is inherently a post-hoc procedure: \(\tilde{p}_k\) depends on every p-value ranked at or below \(k\), so rank \(k\) cannot be evaluated until all \(K\) p-values are in hand. pytest, meanwhile, determines pass/fail for each test as it runs, long before the rest of the suite has completed.
What the plugin does
pytest-familywise adds an assertNotReject fixture. Tests call it with a sampler: a callable handed a random generator, which draws a group of data and returns the p-value for it. The plugin owns the generator, which keeps draws reproducible and lets the same machinery draw from the null distribution when the correction needs it. The fixture stores the p-value but defers pass/fail entirely. After the session ends, the plugin:
- Sorts all registered p-values ascending.
- Applies the step-down procedure, marking each test passed or failed.
- Retroactively moves reports between pytest’s
passedandfailedstat buckets so the final summary line reflects the corrected outcomes. - Sets the session exit code to non-zero if any test failed.
The implementation hooks into pytest_runtest_logreport to capture reports, pytest_sessionfinish to run the correction and update session.exitstatus, and pytest_terminal_summary to rewrite the stat buckets and print the correction table before pytest prints its own N passed, M failed line.
Holm-Bonferroni: independent draws
Suppose I am testing a custom random number generator and want to verify three properties jointly, with FWER controlled at 5%. Each test draws its own fresh data, so there is nothing to exploit about their joint distribution and the union bound is the honest choice.
# tests/test_rng.py
import scipy.stats
def test_uniform_marginals(ks_sample_size, assertNotReject):
n = ks_sample_size(effect_size=0.05) # one-sample KS, ||F−G||∞ ≥ 0.05
assertNotReject(
lambda rng: scipy.stats.kstest(rng.random(n), "uniform").pvalue)
def test_normal_mean_zero(ztest_sample_size, assertNotReject):
n = ztest_sample_size(effect_size=0.3) # Cohen's d = 0.3
assertNotReject(
lambda rng: scipy.stats.ttest_1samp(rng.standard_normal(n), 0.0).pvalue)
def test_discrete_distribution(chisquare_sample_size, assertNotReject):
n = chisquare_sample_size(effect_size=0.2, df=4) # Cohen's w = 0.2
assertNotReject(
lambda rng: scipy.stats.chisquare(rng.multinomial(n, [0.2] * 5)).pvalue)pytest --correction=holm --alpha=0.05 --power=0.8Output after all three tests complete:
============ Holm-Bonferroni correction α=0.05 n=3 =============
PASSED p=0.312541 threshold=0.016667 test_rng.py::test_uniform_marginals
PASSED p=0.487302 threshold=0.025000 test_rng.py::test_normal_mean_zero
PASSED p=0.621088 threshold=0.050000 test_rng.py::test_discrete_distribution
3 passed, 0 failed after Holm-Bonferroni correctionAll three p-values are large (data consistent with H0), so none are rejected and all tests pass. The thresholds tighten for the lowest-ranked p-values (0.017 for rank 1, relaxing to 0.050 for rank 3). If the rank-1 test had returned p = 0.01, it would be rejected (0.01 < 0.017) and that test would fail — but the remaining tests would still pass, since the step-down procedure stops rejecting once it fails to reject.
Westfall-Young: one dataset, several angles
The union bound is needlessly pessimistic when several tests chew on the same dataset from different angles. Their p-values move together, and treating them as three independent chances to be wrong throws away power we already paid for. Simulating the joint null recovers it. Each test supplies a null_sample callable alongside its statistic: something that recomputes the same p-value with \(H_0\) forced onto the observed data.
import numpy as np
import pytest
import scipy.stats
@pytest.fixture(scope="module")
def samples():
"""One dataset, several tests — so their p-values are correlated."""
return np.random.default_rng(0).standard_normal(2000)
def under_h0(rng, data):
"""Impose H0 (mean zero) on the observed data, then bootstrap."""
centered = data - data.mean()
return rng.choice(centered, size=len(data), replace=True)
def test_mean_zero(assertNotReject, samples):
assertNotReject(
lambda rng: scipy.stats.ttest_1samp(samples, 0.0).pvalue,
null_sample=lambda rng:
scipy.stats.ttest_1samp(under_h0(rng, samples), 0.0).pvalue)
def test_median_zero(assertNotReject, samples):
"""A different statistic on the same data — strongly correlated."""
def pvalue(data):
return scipy.stats.binomtest(int((data > 0).sum()), len(data), 0.5).pvalue
assertNotReject(
lambda rng: pvalue(samples),
null_sample=lambda rng: pvalue(under_h0(rng, samples)))pytest --correction=westfall-young --resamples=1000The resample count is not a free knob: as above, it has to exceed \(1/\alpha\) before any rejection is even representable. A thousand at \(\alpha = 0.05\) is comfortable.
Bootstrapping a centered copy is only one way to impose a null. When \(H_0\) says two groups are exchangeable, the transform is a shuffle of the pooled data:
def test_groups_are_identical(assertNotReject, group_a, group_b):
pooled = np.concatenate([group_a, group_b])
n = len(group_a)
def pvalue(x):
return scipy.stats.ks_2samp(x[:n], x[n:]).pvalue
assertNotReject(
lambda rng: pvalue(pooled),
null_sample=lambda rng: pvalue(rng.permutation(pooled)))Sign-flipping (data * rng.choice([-1, 1], size=len(data))) is the third common idiom, for \(H_0\)s that assert symmetry about zero. Each is a different way of saying “here is what this dataset would have looked like if nothing were going on,” and the plugin does not care which one I pick.
Sample-size fixtures
Each test above starts by asking the plugin how many samples are needed to achieve the requested power. This is worth doing explicitly rather than picking a number by hand: a test with insufficient power will produce false negatives silently, and a test with excess samples wastes time.
The plugin provides three fixtures, one per test family. They all read --alpha and --power from the command line, so the targets are set once and the sizing follows automatically.
z-test (ztest_sample_size): closed-form solution
\[n = \left\lceil \left(\frac{z_\alpha + z_\beta}{d}\right)^2 \right\rceil\]
where \(d\) is Cohen’s \(d\) and \(z_\alpha\), \(z_\beta\) are the appropriate normal quantiles. At \(\alpha = 0.05\), power \(= 0.8\), \(d = 0.5\) this gives \(n = 32\).
chi-square (chisquare_sample_size): numerical root-finding on the non-central \(\chi^2\) survival function. The non-centrality parameter is \(\lambda = n w^2\) where \(w\) is Cohen’s \(w\); we find the smallest integer \(n\) such that \(P(\chi^2(df, \lambda) > c_{1-\alpha}) \ge \text{power}\).
Kolmogorov-Smirnov (ks_sample_size): closed form derived from the Dvoretzky-Kiefer-Wolfowitz inequality. Given the maximum CDF discrepancy \(\Delta\) and setting \(\beta = 1 - \text{power}\), the bound is:
\[n \ge \frac{\left(\sqrt{\ln(2/\alpha)} + \sqrt{\ln(2/\beta)}\right)^2}{2\Delta^2}\]
The derivation equates the Type-I and Type-II DKW boundaries at the effect size, then solves for \(n\). For a two-sample test with equal group sizes the effective \(n\) is \(n_\text{each}/2\), so the returned per-group count is doubled.
All three fixtures treat --power as a per-test rate, not a family-wise one. The sizing itself uses the conservative bound \(\alpha/m\), where \(m\) is the number of collected tests that use assertNotReject; a test that uses a sample-size fixture without assertNotReject is not in the family, sizes against the raw alpha, and does not raise \(m\) for anyone else.
Loading the plugin
The package registers itself via a pytest11 entry point, so installing it is sufficient:
uv add --dev pytest-familywise