Multiple testing

Version Compatibility

Reference examples tested with: qvalue 2.34+, IHW 1.30+, R stats (base) p.adjust, statsmodels 0.14+, scipy 1.12+.

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• Python: pip show <package> then help(module.function) to check signatures

If code throws an error, introspect the installed package and adapt to the actual API. Note: statsmodels.stats.multitest.multipletests defaults to method='hs' (Holm-Sidak, an FWER method), NOT Benjamini-Hochberg — always pass method='fdr_bh'/'fdr_by'/'bonferroni'/'holm' explicitly.

Multiple Testing Correction

"Correct p-values for testing thousands of features" -> Choose an error rate appropriate to the regime (FDR for discovery, FWER for confirmatory), apply a procedure whose dependence assumptions match the data, and report the adjusted quantity with its interpretation.

• R: p.adjust(p, method = 'BH'), qvalue::qvalue(), IHW::ihw()
• Python: statsmodels.stats.multitest.multipletests(p, method='fdr_bh')

The Single Most Important Modern Insight -- FDR vs FWER Is a Choice About Which Error Matters

The choice between false-discovery-rate and family-wise-error control is not a technicality; it is a statement about which kind of mistake is costly. In discovery (20,000 genes, thousands of peaks), tolerating a small, controlled fraction of false positives among the rejections buys enormous power — FDR is the right currency, and Bonferroni would discard nearly every true effect. In confirmatory work (a handful of pre-specified endpoints), a single false positive is unacceptable and FWER/closed testing is the standard (that regime lives in clinical-biostatistics/multiplicity-graphical). Two further levers buy back power that plain BH leaves on the table: estimating pi0 (the proportion of true nulls) turns BH into the more powerful q-value (Storey 2002 J R Stat Soc B 64:479; Storey & Tibshirani 2003 PNAS 100:9440), and weighting hypotheses by an independent informative covariate recovers power via IHW (Ignatiadis 2016 Nat Methods 13:577). The dependence structure matters: BH controls FDR under independence or positive regression dependence (PRDS); under arbitrary or negative dependence use BY (Benjamini & Yekutieli 2001 Ann Stat 29:1165).

Algorithmic Taxonomy

Method	Controls	Dependence assumption	When to use	Tool
Bonferroni	FWER	any	tiny families; confirmatory	p.adjust(method='bonferroni')
Holm	FWER	any	uniformly beats Bonferroni	p.adjust(method='holm')
Hochberg / Hommel	FWER	positive dependence	step-up FWER, more power	p.adjust(method='hochberg'/'hommel')
Benjamini-Hochberg	FDR	independence / PRDS	genome-wide discovery default	p.adjust(method='BH')
Benjamini-Yekutieli	FDR	arbitrary (incl. negative)	unknown/negative dependence	p.adjust(method='BY')
Storey q-value	pFDR	independence / weak dependence	many true positives (pi0 << 1)	qvalue::qvalue
Local FDR	posterior null prob	two-groups model	per-feature null probability	qvalue ($lfdr); locfdr
IHW	FDR	covariate independent of null p	informative covariate available	IHW::ihw
IDR	reproducibility	replicate ranks	thresholding by replicate consistency	idr (ENCODE)

Decision Tree by Scenario

Scenario	Recommended	Why
Genome-wide DE / peaks, discovery	BH or q-value at FDR 0.05	controlled false-positive fraction; high power
Many true positives expected	q-value (estimates pi0)	more powerful than BH when pi0 << 1
Strong/unknown/negative dependence	BY	BH guarantee needs PRDS
Informative covariate (mean expr, peak width)	IHW	data-driven weights recover power
Per-feature "is this one real?"	local FDR	posterior null probability, not tail average
Reporting CIs only on significant hits	FCR-adjusted intervals	naive selected CIs under-cover
Small confirmatory gene panel	Bonferroni/Holm	FWER appropriate; power loss acceptable
GWAS	genome-wide threshold ~5e-8	~1M effective independent tests
Confirmatory trial, few endpoints	-> clinical-biostatistics/multiplicity-graphical	closed testing / gatekeeping
Applying padj to a finished DE table	-> differential-expression/de-results	method choice here; application there

