False Discovery Rate — Advanced Methods

Advanced Statistical Methods

Controlling False Positives in Massive Testing

FDR methods control the expected proportion of false discoveries among rejected hypotheses, enabling powerful inference in genome-wide studies, neuroimaging, and other high-throughput applications. Benjamini-Hochberg and Storey's q-value are foundational procedures.

Genomics — Identify differentially expressed genes across thousands of tests while controlling false positives
Neuroimaging — Detect active brain regions from millions of voxels in fMRI data
Astronomy — Classify celestial objects from massive survey data with controlled error rates

FDR methods let you reject thousands of hypotheses while keeping false alarms in check.

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

DfFalse Discovery Rate

Let $V$ be the number of false rejections and $R$ the total number of rejections. The FDR is defined as:

\text{FDR} = \mathbb{E}\left[\frac{V}{R \vee 1}\right]

where $R \vee 1 = \max(R, 1)$ handles the case of no rejections. The positive FDR (pFDR) conditions on $R > 0$ :

\text{pFDR} = \mathbb{E}\left[\frac{V}{R} \,\bigg|\, R > 0\right]

FDR vs FWER

FWER = $\mathbb{P}(V \geq 1)$ : Probability of making at least one false discovery. Controls even one false positive.
FDR = $\mathbb{E}[V/R]$ : Expected proportion of false discoveries. Allows a controlled fraction of errors.
FWER control is appropriate when any false positive is costly (e.g., drug safety). FDR is appropriate when the goal is discovery (e.g., gene selection for follow-up).

Benjamini-Hochberg (BH) Procedure

Benjamini-Hochberg Procedure

Given $m$ p-values $p_{(1)} \leq p_{(2)} \leq \dots \leq p_{(m)}$ (sorted), for a target FDR level $q$ :

Find the largest $k$ such that $p_{(k)} \leq \frac{k}{m}q$
Reject all hypotheses $H_{(1)}, \dots, H_{(k)}$
If no such $k$ exists, reject none

ThBH Control Under Independence

Under independence of all $m$ p-values, the BH procedure controls FDR at level:

\text{FDR} \leq \frac{m_0}{m}q \leq q

where $m_0$ is the number of true null hypotheses. The bound is tight when all nulls are true ( $m_0 = m$ ).

The BH procedure is equivalent to: reject when the adjusted p-value $p_{(i)}^{\text{BH}} = \min\left(\frac{m}{i}p_{(i)}, 1\right) \leq q$ .

Benjamini-Yekutieli (BY) Procedure

BY Adjustment

Under arbitrary dependence among p-values, the BH procedure with threshold $\frac{k}{m}q$ is replaced by:

p_{(k)} \leq \frac{k}{m \cdot c(m)} \cdot q

where $c(m) = \sum_{i=1}^{m} \frac{1}{i} \approx \log(m) + \gamma$ (Euler-Mascheroni constant $\gamma \approx 0.5772$ ).

The BY procedure controls FDR under arbitrary dependence but is conservative. The penalty factor $c(m) \approx \log(m)$ can be severe for large $m$ .

Dependence Structure Matters

The BH procedure's conservativeness under dependence varies:

Positive regression dependence on subsets (PRDS): BH still controls FDR exactly (Benjamini & Yekutieli, 2001)
AR(1) correlation among test statistics: BH controls FDR at level $q$ with small inflation
Arbitrary dependence: Only BY guarantees control, but with the $\log(m)$ penalty

In practice, block-diagonal or weak dependence structures allow BH with minimal inflation.

Storey's q-value

Dfq-value

The q-value of a test statistic $t$ is the minimum FDR at which the test would be called significant:

q(t) = \min_{t' \geq t} \text{pFDR}(t')

Estimating π₀

Storey (2002) estimates the proportion of true nulls via:

