How to adjust p-values with the Benjamini-Hochberg method in Statsmodels (with Code)

Editorial · Published 30 October 2026

The Short Answer (TL;DR)

Pass the complete planned family of p-values to multipletests with method="fdr_bh" to obtain Benjamini-Hochberg adjusted p-values and rejection decisions. Keep is_sorted=False for unsorted input so the results retain each hypothesis in its original position.

from statsmodels.stats.multitest import multipletests
p = [0.001,0.04,0.009,0.2,0.008,0.5]
reject, adjusted, _, _ = multipletests(p,alpha=.05,method="fdr_bh",is_sorted=False); print(adjusted,reject)

Detailed Implementation

This example uses synthetic data and does not need a network connection. Install the recorded package versions in a separate Python environment before running it.

python -m pip install numpy==2.5.4 scipy==1.18.1 pandas==3.0.6 statsmodels==0.15.0

from pathlib import Path
import numpy as np
from statsmodels.stats.multitest import multipletests

# Six illustrative p-values, already computed for one prespecified family.
labels = np.array(["H1", "H2", "H3", "H4", "H5", "H6"])
p = np.array([.001, .04, .009, .2, .008, .5])
alpha = .05
if p.ndim != 1 or p.size == 0 or not np.isfinite(p).all():
    raise ValueError("Use a nonempty one-dimensional vector of finite p-values.")
if np.any((p < 0) | (p > 1)) or labels.size != p.size:
    raise ValueError("Each p-value must lie in [0,1] and match one hypothesis label.")
reject, adjusted, alpha_sidak, alpha_bonferroni = multipletests(
    p, alpha=alpha, method="fdr_bh", is_sorted=False,
)
# The last two return values are Sidak and Bonferroni thresholds, not BH cutoffs.
lines = [f"family size={p.size}; BH FDR target={alpha:.2f}",
         "hypothesis raw_p adjusted_p reject"]
for label, raw, adj, decision in zip(labels,p,adjusted,reject):
    lines.append(f"{label} {raw:.3f} {adj:.3f} {bool(decision)}")
report = "\n".join(lines) + "\n"
Path("benjamini-hochberg.txt").write_text(report, encoding="utf-8")
print(report)

The output rejects H1, H3 and H5; all six raw and adjusted p-values remain in the original hypothesis order.

Which hypotheses are rejected in this example?

The adjusted p-values in input order are 0.006, 0.060, 0.018, 0.240, 0.018 and 0.500. Comparing them with the prespecified 0.05 target rejects H1, H3 and H5. The labels remain attached to their original p-values even though the calculation sorts internally; do not sort only the numbers and then attach the original labels.

Why is multiplying each sorted p-value by m/rank insufficient?

For the sorted values, the uncorrected rank-scaled numbers are 0.006, 0.024, 0.018, 0.060, 0.240 and 0.500. A reverse cumulative minimum changes the second entry to 0.018 so adjusted values cannot decrease as raw values increase. Clip at one and restore the original order. The checks reconstruct this formula independently and cover tied values, zero, one and reordered input.

What does the FDR target mean?

Let R be the number of rejections and V the number of rejected true null hypotheses. The false discovery rate is E[V/max(R,1)] across repeated experiments, with zero false-discovery proportion when there are no rejections. Standard BH control requires valid null p-values and independence or an appropriate positive-dependence condition, commonly expressed as PRDS on the true nulls. It is not assured merely because all sample correlations happen to be positive. This is not control of the probability of any false rejection, nor the probability that a particular rejected hypothesis is false.

How should the comparison family be defined?

Choose the set of hypotheses and their testing procedures before inspecting which p-values are favorable. Include the complete intended family, rather than correcting only results that passed a raw threshold. The six numbers here illustrate the adjustment and do not validate the underlying tests. For arbitrary dependence, the more conservative Benjamini-Yekutieli method is available as fdr_by; choose the procedure for the study design and inferential goal. The third and fourth return values from multipletests concern other methods and must not be reported as BH critical values.

Mathematical Context

Order m raw p-values as p_(1) through p_(m). The BH adjusted value at rank i is the smaller of one and the minimum of m p_(j)/j over all j from i through m. Map those adjusted values back to their original hypotheses. Equivalently, the step-up rule finds the largest k satisfying p_(k) <= k alpha/m and rejects ranks through k; if none qualifies, reject none. Here m=6 and k=3 because 0.009 <= 0.025 while the later ranks do not qualify. P-values, adjusted p-values and the FDR target are dimensionless; an adjusted p-value is not an effect size or a posterior probability.

\[ \widetilde{p}_{(i)}=\min\left(1,\min_{j\geq i}\frac{m p_{(j)}}{j}\right) \]

Reproducibility

Executed with Python 3.12.14; numpy 2.5.4, scipy 1.18.1, pandas 3.0.6, statsmodels 0.15.0. Automated checks confirm the stated output properties.