Multiplicity and Multiple Comparisons
A statistical discipline governing how to set significance thresholds and interpret p-values when a study tests multiple hypotheses simultaneously — whether across several endpoints, subgroups, interim looks, or unreported analytical choices — so that the overall false-positive rate (or false-discovery rate) is controlled at a meaningful level rather than inflated by the accumulation of chances across many tests; the two main tools are family-wise error rate (FWER) control via Bonferroni or Holm for confirmatory analyses and false discovery rate (FDR) control via Benjamini-Hochberg for exploratory scans.
On this page
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 endpoint is expected to survive as a confirmed finding if the null is rejected. Regulatory submissions require FWER control across the confirmatory testing plan.
- Moderate battery of secondary endpoints (5-30): if each is pre-specified and the intention is to publish results for each, prefer Holm; if the battery is exploratory and will be followed up, BH is more appropriate. The critical requirement is prespecification: write the adjustment plan into the protocol or statistical analysis plan before any data access.
- Hypothesis-free scans (hundreds to thousands of outcomes or codes): FDR via BH or BY, or permutation-based FWER via scan statistics (TreeScan). Bonferroni over thousands of correlated codes is so conservative as to be uninformative.
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/10*0.05 = 0.005 (H01 MI, p=0.001, passes); rank 2 threshold 2/10*0.05 = 0.010 (H03 stroke, p=0.008, passes); rank 3 threshold 3/10*0.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 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 outcomes were analyzed and only 10 with "interesting" results are presented, the multiplicity denominator is wrong, and the apparent correction provides false reassurance. The unreported 20 are still part of the implicit testing context.
- Do not treat multiplicity adjustment as a substitute for prespecification: an analyst who runs 50 models and then applies BH to the 10 most significant results has not controlled FDR — the unseen 40 are part of the effective multiple-testing context.
- Do not apply FWER correction to purely exploratory analyses intended to generate hypotheses: doing so creates false confidence that the surviving hypotheses are confirmed when they are actually just the strongest leads from a screen.
- Do not interpret a multiplicity-adjusted p-value as evidence of clinical importance: at very large n, adjusted p-values remain significant for effect sizes that are clinically irrelevant. Report effect sizes and confidence intervals as the primary evidence and treat any adjusted p-value as a decision aid, not a summary of importance.
- Do not confuse multiplicity adjustment with confounding control: adjusting p-values does not correct for uncontrolled confounding in the underlying effect estimates — these are separate problems requiring separate solutions.
Decision diagram
flowchart TD Q[How many hypotheses?<br/>What is the scientific goal?] --> Conf[Few endpoints<br/>1-5, confirmatory] Q --> Expl[Moderate battery<br/>5-30, exploratory] Q --> Scan[Many outcomes<br/>>30, hypothesis-free scan] Conf --> Ordered[Ordered by<br/>clinical priority?] Ordered -->|Yes| Gate["Gatekeeping / fixed-sequence<br/>(test primary first; secondary<br/>only if primary passes)"] Ordered -->|No| Holm["Holm step-down<br/>(default FWER;<br/>more powerful than Bonferroni)"] Expl --> BH["Benjamini-Hochberg<br/>step-up at Q = 0.05-0.10<br/>(FDR; label as exploratory)"] Scan --> Hier[Outcomes hierarchically<br/>structured ICD/MedDRA?] Hier -->|Yes| Tree["TreeScan permutation FWER<br/>(see tree-based-scan-statistic-rwe)"] Hier -->|No| BH2["Benjamini-Hochberg FDR<br/>or Benjamini-Yekutieli<br/>(for arbitrary dependence)"] Q --> Seq[Sequential interim<br/>looks over time?] Seq --> MaxS["MaxSPRT / alpha-spending<br/>(see maxsprt-sequential-safety-surveillance-rwe)"]
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.
| hypothesis | outcome | raw_p_value |
|---|---|---|
| H01 | myocardial_infarction | 0.001 |
| H02 | LDL_reduction | 0.04 |
| H03 | stroke | 0.008 |
| H04 | medication_switch | 0.2 |
| H05 | hospitalization | 0.012 |
| H06 | all_cause_mortality | 0.5 |
| H07 | ER_visit | 0.025 |
| H08 | any_adverse_event | 0.35 |
| H09 | QOL_score | 0.12 |
| H10 | CRP_reduction | 0.051 |
Steps
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/10*0.05 = 0.005 (H01 MI p=0.001 passes), rank 2 threshold 2/10*0.05 = 0.010 (H03 stroke p=0.008 passes), rank 3 threshold 3/10*0.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, first failure). Largest passing rank is 3; BH rejects 3 hypotheses (H01 MI, H03 stroke, H05 hospitalization) versus 1 under Bonferroni.
Trade-offs
Runnable example
Bonferroni, Holm step-down, and Benjamini-Hochberg FDR adjustment using statsmodels.stats.multitest.multipletests. Demonstrates all three methods on the 10-outcome statin cohort dataset from the worked example. Also shows manual BH step-up computation for pedagogical transparency. No dependencies beyond statsmodels.
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)}")Bonferroni, Holm, and Benjamini-Hochberg adjustments via base-R p.adjust(). Demonstrates all three methods on the 10-outcome statin cohort dataset. Also shows the Benjamini-Yekutieli (BY) adjustment for arbitrary dependence. No external packages required.
# ── 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 = ", ")))Bonferroni, Holm, and Benjamini-Hochberg adjustment via PROC MULTTEST with the INPVALUES option, which accepts a SAS dataset of pre-computed raw p-values. Also demonstrates a DATA step manual Bonferroni check. PROC MULTTEST is the SAS standard for multiplicity adjustment in submissions and is accepted by FDA/EMA.
/* ── Step 1: Create dataset of raw p-values from 10 Cox models ── */
data work.raw_p;
length hypothesis $5 outcome $25;
input hypothesis $ outcome $ raw_p;
datalines;
H01 MI 0.001
H02 LDL_reduction 0.040
H03 stroke 0.008
H04 medication_switch 0.200
H05 hospitalization 0.012
H06 all_cause_mortality 0.500
H07 ER_visit 0.025
H08 any_adverse_event 0.350
H09 QOL_score 0.120
H10 CRP_reduction 0.051
;
run;
/* ── Step 2: PROC MULTTEST — read raw p-values and apply adjustments ──
inpvalues(raw_p): the variable "raw_p" in work.raw_p contains the raw p-values.
bon: Bonferroni — adjusted p = min(m * raw_p, 1); controls FWER.
holm: Holm step-down — uniformly more powerful than Bonferroni; controls FWER.
fdr: Benjamini-Hochberg — step-up procedure; controls FDR at Q = 0.05 by default.
out=: stores adjusted p-values alongside the input dataset variables. */
proc multtest inpvalues(raw_p)=work.raw_p
bon
holm
fdr
out=work.adjusted;
run;
/* ── Step 3: Merge labels back and print ── */
data work.results;
merge work.raw_p work.adjusted;
run;
proc print data=work.results noobs;
var hypothesis outcome raw_p;
title "Multiplicity-adjusted p-values — statin cohort 10-outcome battery";
run;
/* ── Step 4: DATA step manual Bonferroni check ──
Useful to verify PROC MULTTEST output and to document the formula in the SAP. */
data work.bonferroni_check;
set work.raw_p;
m = 10;
alpha = 0.05;
bonf_threshold = alpha / m; /* 0.05 / 10 = 0.005 */
bonf_adj_p = min(m * raw_p, 1); /* adjusted p; reject if <= 0.05 */
bonf_reject = (raw_p <= bonf_threshold); /* 1 = reject, 0 = retain */
label bonf_threshold = "Bonferroni threshold (alpha/m)"
bonf_adj_p = "Bonferroni adjusted p"
bonf_reject = "Bonferroni reject (1=yes)";
run;
proc print data=work.bonferroni_check noobs label;
var hypothesis outcome raw_p bonf_threshold bonf_adj_p bonf_reject;
run;
/* ── Step 5: Manual BH step-up in a DATA step ── */
/* Sort by raw p-value, assign ranks, compute thresholds, find largest passing rank */
proc sort data=work.raw_p out=work.raw_p_sorted; by raw_p; run;
data work.bh_check;
set work.raw_p_sorted;
rank_i = _n_; /* rank 1 = smallest p-value */
m = 10;
Q = 0.05;
bh_thresh = rank_i / m * Q; /* i/m * Q */
bh_pass = (raw_p <= bh_thresh); /* 1 = passes at this rank */
label rank_i = "BH rank"
bh_thresh = "BH threshold (i/m*Q)"
bh_pass = "BH pass (1=yes)";
run;
proc print data=work.bh_check noobs label; run;Citations
- [1]Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society, Series B. 1995;57(1):289-300.
- [2]Bender R, Lange S. Adjusting for multiple testing — when and how? Journal of Clinical Epidemiology. 2001;54(4):343-349.