False Discovery Rate — Advanced Methods

When performing thousands of simultaneous hypothesis tests (e.g., genome-wide association studies, neuroimaging, proteomics), controlling the family-wise error rate (FWER) is overly conservative. The false discovery rate (FDR) controls the expected proportion of false discoveries among all rejections, providing a more powerful framework for discovery-oriented science.

FDR Definition

Benjamini-Hochberg (BH) Procedure

The BH procedure is equivalent to: reject when the adjusted p-value .

Benjamini-Yekutieli (BY) Procedure

The BY procedure controls FDR under arbitrary dependence but is conservative. The penalty factor can be severe for large .

Storey's q-value

The q-value procedure then uses the adjusted threshold instead of .

Local False Discovery Rate

The local FDR provides a continuous measure of evidence against each null, unlike the BH procedure which uses ordered p-values. Efron's empirical Bayes approach estimates from the data (via a histogram or density estimate) and from the central peak of the distribution.

Empirical Bayes FDR (EBF)

Adaptive Procedures

This adapts to the fraction of non-nulls, achieving FDR rather than when many hypotheses are non-null. The adaptive procedure is more powerful than standard BH when .

Power Analysis for FDR

Dependence Structures

Python Implementation

import numpy as np
import pandas as pd
from scipy import stats
from statsmodels.stats.multitest import multipletests
import matplotlib.pyplot as plt

np.random.seed(42)

# --- Generate multiple testing scenario ---
m = 5000  # total tests
m1 = 200  # non-null hypotheses
m0 = m - m1

# Null p-values: Uniform(0,1)
p0 = np.random.uniform(0, 1, m0)
# Non-null p-values: Beta (skewed toward 0)
p1 = np.random.beta(0.5, 1, m1)
p_values = np.concatenate([p0, p1])
labels = np.concatenate([np.zeros(m0), np.ones(m1)])

# Shuffle
idx = np.random.permutation(m)
p_values = p_values[idx]
labels = labels[idx]

# --- BH Procedure ---
alpha = 0.05
reject_bh, pvals_bh, _, _ = multipletests(p_values, alpha=alpha, method='fdr_bh')
tp_bh = np.sum(reject_bh & (labels == 1))
fp_bh = np.sum(reject_bh & (labels == 0))
fn_bh = np.sum(~reject_bh & (labels == 1))
fdp_bh = fp_bh / max(np.sum(reject_bh), 1)

print("=== Benjamini-Hochberg ===")
print(f"Rejections: {reject_bh.sum()}")
print(f"True positives: {tp_bh}, False positives: {fp_bh}")
print(f"FDP: {fdp_bh:.4f} (target: {alpha})")
print(f"Power: {tp_bh / m1:.4f}")

# --- BY Procedure ---
reject_by, pvals_by, _, _ = multipletests(p_values, alpha=alpha, method='fdr_by')
tp_by = np.sum(reject_by & (labels == 1))
fp_by = np.sum(reject_by & (labels == 0))
fdp_by = fp_by / max(np.sum(reject_by), 1)

print(f"\n=== Benjamini-Yekutieli ===")
print(f"Rejections: {reject_by.sum()}")
print(f"True positives: {tp_by}, False positives: {fp_by}")
print(f"FDP: {fdp_by:.4f}")
print(f"Power: {tp_by / m1:.4f}")

# --- Storey's q-value ---
# Estimate pi0
lambdas = np.linspace(0.05, 0.95, 19)
pi0_estimates = []
for lam in lambdas:
    pi0_est = np.sum(p_values > lam) / (m * (1 - lam))
    pi0_estimates.append(pi0_est)

# Smooth estimate via cubic spline
from scipy.interpolate import UnivariateSpline
spline = UnivariateSpline(lambdas, pi0_estimates, s=0.1)
pi0_hat = min(spline(1.0), 1.0)
pi0_hat = max(pi0_hat, 0.0)  # ensure non-negative

print(f"\n=== Storey's q-value ===")
print(f"Estimated π₀: {pi0_hat:.4f} (true: {m0/m:.4f})")

# BH with pi0 adjustment
reject_storey, pvals_storey, _, _ = multipletests(p_values, alpha=alpha, method='fdr_bh')
# Adjusted p-values with pi0
pvals_storey_adj = pvals_bh / pi0_hat
reject_storey = pvals_storey_adj <= alpha

