Multiple Comparisons › Practical Implementation with scipy.stats
Michael BrenndoerferJanuary 9, 20263 min read
Machine Learning from ScratchAll chapters Part 2 of 9Hypothesis Testing
Multiple Comparisons › Practical Implementation with scipy.statsOne section · Sections in this chapter7 / 12 ↓
Sections in this chapter12
Reading tools Definitions & inline links
Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.
Article links
Make inline references clickable
Let's work through a complete example using scipy's built-in multiple comparison functions.
In[15]:
Codefrom scipy.stats import false_discovery_control
# Simulate a genomics study: testing 1000 genes
np.random.seed(42)
n_genes = 1000
n_real_effects = 50
# Generate p-values
# Most genes: null hypothesis true (uniform p-values)
p_null = np.random.uniform(0, 1, n_genes - n_real_effects)
# Some genes: real effects (small p-values)
p_alt = np.random.beta(0.5, 5, n_real_effects) # Skewed toward 0
p_values = np.concatenate([p_null, p_alt])
is_null = np.array(
[True] * (n_genes - n_real_effects) + [False] * n_real_effects
)
# Shuffle to mix them up
shuffle_idx = np.random.permutation(n_genes)
p_values = p_values[shuffle_idx]
is_null = is_null[shuffle_idx]
print(f"Testing {n_genes} genes ({n_real_effects} with true effects)")
print(f"Without correction: {np.sum(p_values < 0.05)} significant (α = 0.05)")
print()
# Apply different corrections
# Bonferroni
bonf_sig = p_values < (0.05 / n_genes)
print(f"Bonferroni (FWER): {np.sum(bonf_sig)} significant")
print(f" True positives: {np.sum(bonf_sig & ~is_null)}")
print(f" False positives: {np.sum(bonf_sig & is_null)}")
print()
# Benjamini-Hochberg using scipy
bh_adjusted = false_discovery_control(p_values, method="bh")
bh_sig = bh_adjusted < 0.05
print(f"Benjamini-Hochberg (FDR): {np.sum(bh_sig)} significant")
print(f" True positives: {np.sum(bh_sig & ~is_null)}")
print(f" False positives: {np.sum(bh_sig & is_null)}")
print(
f" Empirical FDR: {np.sum(bh_sig & is_null) / max(np.sum(bh_sig), 1):.2%}"
)
print()
# Benjamini-Yekutieli using scipy
by_adjusted = false_discovery_control(p_values, method="by")
by_sig = by_adjusted < 0.05
print(f"Benjamini-Yekutieli (FDR, conservative): {np.sum(by_sig)} significant")
print(f" True positives: {np.sum(by_sig & ~is_null)}")
print(f" False positives: {np.sum(by_sig & is_null)}")Out[15]:
ConsoleTesting 1000 genes (50 with true effects) Without correction: 72 significant (α = 0.05) Bonferroni (FWER): 0 significant True positives: 0 False positives: 0 Benjamini-Hochberg (FDR): 0 significant True positives: 0 False positives: 0 Empirical FDR: 0.00% Benjamini-Yekutieli (FDR, conservative): 0 significant True positives: 0 False positives: 0
Comparing Methods on Real-World Scale
Out[16]:
Visualization