FDR -- Benjamini-Hochberg and the q-value

# Benjamini-Hochberg adjusted p-values (the genome-wide default)
padj <- p.adjust(pvalues, method = 'BH')
sum(padj < 0.05)                                  # discoveries at FDR 5%

# Storey q-value: estimates pi0 (fraction of true nulls) for more power when pi0 << 1
library(qvalue)
qobj <- qvalue(pvalues)
qobj$pi0                                           # estimated proportion of true nulls
q   <- qobj$qvalues                                # min FDR at which each feature is called
lfdr <- qobj$lfdr                                  # local FDR: posterior P(null | statistic)

Dependence -- When BH Is Not Enough (BY)

# BH controls FDR under independence or positive regression dependence (PRDS).
# Under arbitrary or negative dependence, use Benjamini-Yekutieli (more conservative).
padj_by <- p.adjust(pvalues, method = 'BY')        # valid under any dependence structure

Covariate-Weighted FDR -- IHW

# Weight hypotheses by an INDEPENDENT informative covariate (e.g. mean expression),
# which must be independent of the p-value under the null. Recovers power vs plain BH.
library(IHW)
res <- ihw(pvalue ~ mean_expression, data = de_table, alpha = 0.05)
de_table$padj_ihw <- adj_pvalues(res)
rejections(res)

Independent Filtering -- Power for Free, If the Filter Is Independent

Filtering out features before testing increases power only if the filter statistic is independent of the test statistic under the null (Bourgon, Gentleman & Huber 2010 PNAS 107:9546). Overall mean count is independent and is why DESeq2 filters low-count genes automatically; a pre-test on variance or a preliminary t-test is not independent and biases the FDR. The DE filtering itself is executed in differential-expression; this skill governs whether a proposed filter is legitimate.

Python Equivalent (mind the default)

from statsmodels.stats.multitest import multipletests
# DEFAULT method is 'hs' (Holm-Sidak, FWER) -- ALWAYS pass method explicitly.
rej, padj, _, _ = multipletests(pvalues, alpha=0.05, method='fdr_bh')   # Benjamini-Hochberg
rej_by, padj_by, _, _ = multipletests(pvalues, alpha=0.05, method='fdr_by')  # BY

GWAS and the Family-Definition Problem

The genome-wide significance threshold of ~5e-8 is a Bonferroni-style bound for roughly one million effectively independent common-variant tests; Dudbridge & Gusnanto 2008 (Genet Epidemiol 32:227) derived ~7.2e-8 for European-ancestry data, near the standard 5e-8. The GWAS test machinery lives in population-genetics/association-testing. More broadly, what counts as "the family" of tests is an analyst decision and part of the garden of forking paths: correcting within one contrast, across all contrasts, or across a whole paper are different alpha budgets. Pre-specify the family before seeing results.

Reconciliation: When Methods Disagree

Pattern	Likely cause	Action
q-value finds many more hits than BH	pi0 << 1 (many true positives)	q-value legitimately more powerful; report pi0
BY far more conservative than BH	strong/negative dependence penalty	if dependence is positive, BH is justified; state the assumption
IHW and BH differ substantially	informative, null-independent covariate	IHW gain is real if independence holds; verify the covariate
Filtering changed the hit count	filter not independent of the test statistic	use a null-independent filter (mean count), not variance/preliminary test
Per-feature local FDR high but BH q low	tail-average vs per-feature interpretation	report both; local FDR answers "is THIS one real?"

Per-Method Failure Modes

Bonferroni on a transcriptome

• Trigger: Bonferroni across 20,000 genes in a discovery study.
• Mechanism: FWER control is far too strict for discovery.
• Symptom: almost nothing significant; true effects discarded.
• Fix: BH or q-value at a target FDR.

BH under arbitrary/negative dependence

