Multiple testing correction in Python: Holm, BH and Romano–Wolf

A study that tests ten outcomes, five subgroups or several treatment arms runs many hypothesis tests at once. Each test at the 5% level has a 5% chance of a false rejection when its null is true, but the chance that at least one of many true nulls is rejected is much higher. Multiple testing corrections adjust the p-values so that the error rate refers to the whole family of tests rather than to each one separately.

What multiple testing corrections control

Two error rates are common. The family-wise error rate (FWER) is the probability of rejecting at least one true null hypothesis. The false discovery rate (FDR) is the expected share of false rejections among all rejections. Controlling the FWER is stricter: it is appropriate when any single false claim is costly, such as a headline result of a policy evaluation. Controlling the FDR accepts some false discoveries in exchange for power, which suits screening many hypotheses.

• Bonferroni multiplies each p-value by the number of tests m. It controls the FWER under any dependence between the tests, but it is conservative.
• Holm applies the same idea step by step. It is never less powerful than Bonferroni and makes no assumption about dependence. With the p-values sorted from smallest to largest, the adjusted value of the i-th is

p̃₍ᵢ₎ = max over j ≤ i of min{1, (m − j + 1) · p₍ⱼ₎}

• Šidák, Holm–Šidák, Hochberg and Hommel are slightly more powerful, but each needs independence or a specific dependence condition.
• Benjamini–Hochberg (BH) controls the FDR when the tests are independent (Benjamini and Hochberg 1995) or positively dependent (Benjamini and Yekutieli 2001). Benjamini–Yekutieli (BY) controls it under any dependence, at a cost in power (Benjamini and Yekutieli 2001).
• Romano–Wolf controls the FWER using the joint distribution of the test statistics (Romano and Wolf 2005). When the statistics are correlated, as with related outcomes measured on the same people, the largest of them is less extreme than it would be for independent tests, so the procedure can reject more hypotheses than Holm while keeping the same guarantee. It needs draws of all the statistics under the null, typically from a bootstrap or a re-randomization.

List, Shaikh and Xu (2019) discuss these issues for experiments with several outcomes, subgroups and treatments. Whatever the method, the family has to be fixed before the results are seen. Choosing which tests to correct after looking at the p-values defeats the purpose.

In OpenEconometrics

oe.multipletests(p_values, method="holm", alpha=0.05, labels=None) adjusts a list of p-values. method is "bonferroni", "sidak", "holm", "holm_sidak", "hochberg", "hommel", "bh" or "by". The result is a table with the hypothesis, the raw p_value, the adjusted_p_value and reject, in the order the p-values were given; a hypothesis is rejected when its adjusted p-value is at most alpha. Missing p-values raise an error unless missing="drop" excludes them explicitly. The table's attrs record the family, the method and its dependence assumption.

oe.stepdown(observed, null_draws, method="romano_wolf", null_description=...) runs the Romano–Wolf stepdown. observed holds the m test statistics and null_draws a B × m matrix of the same statistics computed under the null. The function does not generate the draws: you supply them from a design you can justify, and null_description records what that design is. tail is "two-sided", "greater" or "less". calibration="monte_carlo" (the default) uses (1 + count)/(B + 1); calibration="enumerated" is for draws that list every equally likely assignment. method="westfall_young" runs the minP version on supplied p-values. The table adds the raw and adjusted p-values to the statistics.

For a classical OLS fit, oe.ols_stepdown(fit, terms=[...], null_values=[...], error_model="iid_gaussian") generates the joint null itself under an explicitly declared Gaussian fixed-design model, and oe.simultaneous_ci builds simultaneous confidence intervals from estimates and their joint covariance. In Stata, the community-contributed rwolf command computes Romano–Wolf p-values by bootstrapping regressions. The multiple-testing guide states each method's formula and assumptions.

Supported options

oe.multipletests

Stata command: none listed.

oe.stepdown

Stata command: none listed.

The inventory section multiple_testing for oe.multipletests and oe.stepdown lists adjustments: bonferroni, sidak, holm, holm_sidak, hochberg, hommel, bh, by; joint null: romano_wolf_maxT, westfall_young_minP.

Generated from oe.capabilities() of openecon 0.3.20a1. Stata names identify the corresponding workflow, not numerical parity.

Example

The script simulates a randomized experiment with 400 participants, exactly 200 of them treated, and ten outcomes measured in standard-deviation units. The outcomes share a common factor, so any two have a correlation of 0.5. Treatment raises the first three outcomes by 0.25 standard deviations and leaves the other seven unchanged. Each outcome is tested with a Welch t statistic. The null draws come from 1,999 re-randomizations of the treatment labels: for an outcome that treatment does not affect, re-assigning the labels reproduces the distribution of its statistic, jointly with the other unaffected outcomes.

import numpy as np
import pandas as pd
import openecon as oe


