GPTomics/bio-experimental-design-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