P-Value Correction for Multiple Tests

P-Value Correction for Multiple Tests#

Family-Wise Error Rate (FWER): the probability of rejecting the null hypothesis when it is actually true while performing multiple hypotheses tests.

We often need do account for multiple tests by adjusting the \(\alpha\) to keep the type 1 error rate down. Common methods to adjusting \(\alpha\) include:

Bonferroni: When strict control over false positives is paramount or when number of tests is relatively small.
Holm-Bonferroni: Balances power and strictness. Default for most business/scientific scenarios.
False Discovery Rate (FDR): For exploratory analyses or many tests. Useful when maximizing discovery outweighs risk of false positives.

How Each Method is Calculated#

For each of the examples, assume we have a family of three tests. The test is considered statistically significant when \(\text{p-value} \leq \alpha \)

Test	P-value	\( \alpha \)	Statistically Significant?
1	0.01	0.05	Yes ✅
2	0.02	0.05	Yes ✅
3	0.05	0.05	Yes ✅

Bonferroni Correction#

Divide \(\alpha\) by the number of tests (\(m\)) to get a \(\alpha_{adjusted}\) :

\[ \alpha_{adjusted} = \frac{\alpha}{m} \]

Test	P-value	Adjusted \(\alpha\)	Statistically Significant?
1	0.01	0.05/3 = 0.0167	Yes ✅
2	0.02	0.05/3 = 0.0167	No ❌
3	0.05	0.05/3 = 0.0167	No ❌

Holm-Bonferroni Correction#

Compare each p-value against decreasing alpha thresholds.

Order p-values smallest to largest: \( p_{(1)} \le p_{(2)} \le ... \le p_{(m)} \)
For the ith test, the adjusted alpha is:

\[ \alpha_{i} = \frac{\alpha}{m - i + 1} \]

Test	P-value	Adjusted \(\alpha\)	Statistically Significant?
1	0.01	0.05/3 = 0.0167	Yes ✅
2	0.02	0.05/2 = 0.025	Yes ✅
3	0.05	0.05/1 = 0.05	No ❌

False discovery rate (FDR)#

There are a few different methods to control the FDR. Here is the BH step-up procedure:

Order p-values smallest to largest: \( p_{(1)} \le p_{(2)} \le ... \le p_{(m)} \)
For each p-value, compute the threshold:

\[ \alpha_i = \frac{i}{m}\alpha \]

Test	P-value	Adjusted \(\alpha\)	Statistically Significant?
1	0.01	(1/3)*0.05 = 0.0167	Yes ✅
2	0.02	(2/3)*0.05 = 0.0333	Yes ✅
3	0.05	(3/3)*0.05 = 0.05	Yes ✅

There are several methods available on statsmodels.

Benjamini–Hochberg: Proven for independent test statistics, and has been shown to be robust to violations of the independence assumption.
Benjamini–Yekutieli: A more conservative extension that controls for correlated tests.

Matti Pirinen has great lecture notes online that dive deeper into the details of FDR.

Multiple Tests Example#

Three independent two-sample t-tests are performed comparing random data for three groups

import numpy as np
import scipy.stats as stats
from statsmodels.stats.multitest import multipletests

rng = np.random.default_rng()

alpha = 0.05
sample_size = 100
n_tests = 3

data1 = rng.standard_normal(sample_size)
data2 = rng.standard_normal(sample_size)
data3 = rng.standard_normal(sample_size)

t_stat1, p_val1 = stats.ttest_ind(data1, data2)
t_stat2, p_val2 = stats.ttest_ind(data1, data3)
t_stat3, p_val3 = stats.ttest_ind(data2, data3)

p_values = [p_val1, p_val2, p_val3]
print(f"Original p-values: {p_val1:.4f}, {p_val2:.4f}, {p_val3:.4f}")

print("\nBonferroni Correction Results:")
rejected, pvals_corrected, _, _ = multipletests(p_values, alpha=0.05, method='bonferroni')
for i, (orig_p, adj_p, rej) in enumerate(zip(p_values, pvals_corrected, rejected), start=1):
    status = 'Rejected' if rej else 'Not rejected'
    print(f"\tTest {i}: Original p = {orig_p:.4f}, Adjusted p = {adj_p:.4f}, {status}")

print("\nFDR with Benjamini-Hochberg Procedure:")
rejected, pvals_corrected, _, _ = multipletests(p_values, alpha=0.05, method='fdr_bh')
for i, (orig_p, adj_p, rej) in enumerate(zip(p_values, pvals_corrected, rejected), start=1):
    status = 'Rejected' if rej else 'Not rejected'
    print(f"\tTest {i}: Original p = {orig_p:.4f}, Adjusted p = {adj_p:.4f}, {status}")

print("\nFDR with Benjamini-Yekutieli Procedure:")
rejected, pvals_corrected, _, _ = multipletests(p_values, alpha=0.05, method='fdr_by')
for i, (orig_p, adj_p, rej) in enumerate(zip(p_values, pvals_corrected, rejected), start=1):
    status = 'Rejected' if rej else 'Not rejected'
    print(f"\tTest {i}: Original p = {orig_p:.4f}, Adjusted p = {adj_p:.4f}, {status}")

Original p-values: 0.5700, 0.6228, 0.9499

Bonferroni Correction Results:
	Test 1: Original p = 0.5700, Adjusted p = 1.0000, Not rejected
	Test 2: Original p = 0.6228, Adjusted p = 1.0000, Not rejected
	Test 3: Original p = 0.9499, Adjusted p = 1.0000, Not rejected

FDR with Benjamini-Hochberg Procedure:
	Test 1: Original p = 0.5700, Adjusted p = 0.9342, Not rejected
	Test 2: Original p = 0.6228, Adjusted p = 0.9342, Not rejected
	Test 3: Original p = 0.9499, Adjusted p = 0.9499, Not rejected

FDR with Benjamini–Yekutieli Procedure:
	Test 1: Original p = 0.5700, Adjusted p = 1.0000, Not rejected
	Test 2: Original p = 0.6228, Adjusted p = 1.0000, Not rejected
	Test 3: Original p = 0.9499, Adjusted p = 1.0000, Not rejected