def simulate(rng, n=400, m=10):
    """An experiment with ten correlated outcomes; treatment moves only the first three."""
    effects = np.array([0.25, 0.25, 0.25] + [0.0] * (m - 3))
    treat = rng.permutation(np.repeat([1, 0], n // 2))  # 200 of 400 treated
    common = rng.normal(size=(n, 1))
    outcomes = np.sqrt(0.5) * common + np.sqrt(0.5) * rng.normal(size=(n, m))  # correlation 0.5
    return treat, outcomes + np.outer(treat, effects), effects


def welch_t(treat, y):
    """Difference in means divided by its unequal-variance standard error, one per outcome."""
    a, b = y[treat == 1], y[treat == 0]
    se = np.sqrt(a.var(axis=0, ddof=1) / len(a) + b.var(axis=0, ddof=1) / len(b))
    return (a.mean(axis=0) - b.mean(axis=0)) / se


if __name__ == "__main__":
    rng = np.random.default_rng(31)
    treat, y, effects = simulate(rng)
    names = [f"y{j + 1}" for j in range(y.shape[1])]
    observed = welch_t(treat, y)
    null_draws = np.array([welch_t(rng.permutation(treat), y) for _ in range(1999)])
    rw = oe.stepdown(observed, null_draws, method="romano_wolf", labels=names,
                     null_description="1,999 re-randomizations of the treatment labels")
    table = pd.DataFrame({"effect": effects, "t": observed,
                          "raw p": rw["p_value"].to_numpy()}, index=names)
    for method in ["bonferroni", "holm", "bh"]:
        adjusted = oe.multipletests(list(rw["p_value"]), method=method, labels=names)
        table[method] = list(adjusted["adjusted_p_value"])
    table["romano_wolf"] = list(rw["adjusted_p_value"])
    print(table.round(3).to_string())
    print("Rejections at 5%:", {col: int((table[col] <= 0.05).sum()) for col in table.columns[2:]})

With openecon 0.3.20a1 from PyPI the script prints:

     effect      t  raw p  bonferroni   holm     bh  romano_wolf
y1     0.25  3.550  0.001       0.010  0.010  0.010        0.002
y2     0.25  2.110  0.034       0.345  0.276  0.115        0.182
y3     0.25  3.031  0.004       0.040  0.036  0.020        0.019
y4     0.00  0.884  0.358       1.000  1.000  0.598        0.832
y5     0.00 -0.507  0.621       1.000  1.000  0.776        0.945
y6     0.00  0.983  0.342       1.000  1.000  0.598        0.831
y7     0.00  0.010  0.994       1.000  1.000  0.994        0.994
y8     0.00 -0.574  0.580       1.000  1.000  0.776        0.945
y9     0.00 -0.389  0.706       1.000  1.000  0.785        0.945
y10    0.00  1.861  0.060       0.605  0.424  0.151        0.276
Rejections at 5%: {'raw p': 3, 'bonferroni': 2, 'holm': 2, 'bh': 2, 'romano_wolf': 2}

Reading the output

• Without correction, three outcomes are significant. The raw re-randomization p-values are 0.001, 0.034 and 0.004 for the three affected outcomes, and 0.060 for the unaffected y10. In this sample the uncorrected tests happen to make no false rejection, but with seven true nulls that is luck: if the seven tests were independent, the chance of at least one false rejection would be 1 − 0.95⁷ ≈ 30%.
• Holm improves on Bonferroni, slightly. Both reject y1 and y3. Holm's adjusted p-value for y3 is 0.036 against Bonferroni's 0.040, because Holm multiplies the second-smallest p-value by 9 rather than 10.
• Benjamini–Hochberg is less strict but answers a different question. It adjusts y2 to 0.115 and y10 to 0.151. At an FDR of 5% it still rejects only y1 and y3; a 12% FDR would add y2, and y10 would follow only above 15%.
• Romano–Wolf uses the correlation between the outcomes. Its adjusted p-values for the affected outcomes, 0.002, 0.182 and 0.019, are well below Holm's 0.010, 0.276 and 0.036. With correlated outcomes, the largest of ten t statistics under the null is less extreme than for ten independent tests, so the critical value is lower.
• Every correction misses y2. Its true effect is as large as the others, but its t statistic is only 2.11. Controlling the error rate across ten tests costs power, and with 200 participants per arm an effect of 0.25 standard deviations is not always detected.

Limitations

• oe.multipletests and oe.stepdown work on in-memory data only.
• oe.stepdown does not generate null draws. The validity of the result rests on the supplied draws, and the table's attrs record null_design_verified=False; a plus-one Monte Carlo correction cannot make an invalid bootstrap valid.
• Strong FWER control by Romano–Wolf requires the null draws to approximate every subset of true nulls, with monotone critical values; Westfall–Young requires valid marginal null distributions and subset pivotality, and its tail must stay "two-sided".
• Šidák, Holm–Šidák, Hochberg, Hommel and Benjamini–Hochberg rely on independence or specific dependence conditions that the function cannot check; only Bonferroni, Holm and Benjamini–Yekutieli hold under arbitrary dependence.
• missing="drop" changes the tested family; it is not a remedy for outcomes that are missing because of their values.
• oe.ols_stepdown covers classical unweighted or known-precision weighted OLS only; robust, cluster and HAC covariances and model-specific bootstraps are not supported.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Benjamini, Y. and Hochberg, Y. (1995). Controlling the false discovery rate: A practical and powerful approach to multiple testing. Journal of the Royal Statistical Society: Series B (Methodological), 57(1), 289–300. doi:10.1111/j.2517-6161.1995.tb02031.x
• Benjamini, Y. and Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under dependency. Annals of Statistics, 29(4), 1165–1188. doi:10.1214/aos/1013699998
• List, J. A., Shaikh, A. M. and Xu, Y. (2019). Multiple hypothesis testing in experimental economics. Experimental Economics, 22(4), 773–793. doi:10.1007/s10683-018-09597-5
• Romano, J. P. and Wolf, M. (2005). Exact and approximate stepdown methods for multiple hypothesis testing. Journal of the American Statistical Association, 100(469), 94–108. doi:10.1198/016214504000000539