tp_storey = np.sum(reject_storey & (labels == 1))
fp_storey = np.sum(reject_storey & (labels == 0))
fdp_storey = fp_storey / max(np.sum(reject_storey), 1)

print(f"Rejections: {reject_storey.sum()}")
print(f"True positives: {tp_storey}, False positives: {fp_storey}")
print(f"FDP: {fdp_storey:.4f}")
print(f"Power: {tp_storey / m1:.4f}")

# --- Local FDR ---
# Estimate marginal density via KDE
from scipy.stats import gaussian_kde
kde = gaussian_kde(p_values, bw_method=0.1)
p_grid = np.linspace(0, 1, 500)
f_hat = kde(p_grid)

# Null density: Uniform(0,1)
f0 = np.ones_like(p_grid)

# Local FDR
fdr_local = pi0_hat * f0 / np.maximum(f_hat, 1e-10)
fdr_local = np.minimum(fdr_local, 1.0)

# Assign local FDR to each p-value
from scipy.interpolate import interp1d
fdr_interp = interp1d(p_grid, fdr_local, kind='linear', fill_value='extrapolate')
local_fdr_values = fdr_interp(p_values)
local_fdr_values = np.clip(local_fdr_values, 0, 1)

reject_lfdr = local_fdr_values <= alpha
tp_lfdr = np.sum(reject_lfdr & (labels == 1))
fp_lfdr = np.sum(reject_lfdr & (labels == 0))
fdp_lfdr = fp_lfdr / max(np.sum(reject_lfdr), 1)

print(f"\n=== Local FDR ===")
print(f"Rejections: {reject_lfdr.sum()}")
print(f"True positives: {tp_lfdr}, False positives: {fp_lfdr}")
print(f"FDP: {fdp_lfdr:.4f}")
print(f"Power: {tp_lfdr / m1:.4f}")

# --- Visualization ---
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# P-value histogram
axes[0, 0].hist(p_values[labels == 0], bins=50, alpha=0.6, label='Null', density=True)
axes[0, 0].hist(p_values[labels == 1], bins=50, alpha=0.6, label='Non-null', density=True)
axes[0, 0].set_xlabel('p-value')
axes[0, 0].set_ylabel('Density')
axes[0, 0].set_title('P-value Distribution')
axes[0, 0].legend()

# BH threshold
sorted_p = np.sort(p_values)
bh_line = np.arange(1, m + 1) / m * alpha
axes[0, 1].plot(sorted_p, bh_line, 'r-', lw=2, label='BH threshold')
axes[0, 1].plot(sorted_p, sorted_p, 'k--', lw=1, label='y = x')
axes[0, 1].set_xlabel('Sorted p-values')
axes[0, 1].set_ylabel('Threshold')
axes[0, 1].set_title('Benjamini-Hochberg Procedure')
axes[0, 1].legend()
axes[0, 1].set_xlim([0, 0.2])

# Local FDR curve
axes[1, 0].plot(p_grid, fdr_local, 'b-', lw=2, label='Local FDR')
axes[1, 0].axhline(y=alpha, color='red', ls='--', label=f'α = {alpha}')
axes[1, 0].set_xlabel('p-value')
axes[1, 0].set_ylabel('Local FDR')
axes[1, 0].set_title('Local False Discovery Rate')
axes[1, 0].legend()
axes[1, 0].set_ylim([0, 1.05])

# Power comparison
methods = ['BH', 'BY', 'Storey', 'Local FDR']
powers = [tp_bh/m1, tp_by/m1, tp_storey/m1, tp_lfdr/m1]
fdps = [fdp_bh, fdp_by, fdp_storey, fdp_lfdr]

x_pos = np.arange(len(methods))
width = 0.35
axes[1, 1].bar(x_pos - width/2, powers, width, label='Power', color='steelblue')
axes[1, 1].bar(x_pos + width/2, fdps, width, label='FDP', color='coral')
axes[1, 1].axhline(y=alpha, color='red', ls='--', lw=1, label=f'Target α = {alpha}')
axes[1, 1].set_xticks(x_pos)
axes[1, 1].set_xticklabels(methods)
axes[1, 1].set_ylabel('Rate')
axes[1, 1].set_title('Power and FDP Comparison')
axes[1, 1].legend()

plt.tight_layout()
plt.savefig("fdr_advanced.png", dpi=150, bbox_inches="tight")
plt.show()

Summary