​Comparing Multiple Testing Corrections: Bonferroni vs. Benjamini-Hochberg Procedures​

Bonferroni Correction

The Bonferroni correction adjusts the significance level (alpha) for each individual test by dividing the original alpha by the number of tests performed.

α' = α / m

where:

α' is the adjusted significance level for each individual test. α is the original significance level (e.g., 0.05). m is the total number of tests performed.

Benjamini-Hochberg Correction

The Benjamini-Hochberg correction controls the False Discovery Rate (FDR) by adjusting the p-values of each individual test.

pᵢ ≤ (i/m)α

where:

pᵢ is the p-value for the iᵗʰ test (sorted in ascending order). i is the rank of the iᵗʰ test. m is the total number of tests performed. α is the desired false discovery rate (e.g., 0.05). The Benjamini-Hochberg procedure rejects the null hypothesis for all tests with a p-value less than or equal to the largest k such that pₖ ≤ (k/m)α.

In simpler terms:

Bonferroni Correction: Divide the significance level by the number of tests. Benjamini-Hochberg Correction: Compare each p-value to a scaled significance level based on its rank.

In [9]:

# an example to demonstrate the Bonferroni correction and Benjamini-Hochberg correction method

import numpy as np
from scipy import stats

# Example p-values (replace with your actual p-values)
p_values = np.array([0.01, 0.03, 0.005, 0.08, 0.02, 0.05, 0.12, 0.06, 0.04, 0.09])
alpha = 0.05  # Significance level

# Bonferroni correction
bonferroni_threshold = alpha / len(p_values)
bonferroni_rejected = p_values < bonferroni_threshold
print("Bonferroni-corrected p-values:", p_values)
print("Bonferroni threshold:", bonferroni_threshold)
print("Bonferroni rejected hypotheses:", bonferroni_rejected)


# Benjamini-Hochberg correction
sorted_p_values = np.sort(p_values)
ranks = np.arange(1, len(p_values) + 1)
bh_threshold = (ranks / len(p_values)) * alpha
print(sorted_p_values)
print(bh_threshold)
print(sorted_p_values <= bh_threshold)
bh_rejected = sorted_p_values <= bh_threshold
# Find the largest k such that p_k <= k/m * alpha
k = np.max(np.where(sorted_p_values <= bh_threshold))
bh_rejected = sorted_p_values <= bh_threshold[k]
print(k)
bh_rejected_indices = np.argsort(p_values)[np.where(bh_rejected)]
bh_rejected_original = np.zeros_like(p_values, dtype=bool)
bh_rejected_original[bh_rejected_indices] = True

print("\nBenjamini-Hochberg-corrected p-values:", p_values)
print("Benjamini-Hochberg rejected hypotheses:", bh_rejected_original)
print("Benjamini-Hochberg threshold:", bh_threshold[k])

Bonferroni-corrected p-values: [0.01  0.03  0.005 0.08  0.02  0.05  0.12  0.06  0.04  0.09 ]
Bonferroni threshold: 0.005
Bonferroni rejected hypotheses: [False False False False False False False False False False]
[0.005 0.01  0.02  0.03  0.04  0.05  0.06  0.08  0.09  0.12 ]
[0.005 0.01  0.015 0.02  0.025 0.03  0.035 0.04  0.045 0.05 ]
[ True  True False False False False False False False False]
1

Benjamini-Hochberg-corrected p-values: [0.01  0.03  0.005 0.08  0.02  0.05  0.12  0.06  0.04  0.09 ]
Benjamini-Hochberg rejected hypotheses: [ True False  True False False False False False False False]
Benjamini-Hochberg threshold: 0.010000000000000002