Multiplicity and Multiple Comparisons

In plain language

When a study tests many different outcomes or research questions at once, the chance of finding at least one false positive just by chance grows rapidly — running 20 tests each at the standard 5% threshold gives roughly a 64% chance of at least one spurious result even when nothing is truly different. Multiplicity adjustment is a set of statistical rules for raising the bar when many tests are run together, so the overall false-alarm rate stays controlled. The two main schools are controlling the family-wise error rate (keeping the chance of any false positive below 5%, using Bonferroni or Holm corrections, required for regulatory submissions) and controlling the false discovery rate (keeping the expected share of falsely flagged results below 5%, using Benjamini-Hochberg, suited for exploratory outcome scans). The right choice depends on whether the analysis is confirmatory or exploratory, and the most important protection of all is writing down exactly which tests will be run — and how they will be adjusted — before the data are analyzed.

The multiplicity problem in RWE

Every additional hypothesis test gives a study one more chance to declare a false positive by chance alone. At a 5% significance threshold and one test, the type-I error rate is 5%. With ten independent tests all truly null, the probability of at least one spurious rejection rises to 1 - 0.95^10, which is approximately 40%. At 100 truly null tests it approaches 99%. Real-world evidence studies routinely generate exactly this inflation: secondary endpoint batteries of 10-30 outcomes, subgroup analyses across age, sex, and comorbidity strata, multiple sensitivity analyses, and the silent flexibility of choosing which results to highlight after seeing the data. The discipline of multiplicity adjustment provides a coherent framework for deciding — before looking at the data — how to allocate the permitted false-positive rate across the tests being run.

Two philosophies: family-wise error rate versus false discovery rate

The central conceptual divide is between two definitions of what "error" means when many tests are run simultaneously.

The family-wise error rate (FWER) is the probability of making at least one false positive rejection in the entire family of hypotheses. This is the standard for confirmatory clinical trials and regulatory submissions: the FDA and EMA require that the chance of erroneously declaring any efficacy endpoint effective is controlled at alpha = 0.05 across the full testing plan. FWER is the right criterion when the consequence of any single wrong conclusion is unacceptable — for example, approving a drug on a secondary endpoint that is actually null.

The false discovery rate (FDR) is the expected proportion of rejected hypotheses that are false positives. If a study rejects 20 hypotheses and the FDR is controlled at 5%, the analyst expects at most 1 of those 20 rejections to be wrong on average across many hypothetical replications of the same procedure. FDR is the right criterion for exploratory science — pharmacovigilance scans of thousands of drug-event pairs, phenome-wide association studies across thousands of disease codes, or metabolomics screens — where the goal is to nominate a tractable shortlist of hypotheses for follow-up, not to make a definitive claim about any single one. Accepting a few false leads in exchange for better sensitivity is a reasonable trade-off when the rejected hypotheses will be independently validated.

The conceptual consequence: FWER control becomes harder as the number of tests m grows — Bonferroni sets the per-test threshold at alpha/m, which approaches zero. FDR control via Benjamini-Hochberg scales more gracefully; the threshold for the k-th-ranked test is k/(m*Q), which remains non-trivial even for large m. For hypothesis-free scans across hundreds of outcomes — where Bonferroni would reject almost nothing — FDR or permutation-based FWER (as in TreeScan) are the appropriate tools.

Bonferroni and Holm (FWER methods)

The Bonferroni correction is the simplest FWER method: divide alpha by the total number of tests m and reject hypothesis i if its raw p-value p_i is at most alpha/m. With ten tests and alpha = 0.05, the Bonferroni threshold is 0.005. Bonferroni is valid under any correlation structure among the tests — it relies only on the union bound (the Boole inequality) and does not assume independence. Its weakness is conservatism: when the tests are positively correlated (as they often are across related clinical outcomes or overlapping diagnostic code families), Bonferroni over-adjusts, and many genuine signals are missed.