• Trigger: BH on strongly/negatively correlated statistics.
• Mechanism: BH guarantee requires independence or PRDS (Benjamini-Yekutieli 2001).
• Symptom: realized FDR exceeds nominal.
• Fix: BY when dependence is unknown or negative.

statsmodels default is not BH

• Trigger: multipletests(p) expecting Benjamini-Hochberg.
• Mechanism: default method='hs' (Holm-Sidak, FWER).
• Symptom: far fewer significant calls than expected.
• Fix: pass method='fdr_bh' explicitly.

Non-independent filtering

• Trigger: filter on variance or a preliminary test before the main test.
• Mechanism: filter statistic correlated with the test statistic under the null (Bourgon 2010).
• Symptom: anti-conservative FDR.
• Fix: filter only on a null-independent statistic (overall mean count).

Selected CIs without FCR adjustment

• Trigger: reporting unadjusted CIs only for significant features.
• Mechanism: selection induces under-coverage (false coverage rate).
• Symptom: intervals too narrow; replication misses.
• Fix: FCR-adjusted intervals for the selected set.

Quantitative Thresholds

Threshold	Source	Rationale
FDR < 0.05 discovery default	Benjamini-Hochberg 1995 JRSS-B 57:289	5% of calls expected false
FDR < 0.10 exploratory	common practice	more leads at higher false fraction
q-value uses estimated pi0	Storey 2002 JRSS-B 64:479	power gain when pi0 << 1
BH valid under independence/PRDS; else BY	Benjamini-Yekutieli 2001 Ann Stat 29:1165	dependence governs validity
GWAS ~5e-8 (7.2e-8 derived)	Dudbridge-Gusnanto 2008 Genet Epidemiol 32:227	~1M effective tests
Filter must be null-independent	Bourgon 2010 PNAS 107:9546	otherwise FDR is biased

Common Errors

Error / symptom	Cause	Solution
Almost nothing significant genome-wide	Bonferroni in a discovery study	BH or q-value
Realized FDR exceeds nominal	BH under negative dependence	BY
Far fewer hits than expected in Python	statsmodels default 'hs'	method='fdr_bh'
FDR biased after pre-filtering	non-independent filter	filter on mean count only
Replication misses "significant" effects	unadjusted selected CIs	FCR-adjusted intervals

References

• Benjamini Y, Hochberg Y. 1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B 57:289-300.
• Benjamini Y, Yekutieli D. 2001. The control of the false discovery rate in multiple testing under dependency. Ann Stat 29:1165-1188.
• Storey JD. 2002. A direct approach to false discovery rates. J R Stat Soc B 64:479-498.
• Storey JD, Tibshirani R. 2003. Statistical significance for genomewide studies. PNAS 100:9440-9445.
• Efron B. 2008. Microarrays, empirical Bayes and the two-groups model. Stat Sci 23:1-22.
• Bourgon R, Gentleman R, Huber W. 2010. Independent filtering increases detection power for high-throughput experiments. PNAS 107:9546-9551.
• Ignatiadis N, Klaus B, Zaugg JB, Huber W. 2016. Data-driven hypothesis weighting increases detection power in genome-scale multiple testing. Nat Methods 13:577-580.
• Li Q, Brown JB, Huang H, Bickel PJ. 2011. Measuring reproducibility of high-throughput experiments. Ann Appl Stat 5:1752-1779.
• Dudbridge F, Gusnanto A. 2008. Estimation of significance thresholds for genomewide association scans. Genet Epidemiol 32:227-234.

Related Skills

• power-analysis - The FDR target feeds the power/EDR calculation
• sample-size - Replicate number depends on the FDR threshold chosen here
• batch-design - Surrogate variables change the effective number of tests
• differential-expression/de-results - Where the padj column is applied to a DE table
• population-genetics/association-testing - GWAS genome-wide significance machinery
• pathway-analysis/go-enrichment - Correcting enrichment p-values
• clinical-biostatistics/multiplicity-graphical - Confirmatory FWER / closed testing for trials with few endpoints