Michael Brenndoerfer | Data & AI, Private Equity, Technology

Michael Brenndoerfer

Multiple Comparisons › FDR-Controlling Methods

Michael BrenndoerferJanuary 9, 20265 min read

Machine Learning from ScratchAll chapters

Part 2 of 9Hypothesis Testing

Chapter 21Multiple Comparisons › FDR-Controlling Methods

Reading section 5 of 12 · FDR-Controlling Methods · 5 minSection 5 of 12 · 5 min

Section 5 of 12 · 5 min read

Multiple Comparisons › FDR-Controlling MethodsOne section · Sections in this chapter5 / 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

Benjamini-Hochberg ProcedureLink Copied

The Benjamini-Hochberg (BH) procedure (1995) is the most widely used method for controlling the false discovery rate. It's less conservative than FWER methods, allowing more discoveries while controlling the expected proportion of false positives among rejected hypotheses.

The algorithm:

Order p-values: $p_{(1)} \leq p_{(2)} \leq \ldots \leq p_{(m)}$
Find the largest $k$ such that $p_{(k)} \leq \frac{k}{m} \cdot q$ (where $q$ is the desired FDR level)
Reject all hypotheses with p-values $\leq p_{(k)}$

Mathematical justification: The BH procedure controls FDR at level $q$ under the assumption that the test statistics are independent (or positively dependent under some conditions). The key insight is that the threshold $\frac{k}{m} \cdot q$ increases with rank, allowing more rejections among the smaller p-values.

Adjusted p-values: The BH-adjusted p-value for the hypothesis with rank $k$ is:

\tilde{p}_{(k)} = \min_{j \geq k} \left\{ \min\left(\frac{m \cdot p_{(j)}}{j}, 1\right) \right\}

This ensures monotonicity: if $p_{(i)} < p_{(j)}$ , then $\tilde{p}_{(i)} \leq \tilde{p}_{(j)}$ .

In[10]:

Code

def benjamini_hochberg(p_values, q=0.05):
    """
    Apply Benjamini-Hochberg procedure for FDR control.

    Parameters
    ----------
    p_values : array-like
        Raw p-values from multiple tests
    q : float
        Desired FDR level

    Returns
    -------
    adjusted_p : ndarray
        BH-adjusted p-values
    rejected : ndarray
        Boolean array indicating which hypotheses to reject
    """
    p = np.asarray(p_values)
    m = len(p)

    # Sort p-values
    sorted_idx = np.argsort(p)
    p_sorted = p[sorted_idx]

    # BH thresholds
    thresholds = q * np.arange(1, m + 1) / m

    # Find largest k where p(k) <= threshold(k)
    significant = p_sorted <= thresholds
    if np.any(significant):
        k_max = np.max(np.where(significant)[0])
        rejected_sorted = np.arange(m) <= k_max
    else:
        rejected_sorted = np.zeros(m, dtype=bool)

    # Map back to original order
    rejected = np.zeros(m, dtype=bool)
    rejected[sorted_idx] = rejected_sorted

    # Calculate adjusted p-values (ensuring monotonicity)
    adjusted_sorted = np.zeros(m)
    adjusted_sorted[m - 1] = p_sorted[m - 1]
    for i in range(m - 2, -1, -1):
        adjusted_sorted[i] = min(
            adjusted_sorted[i + 1], m * p_sorted[i] / (i + 1)
        )
    adjusted_sorted = np.minimum(adjusted_sorted, 1.0)

    adjusted_p = np.zeros(m)
    adjusted_p[sorted_idx] = adjusted_sorted

    return adjusted_p, rejected

Benjamini-Yekutieli ProcedureLink Copied

The Benjamini-Yekutieli (BY) procedure (2001) extends BH to handle arbitrary dependence structures among tests. It's more conservative than BH but valid under any dependence.

The method: Replace the BH threshold $\frac{k}{m} \cdot q$ with $\frac{k}{m \cdot c(m)} \cdot q$ , where:

c(m) = \sum_{i=1}^{m} \frac{1}{i} \approx \ln(m) + 0.577

This correction factor accounts for potential negative correlations among tests.

In[11]:

Code

def benjamini_yekutieli(p_values, q=0.05):
    """
    Apply Benjamini-Yekutieli procedure for FDR control under arbitrary dependence.
    """
    p = np.asarray(p_values)
    m = len(p)

    # Correction factor for arbitrary dependence
    c_m = np.sum(1 / np.arange(1, m + 1))

    # Sort p-values
    sorted_idx = np.argsort(p)
    p_sorted = p[sorted_idx]

    # BY thresholds
    thresholds = q * np.arange(1, m + 1) / (m * c_m)

    # Find largest k where p(k) <= threshold(k)
    significant = p_sorted <= thresholds
    if np.any(significant):
        k_max = np.max(np.where(significant)[0])
        rejected_sorted = np.arange(m) <= k_max
    else:
        rejected_sorted = np.zeros(m, dtype=bool)

    rejected = np.zeros(m, dtype=bool)
    rejected[sorted_idx] = rejected_sorted

    return rejected

Visualizing the BH ProcedureLink Copied

Out[12]:

Visualization