The Holm step-down procedure uniformly improves on Bonferroni while still controlling FWER. Sort p-values ascending p_(1) ≤ ... ≤ p_(m); compare p_(i) to alpha/(m - i + 1) starting from rank 1; continue rejecting as long as the condition holds; stop at the first failure. Holm is always at least as powerful as Bonferroni and often substantially more so when the smallest p-values are far below the others. Holm should be the default wherever Bonferroni is considered. Both Bonferroni and Holm are conservative when tests are positively correlated; procedures that account for correlation (e.g., permutation-based FWER, the Westfall-Young method) exist but require access to the joint distribution of the test statistics.

Benjamini-Hochberg step-up (FDR method)

The Benjamini-Hochberg (BH) procedure (1995) is the workhorse FDR method. Sort p-values ascending p_(1) ≤ ... ≤ p_(m). For each rank i, compute the BH threshold (i/m) Q where Q is the desired FDR level (typically 0.05 or 0.10). Find the largest rank k such that p_(k) ≤ (k/m) Q. Reject all hypotheses with rank 1 through k. This is a step-up procedure: it looks for the largest k where the condition holds and rejects everything at or below that rank.

BH provably controls FDR at level Q * (true nulls / m) ≤ Q under independent tests, and under certain positive-dependence structures (positive regression dependence on a subset, or PRDS). Under arbitrary dependence, the Benjamini-Yekutieli (BY) extension controls FDR at the cost of dividing each threshold by an additional log(m) factor. In most claims and EHR applications, where outcomes are positively correlated (patients with one cardiovascular event tend to have elevated risk of others), the standard BH procedure is conservative relative to the true FDR, making it a defensible choice.

When each method fits in RWE

The decision maps onto the confirmatory-versus-exploratory spectrum:

• Few pre-specified confirmatory endpoints (1-5): use Holm or Bonferroni, because each
• Moderate battery of secondary endpoints (5-30): if each is pre-specified and the intention
• Hypothesis-free scans (hundreds to thousands of outcomes or codes): FDR via BH or BY, or

Gatekeeping and hierarchical testing

When endpoints have a natural priority order — overall survival is the primary endpoint and progression-free survival is secondary — gatekeeping procedures allocate alpha in a structured hierarchy: the secondary endpoint is only tested if the primary endpoint passes. This controls FWER by construction because no alpha is spent on the secondary if the primary fails. Variants include fixed-sequence testing (strict ordering), parallel gatekeeping (multiple primaries must pass before any secondary is tested), and tree gatekeeping (a hierarchy of families). These procedures are common in oncology and cardiovascular regulatory submissions, where the ordering reflects clinical importance and partial approvals have well-defined regulatory implications.

Subgroup multiplicity

Every subgroup analysis is an additional test in the multiplicity family. If an analyst tests a treatment comparison in males, females, the elderly, and those with baseline comorbidity, the chance of finding a spurious interaction in at least one subgroup at alpha = 0.05 is approximately 19%, not 5%. Standard pharmacoepidemiology practice treats subgroup analyses as exploratory and hypothesis-generating unless pre-specified in the protocol with an explicit correction plan. Route to subgroup-analysis-hte for the full treatment of heterogeneous treatment effect estimation and the multiplicity implications of multiple subgroup comparisons.

Sequential-look multiplicity