\hat{\pi}_0 = \frac{\#\{p_i > \lambda\}}{m(1 - \lambda)}

for tuning parameter $\lambda \in (0, 1)$ (typically $\lambda = 0.5$ ). A smoother estimate uses cubic spline regression:

\hat{\pi}_0(\lambda) = \text{cubic spline}\left(\{(\lambda_j, \hat{\pi}_0(\lambda_j))\}_{j=1}^{L}\right)

The q-value procedure then uses the adjusted threshold $\frac{p_{(i)}}{i/m} \cdot \hat{\pi}_0$ instead of $\frac{p_{(i)}}{i/m}$ .

Thq-value Optimality

The q-value procedure is optimal in the sense that it rejects exactly those hypotheses that would be rejected by the BH procedure at some FDR level, while providing an estimate of the minimum FDR for each rejection.

Local False Discovery Rate

DfLocal FDR (Efron, 2001)

The local FDR is the posterior probability that a hypothesis is true given its test statistic $z$ :

\text{fdr}(z) = P(H_0 \mid Z = z) = \frac{\pi_0 f_0(z)}{f(z)}

where $f_0$ is the null density, $f$ is the marginal density, and $\pi_0 = P(H_0)$ .

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 $f(z)$ from the data (via a histogram or density estimate) and $f_0(z)$ from the central peak of the distribution.

Local FDR vs BH

Local FDR provides individual-level evidence for each test (continuous, not step function)
BH provides a binary decision at a given FDR level
Local FDR can be more powerful when the null distribution is accurately estimated
Local FDR is sensitive to misspecification of $f_0$ ; BH is more robust

Empirical Bayes FDR (EBF)

Empirical Bayes Framework

Model the mixture:

f(z) = \pi_0 f_0(z) + (1 - \pi_0) f_1(z)

The posterior probability of the null is:

\text{EBF}(z) = \frac{\pi_0 f_0(z)}{f(z)}

In practice, $f_0$ is estimated from the p-values near 1 (conservative approach), or from a known null distribution. $\pi_0$ is estimated as described above.

EBF Estimation Steps

From observed p-values $p_1, \dots, p_m$ , estimate $\hat{\pi}_0$ using the slope of $1 - \hat{F}(p)$ near $p = 1$ .
Under the null, $f_0(p) = 1$ (uniform on $[0,1]$ ).
The marginal density is $\hat{f}(p) = \hat{F}'(p)$ estimated via kernel smoothing.
Compute $\widehat{\text{EBF}}(p) = \hat{\pi}_0 / \hat{f}(p)$ .
Reject when $\widehat{\text{EBF}}(p) \leq q$ .

Adaptive Procedures

DfAdaptive FDR

The adaptive BH procedure (Storey et al., 2004) uses an estimate $\hat{\pi}_0$ to sharpen the threshold:

\text{Reject when } p_{(i)} \leq \frac{i}{m} \cdot \frac{q}{1 - (1-q)^{\hat{\pi}_0}} \approx \frac{i}{m} \cdot \frac{q}{\hat{\pi}_0}

This adapts to the fraction of non-nulls, achieving FDR $\approx \hat{\pi}_0 q$ rather than $q$ when many hypotheses are non-null. The adaptive procedure is more powerful than standard BH when $\pi_0 < 1$ .

ThAdaptive BH Power Gain

For sparse signals ( $\pi_0$ close to 1), the adaptive and standard BH procedures have similar power. For dense signals ( $\pi_0 \ll 1$ ), the adaptive procedure can have substantially higher power, approaching the oracle procedure that knows $\pi_0$ .

Power Analysis for FDR

FDR Power

The power of an FDR procedure at level $q$ is:

\text{Power}(q) = \mathbb{E}\left[\frac{R - V}{m_1}\right]

where $m_1 = m - m_0$ is the number of non-null hypotheses. For the BH procedure with $m$ tests:

\text{Power}_{\text{BH}}(q) \approx 1 - \Phi\left(z_{1-q\pi_0/m_1} - \sqrt{n}\,\Delta\right)

where $\Delta$ is the non-centrality parameter (effect size × $\sqrt{n}$ ) and $z_\alpha$ is the standard normal quantile.

Factors Affecting FDR Power

Effect size: Larger effects are detected with higher power
Sample size: Power increases with $\sqrt{n}$
Proportion of non-nulls ( $\pi_1 = 1 - \pi_0$ ): More non-nulls → more discoveries
FDR level ( $q$ ): Higher $q$ → more discoveries but more false positives
Dependence structure: Positive dependence can increase or decrease power depending on the configuration
p-value distribution under $H_1$ : Concentrated near 0 → higher power

Dependence Structures

DfPRDS (Positive Regression Dependence on Subsets)

Random variables $X_1, \dots, X_m$ satisfy PRDS if for any increasing set $A$ and any $j$ :

P(X_j \in A \mid X_i = x_i \text{ for } i \neq j) \leq P(X_j \in A \mid X_i = x_i \text{ for } i \in S)

for any subset $S \subseteq \{1, \dots, m\} \setminus \{j\}$ and all $x_i$ in the support.

ThBH Under PRDS

If the p-values satisfy PRDS (which includes the case of independent p-values and many common correlation structures), the BH procedure controls FDR at level $\frac{m_0}{m}q \leq q$ .

Anti-Conservative Dependence

The BH procedure can be anti-conservative (FDR > $q$ ) under certain negative dependence structures:

Martingale difference dependence: BH can exceed $q$ by a factor up to $\sqrt{\log m}$
Tail dependence: When extreme p-values are correlated, false discoveries cluster

For arbitrary dependence, the BY procedure or permutation-based FDR estimation is required.

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

Key Takeaways: FDR Advanced Methods

The BH procedure controls FDR at level $\frac{m_0}{m}q$ under independence (and PRDS); the BY procedure extends to arbitrary dependence at the cost of a $\log(m)$ penalty.
Storey's q-value estimates $\pi_0$ to sharpen BH thresholds, achieving higher power when many hypotheses are non-null. The adaptive BH procedure provides a practical approximation.
Local FDR provides continuous, per-test evidence via empirical Bayes; it is more powerful than BH when the null density is accurately estimated but sensitive to misspecification.
FDR power analysis depends on effect size, $\sqrt{n}$ , $\pi_1$ , the FDR level, and the p-value distribution under $H_1$ .
Dependence structure is critical: PRDS guarantees BH control; negative dependence can be anti-conservative. In practice, validate with permutation-based FDR estimation.

False Discovery Rate — Advanced Methods

False Discovery Rate — Advanced Methods

Controlling False Positives in Massive Testing

FDR Definition

DfFalse Discovery Rate

Benjamini-Hochberg (BH) Procedure

Benjamini-Hochberg Procedure

ThBH Control Under Independence

Benjamini-Yekutieli (BY) Procedure

BY Adjustment

Storey's q-value

Dfq-value

Estimating π₀

Thq-value Optimality

Local False Discovery Rate

DfLocal FDR (Efron, 2001)

Empirical Bayes FDR (EBF)

Empirical Bayes Framework

EBF Estimation Steps

Adaptive Procedures

DfAdaptive FDR

ThAdaptive BH Power Gain

Power Analysis for FDR

FDR Power

Dependence Structures

DfPRDS (Positive Regression Dependence on Subsets)

ThBH Under PRDS

Python Implementation

Summary

Key Takeaways: FDR Advanced Methods

Premium Content

Need Expert Statistics Help?