Repeated interim analyses of accumulating data are a temporal form of the multiple-comparisons problem. Looking at the data every month and applying a fixed alpha = 0.05 each time inflates the family-wise false-alarm rate just as badly as running many simultaneous hypothesis tests. Sequential analysis via alpha-spending (MaxSPRT, O'Brien-Fleming, Lan-DeMets) provides analogous boundaries: the total probability of ever falsely signaling across all planned looks is controlled at the pre-specified alpha. The boundary is computed before the first look and governs every subsequent one. See maxsprt-sequential-safety-surveillance-rwe for the full mechanics of sequential alpha spending in pharmacovigilance contexts.

The garden of forking paths: silent multiplicity

The most insidious multiplicity in observational RWE is the kind that never appears in the analysis plan: the analyst who reports the washout period that gave the "cleanest" result, the covariate set selected after seeing the coefficient, the subgroup noted because the p-value looked interesting. Each decision that could have gone differently — and that is made after the analyst has some knowledge of the data — creates a silent test that is never counted in the multiplicity denominator. Gelman and Loken (2014) called this the "garden of forking paths." Post-hoc Bonferroni or BH adjustment applied to a selectively reported set of results cannot remedy multiplicity that was incurred during data exploration, because the denominator m reflects only reported tests, not all considered paths. The only real fix is prespecification: an SAP written and locked before database access or unblinding, with a complete multiplicity plan covering the primary test, all secondary tests, and the rules for subgroup and sensitivity analyses.

Large-n datasets and the limits of multiplicity adjustment

In very large claims databases with millions of patients, virtually every test is statistically significant at any conventional alpha. Adjusting p-values for multiplicity does not rescue inference when the root problem is that p-values are meaninglessly small for clinically trivial effects. A hazard ratio of 1.003 with a Bonferroni-adjusted p-value of 0.0001 is not a meaningful finding. The discipline of multiplicity adjustment was developed to protect against false positives when true signal is sparse relative to noise; it does not substitute for interpreting effect sizes and clinical relevance. Report hazard ratios, risk differences, and confidence intervals alongside any adjusted p-value — especially at the large sample sizes typical in administrative claims research.

Interpreting the output

From the worked example: 10 secondary endpoints tested in a statin initiation cohort, Bonferroni threshold alpha/m = 0.05/10 = 0.005, and BH at Q = 0.05. Bonferroni rejects 1 hypothesis (H01 MI, p = 0.001). BH step-up finds the largest rank k where p_(k) ≤ (k/10)0.05: rank 1 threshold 1/100.05 = 0.005 (H01 MI, p=0.001, passes); rank 2 threshold 2/100.05 = 0.010 (H03 stroke, p=0.008, passes); rank 3 threshold 3/100.05 = 0.015 (H05 hospitalization, p=0.012, passes); rank 4 threshold 4/10*0.05 = 0.020 (H07 ER visit, p=0.025, fails). BH rejects 3 hypotheses: H01, H03, H05.

(1) Formal interpretation. The Bonferroni procedure rejects H01 (MI, p = 0.001) because 0.001 ≤ 0.005 = alpha/m. Under the Boole inequality, the family-wise error rate across all 10 tests is controlled at ≤ 0.05 regardless of the correlation structure among the tests. H03 (stroke, p = 0.008) is not rejected: although p = 0.008 is below the unadjusted alpha = 0.05, it exceeds the Bonferroni threshold. The BH procedure rejects H01, H03, and H05 by identifying k = 3 as the largest rank satisfying p_(k) ≤ (k/m)*Q. The FDR guarantee for these 3 rejections is an expected bound: across many hypothetical replications of the same study and procedure, the expected proportion of false positives among the rejected hypotheses is at most Q = 0.05. This is a statement about the long-run average over replications, not a claim that exactly 5% of these three specific rejections are false.

(2) Practical interpretation. For a regulatory submission where MI is the confirmatory primary endpoint, Bonferroni is appropriate: the chance of any false confirmatory claim is controlled at 5%, and the stroke result (p = 0.008) is appropriately classified as exploratory. For an internal clinical evidence review — deciding which outcomes to carry into further study — BH is more appropriate: MI, stroke, and hospitalization are nominated together as a cluster warranting independent confirmatory analysis, with the understanding that roughly 5% of nominated findings are expected to be false leads. The analyst must communicate clearly which procedure governed the primary inference and treat the other procedure's results as supportive or exploratory.

Pros, cons, and trade-offs

Bonferroni/Holm (FWER): - Pros: valid under any correlation structure; conservative and transparent; required by FDA/EMA for confirmatory analyses; straightforward to apply and communicate. - Cons: increasingly conservative as m grows; at m > 100, almost no hypothesis survives; ignores positive correlation among tests (which could safely allow a less strict threshold). - When to prefer: confirmatory regulatory submissions; few pre-specified primary endpoints; settings where any single false positive has serious consequences (approval, label change).

Benjamini-Hochberg (FDR): - Pros: substantially more powerful than Bonferroni at the same nominal error level; scales gracefully to hundreds or thousands of tests; standard for genomics, metabolomics, and pharmacovigilance outcome scans. - Cons: controls an expected average error, not the probability of any false positive — in any single study more than Q fraction of rejections could be false; assumes independence or PRDS (BY extension handles arbitrary dependence at a power cost). - When to prefer: hypothesis-generating screens; pharmacovigilance scans; phenome-wide studies; any setting where follow-up validation is planned and some false leads are acceptable.

Hierarchical/gatekeeping procedures: - Pros: FWER-controlling; structure testing to match clinical priority; allows secondary endpoints to be tested without alpha inflation when the primary succeeds. - Cons: if the primary fails, all alpha is spent and no secondary information is extracted; requires pre-specification of the full hierarchy before any data access. - When to prefer: regulatory submissions with a primary-secondary endpoint hierarchy; oncology trials with OS as gatekeeper for PFS or OS, PFS, and QoL in a fixed sequence.

When to use

Apply formal multiplicity adjustment whenever: - Two or more hypotheses are tested in a single analysis and both will be reported as confirmatory findings, not exploratory observations. - A regulatory or HTA body requires control of the family-wise error rate across the study's testing plan. - A pharmacovigilance or outcomes scan tests hundreds of outcomes simultaneously — use FDR rather than FWER; consider permutation-based methods (TreeScan) when outcomes are hierarchically structured. - An interim analysis plan involves more than one planned look at accumulating data — use sequential alpha-spending, not repeated fixed-alpha tests. - Subgroup analyses are pre-specified and must be reported alongside the primary analysis: apply a correction (Holm across subgroups) or commit in the protocol to treating them as exploratory.

When NOT to use — and when adjustment is actively misleading

• Do not apply multiplicity adjustment post-hoc to a selectively reported set of outcomes: if 30
• Do not treat multiplicity adjustment as a substitute for prespecification: an analyst who runs
• Do not apply FWER correction to purely exploratory analyses intended to generate hypotheses:
• Do not interpret a multiplicity-adjusted p-value as evidence of clinical importance: at very
• Do not confuse multiplicity adjustment with confounding control: adjusting p-values does not

Worked example

Scenario

An observational claims study of statin initiation tests 10 outcomes simultaneously — myocardial infarction, stroke, all-cause hospitalization, ER visit, LDL reduction, CRP reduction, QOL score, medication switch, any adverse event, and all-cause mortality — all pre-specified as secondary endpoints. The analyst runs 10 separate Cox models and collects the raw p-values. She applies both Bonferroni and Benjamini-Hochberg corrections at alpha/Q = 0.05 and compares how many hypotheses each procedure rejects.

Dataset

Raw p-values from 10 Cox models (one per outcome) in a statin initiation observational cohort. P-values appear in their original analysis order; multiplicity procedures sort them internally by rank.

Steps

• Sort the 10 raw p-values from smallest to largest to assign ranks 1 through 10. Sorted order: H01=0.001 (rank 1), H03=0.008 (rank 2), H05=0.012 (rank 3), H07=0.025 (rank 4), H02=0.040 (rank 5), H10=0.051 (rank 6), H09=0.120 (rank 7), H04=0.200 (rank 8), H08=0.350 (rank 9), H06=0.500 (rank 10). Label these p_(1) through p_(10).

• Bonferroni correction: the adjusted threshold is 0.05/10 = 0.005. Compare each raw p-value to 0.005. Only H01 (MI, p=0.001) satisfies p at most 0.005. All others, including H03 (stroke, p=0.008), exceed the threshold. Bonferroni rejects exactly 1 hypothesis.

• Benjamini-Hochberg step-up: for each rank i, the BH threshold is i/100.05. Compute thresholds for ranks 1 through 4. Rank 1: 1/100.05 = 0.005. Rank 2: 2/100.05 = 0.010. Rank 3: 3/100.05 = 0.015. Rank 4: 4/10*0.05 = 0.020.

• Compare each sorted p-value to its BH threshold going upward from rank 1. p_(1) = 0.001 is at most 0.005 (rank 1 passes). p_(2) = 0.008 is at most 0.010 (rank 2 passes). p_(3) = 0.012 is at most 0.015 (rank 3 passes). p_(4) = 0.025 exceeds 0.020 (rank 4 fails, first failure going upward).

• BH step-up rule: find the largest rank k such that p_(k) is at most k/10*0.05, then reject all hypotheses with rank 1 through k. The largest such k is 3, because rank 3 passes and rank 4 fails. Reject H01 (MI), H03 (stroke), H05 (hospitalization). BH rejects 3 hypotheses.

Result

Bonferroni threshold: 0.05/10 = 0.005. Only H01 (MI, p=0.001) satisfies p at most 0.005; Bonferroni rejects 1 hypothesis. BH step-up: rank 1 threshold 1/100.05 = 0.005 (H01 MI p=0.001 passes), rank 2 threshold 2/100.05 = 0.010 (H03 stroke p=0.008 passes), rank 3 threshold 3/100.05 = 0.015 (H05 hospitalization p=0.012 passes), rank 4 threshold 4/100.05 = 0.020 (H07 ER visit p=0.025 fails, first failure). Largest passing rank is 3; BH rejects 3 hypotheses (H01 MI, H03 stroke, H05 hospitalization) versus 1 under Bonferroni.

Runnable example

from statsmodels.stats.multitest import multipletests
import pandas as pd

# ── Worked-example dataset: 10 raw p-values from statin initiation cohort ──
outcomes = [
    "MI", "LDL_reduction", "stroke", "medication_switch",
    "hospitalization", "all_cause_mortality", "ER_visit",
    "any_adverse_event", "QOL_score", "CRP_reduction",
]
p_raw = [0.001, 0.040, 0.008, 0.200, 0.012, 0.500, 0.025, 0.350, 0.120, 0.051]

# ── 1. Bonferroni (FWER): threshold = alpha / m = 0.05 / 10 = 0.005 ──
rej_bon, p_bon, _, _ = multipletests(p_raw, alpha=0.05, method="bonferroni")

# ── 2. Holm step-down (FWER; uniformly more powerful than Bonferroni) ──
rej_holm, p_holm, _, _ = multipletests(p_raw, alpha=0.05, method="holm")

# ── 3. Benjamini-Hochberg step-up (FDR at Q = 0.05) ──
rej_bh, p_bh, _, _ = multipletests(p_raw, alpha=0.05, method="fdr_bh")

# ── 4. Display results ──
df = pd.DataFrame({
    "outcome": outcomes,
    "raw_p":   p_raw,
    "bon_p":   [round(p, 4) for p in p_bon],
    "holm_p":  [round(p, 4) for p in p_holm],
    "bh_p":    [round(p, 4) for p in p_bh],
    "bon_rej": rej_bon,
    "holm_rej": rej_holm,
    "bh_rej":  rej_bh,
})
print("Multiplicity adjustment — statin cohort 10-outcome battery")
print(df.to_string(index=False))
print(f"\nBonferroni rejects: {rej_bon.sum()} / 10")
print(f"Holm rejects:       {rej_holm.sum()} / 10")
print(f"BH (FDR) rejects:   {rej_bh.sum()} / 10")

# ── 5. Manual BH step-up for transparency ──
# Sort indices by raw p-value
m, Q = 10, 0.05
sorted_idx = sorted(range(m), key=lambda i: p_raw[i])
print("\nManual BH step-up (rank, outcome, raw_p, threshold, pass):")
k_max = 0
for rank, idx in enumerate(sorted_idx, start=1):
    thresh = rank / m * Q   # i/m * Q
    passed = p_raw[idx] <= thresh
    if passed:
        k_max = rank
    print(f"  rank {rank}: {outcomes[idx]:25s} raw_p={p_raw[idx]:.3f}  "
          f"thresh={thresh:.4f}  {'PASS' if passed else 'FAIL'}")
print(f"Largest passing rank k = {k_max}; reject ranks 1-{k_max}.")
reject_manual = {outcomes[sorted_idx[i]] for i in range(k_max)}
print(f"Rejected outcomes: {sorted(reject_manual)}")

# ── Worked-example dataset: 10 raw p-values from statin initiation cohort ──
outcomes <- c("MI", "LDL_reduction", "stroke", "medication_switch",
              "hospitalization", "all_cause_mortality", "ER_visit",
              "any_adverse_event", "QOL_score", "CRP_reduction")
p_raw <- c(0.001, 0.040, 0.008, 0.200, 0.012, 0.500, 0.025, 0.350, 0.120, 0.051)

# ── 1. Bonferroni (FWER): each adjusted p = min(m * raw_p, 1) ──
p_bon  <- p.adjust(p_raw, method = "bonferroni")

# ── 2. Holm step-down (FWER; always at least as powerful as Bonferroni) ──
p_holm <- p.adjust(p_raw, method = "holm")

# ── 3. Benjamini-Hochberg step-up (FDR at Q = 0.05) ──
p_bh   <- p.adjust(p_raw, method = "BH")

# ── 4. Benjamini-Yekutieli (FDR under arbitrary dependence; more conservative) ──
p_by   <- p.adjust(p_raw, method = "BY")

# ── 5. Collect results and display ──
results <- data.frame(
  outcome   = outcomes,
  raw_p     = p_raw,
  bon_adj   = round(p_bon,  4),
  holm_adj  = round(p_holm, 4),
  bh_adj    = round(p_bh,   4),
  by_adj    = round(p_by,   4),
  bon_rej   = p_bon  <= 0.05,
  holm_rej  = p_holm <= 0.05,
  bh_rej    = p_bh   <= 0.05,
  by_rej    = p_by   <= 0.05
)
print(results)
cat(sprintf("\nBonferroni rejects: %d / 10\n", sum(p_bon  <= 0.05)))
cat(sprintf("Holm rejects:       %d / 10\n", sum(p_holm <= 0.05)))
cat(sprintf("BH (FDR) rejects:   %d / 10\n", sum(p_bh   <= 0.05)))
cat(sprintf("BY (FDR) rejects:   %d / 10\n", sum(p_by   <= 0.05)))

# ── 6. Manual BH step-up for transparency ──
m <- 10L; Q <- 0.05
ord       <- order(p_raw)               # ascending rank indices
p_sorted  <- p_raw[ord]
thresholds <- seq_len(m) / m * Q        # i/m * Q for i = 1 ... m
pass       <- p_sorted <= thresholds
k_max      <- max(which(pass), 0)       # largest rank where condition holds

cat("\nManual BH step-up:\n")
for (i in seq_len(m)) {
  cat(sprintf("  rank %2d: %-25s raw_p=%.3f  thresh=%.4f  %s\n",
              i, outcomes[ord[i]], p_sorted[i], thresholds[i],
              if (pass[i]) "PASS" else "FAIL"))
}
cat(sprintf("Largest passing rank k = %d; reject ranks 1 to %d.\n", k_max, k_max))
cat(sprintf("BH-rejected outcomes: %s\n",
            paste(outcomes[ord[seq_len(k_max)]], collapse = ", ")))