跳转至

The Benjamini--Hochberg Procedure Can Fail to Control the FDR for Correlated Two-Sided Gaussian Tests

作者: Edgar Dobriban
主题: 数理统计 / 假设检验
相关性: 8/10
链接: https://arxiv.org/abs/2607.12208

一、领域脉络与小综述

这个方向是什么

这个子方向研究的是多重假设检验中错误发现率（FDR）的控制问题，特别是当检验统计量之间存在相关性时，经典的Benjamini-Hochberg（BH）程序是否仍然能保证FDR不超过名义水平α。这是一个从1995年BH程序提出以来就持续活跃的核心问题，其成熟度表现为：在独立性和特定正相关结构下已有严格证明，但在更一般的相关结构（尤其是两尾检验）下，理论结果长期滞后于实践信念。

发展脉络（history）

  1. 奠基工作：Benjamini和Hochberg (1995) 证明了在独立p值下，BH程序控制FDR ≤ π₀α（π₀是真零假设比例）。Benjamini和Yekutieli (2001) 将这一结果推广到正回归依赖性（PRDS）条件，并给出了任意依赖下的保守校正（BY程序，用调和数Hₘ放大临界值）。这两篇奠定了FDR控制的理论基础。

  2. 主要进展：Genovese和Wasserman (2002, 2004) 发展了渐近框架，将BH阈值与p值经验CDF和直线t/α的交叉点联系起来，为后续渐近分析提供了工具。Storey等人 (2004) 和Finner等人 (2007) 进一步研究了混合模型和交换相关结构下的渐近行为。Blanchard和Roquain (2008) 给出了FDR控制的“自洽性”和“依赖控制”两个充分条件，统一了多种程序的分析。

  3. 当前frontier：对于两尾高斯检验这一具体且重要的设定，情况变得复杂。Benjamini和Yekutieli (2001) 指出，非负相关对单尾检验是PRDS，但两尾检验的折叠性质（{Pᵢ ≤ t} = {Xᵢ ≥ cₜ} ∪ {Xᵢ ≤ -cₜ}）破坏了单调性。Reiner-Benaim (2007) 和Benjamini (2010) 提供了模拟和部分理论证据，支持“BH在任意相关两尾高斯检验下仍控制FDR”这一猜想。Sarkar (2023) 明确将此称为“广泛相信的猜想”，并给出了在已知协方差矩阵下的部分证明（加权检验统计量）。Fithian和Lei (2022)、Sarkar和Zhang (2025)、Ghosh和Sarkar (2025) 则转向开发依赖调整的替代方法（如dBH、shifted BH），而非证明普通BH的控制性。

  4. 本文的位置：本文直接反驳了这个持续二十年的猜想。通过构造一个特定的高斯因子模型，并利用严格的区间算术证书，证明在名义水平α=0.01下，对于足够大的假设数量，BH程序的FDR严格超过0.0104。这是一个反例构造，而非正面证明。

子线索聚类

  • 线索一：独立/正相关下的FDR控制。核心工作：Benjamini和Hochberg (1995)、Benjamini和Yekutieli (2001)、Sarkar (2002)、Blanchard和Roquain (2008)。这一簇建立了BH程序在理想依赖结构下的理论保证。
  • 线索二：任意依赖下的保守控制与渐近分析。核心工作：Benjamini和Yekutieli (2001)（BY校正）、Genovese和Wasserman (2002, 2004)、Storey等人 (2004)、Finner等人 (2007)。这一簇研究在更弱假设下的渐近行为或保守校正。
  • 线索三：依赖调整方法。核心工作：Fithian和Lei (2022)（dBH）、Sarkar (2023)、Sarkar和Zhang (2025)、Ghosh和Sarkar (2025)。这一簇放弃证明普通BH的控制性，转而开发利用相关结构信息的新程序。
  • 线索四：两尾检验的特殊性。核心工作：Reiner-Benaim (2007)、Benjamini (2010)、Sarkar (2023)、Ghosh和Sarkar (2025)。这一簇专门关注两尾高斯检验，其核心困难在于两尾折叠破坏了PRDS所需的单调性。

这个方向在追问的核心问题

  1. BH程序在相关两尾高斯检验下是否控制FDR？ 这是本文直接回答的问题，答案是否定的。
  2. 如果控制失败，FDR的偏离有多大？是否存在一个通用上界？ 本文的反例显示偏离很小（约3.59%），但未给出通用上界。
  3. 对于有限样本，BH程序的行为如何？ 本文的模拟显示，在N=50和100时，FDR估计值仍与名义水平重叠，仅在N=200时才显著偏离。
  4. 是否存在更简单的依赖结构（如等相关）也能导致失败？ 本文的因子模型是精心构造的，其一般性有待探索。

⚠️ 作者的framing

作者将缺口frame成“一个被广泛相信但错误的猜想”，通过构造一个反例来“推翻”它。他淡化了以下竞争路线： - 依赖调整方法（Fithian和Lei, 2022; Sarkar, 2023等）被描述为“转向开发替代方法”，而非“对原问题的正面攻击”。作者暗示这些工作回避了核心问题。 - Sarkar (2023) 的部分证明（在已知协方差矩阵下）被提及，但作者的反例是否满足Sarkar的设定？本文的因子模型协方差是已知的（由构造给出），但Sarkar的证明是否覆盖了这种特定结构？作者未明确讨论这一点，值得研究者去查Sarkar (2023) 的原文，看其假设是否排除了本文的反例。

值得查的问题：为什么之前的模拟（Reiner-Benaim, 2007; Kim和van de Wiel, 2008）没有发现这个反例？是因为模拟的样本量不够大，还是因为反例的构造过于特殊？作者在讨论中承认“违反很小”，这暗示了模拟可能缺乏检测这种微小偏离的统计功效。

张力

未见明显对立引用。所有被引工作都承认两尾检验的困难，只是对“BH是否仍能控制”持有不同信念。本文是第一个提供严格反例的。

二、最核心、最简单的例子 / 数学问题

第一步：把符号、模型、可观测数据交代清楚

  • 符号：
  • \( m \)：假设总数（检验个数）。
  • \( \alpha \)：名义FDR水平（本文取0.01）。
  • \( H_i: \mu_i = 0 \)：第i个零假设。
  • \( X_i \)：第i个检验统计量，服从 \( N(\mu_i, 1) \)（边际上）。
  • \( P_i = 2\bar{\Phi}(|X_i|) \)：两尾高斯p值，其中 \( \bar{\Phi} = 1 - \Phi \) 是标准正态生存函数。
  • \( t_r = \alpha r / m \)：BH程序的第r个临界值。
  • \( R = \max\{ r : \#\{i: P_i \leq t_r\} \geq r \} \)：BH拒绝数。
  • \( V = \#\{i \in I_0: P_i \leq t_R\} \)：错误拒绝数（真零假设中被拒绝的个数）。
  • \( \text{FDP} = V / \max(R, 1) \)：错误发现比例。
  • \( \text{FDR} = \mathbb{E}[\text{FDP}] \)：错误发现率。
  • \( Z \)：公共潜因子（标准正态）。
  • \( \varepsilon_i, \eta_j, \xi_k \)：独立标准正态噪声。
  • \( N \)：控制问题规模的整数参数。
  • \( m_N = 100N \)：总假设数。
  • \( \pi_0 = 0.96 \)：真零假设比例。
  • \( \hat{G}_N(t) \)：所有p值的经验CDF。
  • \( \hat{F}_{0,N}(t) \)：真零假设p值的经验CDF。
  • \( G_z(t) \)：条件于Z=z的极限p值CDF。
  • \( F_{0,z}(t) \)：条件于Z=z的极限真零p值CDF。
  • \( u(c) = 2\bar{\Phi}(c) \)：从临界值c到p值的映射。

  • \( Q(c; a, s) = \bar{\Phi}((c-a)/s) + \bar{\Phi}((c+a)/s) \)：\( N(a, s^2) \) 变量的两尾概率。

  • 模型：数据来自一个单因子高斯模型：

    \[X = a Z + \varepsilon,\]
    其中 \( Z \sim N(0,1) \) 是公共因子，\( \varepsilon \sim N(0, D) \) 是独立噪声，\( D \) 是对角矩阵。具体地，三个坐标块有不同的载荷a和均值μ：

  • 块0（96N个真零）：\( X_i^{(0)} = \frac{3}{10}Z + \frac{\sqrt{91}}{10}\varepsilon_i \)，均值0。
  • 块1（N个非零）：\( X_j^{(1)} = \frac{12}{5} - \frac{3}{10}Z + \frac{\sqrt{91}}{10}\eta_j \)，均值12/5。

  • 块2（3N个非零）：\( X_k^{(2)} = \frac{22}{5} - \frac{18}{25}Z + \frac{\sqrt{301}}{25}\xi_k \)，均值22/5。 所有噪声项和Z独立。协方差矩阵 \( \Sigma = aa^T + D \) 是对角线为1的相关系数矩阵。

  • 可观测数据：研究者观测到 \( X_1, \ldots, X_{m_N} \) 的样本，并据此计算p值 \( P_i = 2\bar{\Phi}(|X_i|) \)。不可观测的是潜因子Z和噪声项的具体实现。Z的存在使得所有X_i之间产生相关性，但这种相关性结构（因子模型）是已知的。

第二步：讲最小内核

本文的核心思路可以简化为以下最简特例：考虑一个只有两个块的因子模型： - 块A（大量真零）：\( X_i = \rho Z + \sqrt{1-\rho^2} \varepsilon_i \)，均值0。 - 块B（少量非零）：\( X_j = \mu + \rho Z + \sqrt{1-\rho^2} \eta_j \)，均值μ > 0。

所有变量独立标准正态。当Z取某个较大的正值时（例如z=2），条件于Z=z： - 块A的X_i条件均值为ρz，非零。这意味着真零假设的p值分布不再均匀——它们倾向于产生较小的p值（因为|X_i|倾向于变大），从而更容易被BH拒绝。 - 块B的X_j条件均值为μ + ρz，更大，产生更小的p值，更容易被拒绝。

现在，BH程序在条件于Z=z下运行。由于块A的p值分布向左偏移（更小），经验CDF \( \hat{G}_N(t) \) 在t较小时就超过了t/α，导致BH阈值τ_N被推高。同时，块B的非零p值进一步推高τ_N。结果，τ_N可能超过一个临界点，使得真零假设中被拒绝的比例（即FDP）超过α。

本文的构造比这个特例更精细：它用了三个块（两个信号块，不同强度），并精心调整了载荷和均值，使得在Z的某个取值范围内，条件FDP的下界积分后超过α。核心数学困难在于：需要证明对于所有足够大的N，FDR > α，而不是仅仅在某个极限下。这需要处理有限样本的随机性，并给出一个严格的数值下界。

三、这篇论文做了什么

三句话

  1. 研究了什么问题：Benjamini-Hochberg程序在相关两尾高斯p值下是否控制FDR？具体地，反驳了“BH在任意相关两尾高斯检验下仍控制FDR”这一持续二十年的猜想。
  2. 核心工具/方法：构造了一个特定的高斯因子模型，利用经验CDF crossing表示将BH阈值与极限p值CDF联系起来，并通过区间算术（Arb库） 给出严格的数值证书，证明FDR的下界超过名义水平。
  3. 主要结论：在名义水平α=0.01下，对于足够大的假设数量（m=100N, N足够大），BH程序的FDR严格大于0.0104，即FDR > α。

关键设定与假设

  • 设定：\( X \sim N_m(\mu, \Sigma) \)，其中Σ是相关系数矩阵。p值为两尾高斯p值 \( P_i = 2\bar{\Phi}(|X_i|) \)。BH程序在名义水平α=0.01下运行。
  • 假设：
  • 因子模型结构：Σ = aaᵀ + D，其中a是载荷向量，D是对角正定矩阵。这保证了Σ是相关系数矩阵（对角线为1）。
  • 块结构：三个块，比例分别为0.96（真零）、0.01（非零）、0.03（非零），每个块内条件独立。
  • 均值非零：两个信号块的均值均为正（12/5和22/5）。
  • 载荷符号：真零块载荷为正（3/10），信号块载荷为负（-3/10和-18/25）。这导致当Z增大时，真零块的条件均值增大（p值变小），而信号块的条件均值减小（p值变大），形成一种“对抗”关系。
  • 相比已有文献：本文的假设是构造性的，而非一般性的。它不试图覆盖所有相关结构，而是精心设计一个使BH失败的具体例子。这与Sarkar (2023) 的“已知协方差矩阵”设定不同——本文的协方差矩阵是已知的，但Sarkar的证明是否覆盖了这种特定结构？作者未明确讨论。

主要结果

定理1：对于每个整数N≥1，令m_N = 100N，前96N个假设为真零，其余4N个为非零（按上述因子模型）。则BH程序在α=0.01下满足：对所有足够大的N，FDR_N > 0.0104 > α。

  • 直觉：条件于Z=z，BH阈值τ_N被两个信号块产生的极小p值推高，同时真零块的p值分布因Z而左移（更易被拒绝）。在Z的某些取值下，这种效应导致FDP超过α。通过对Z的分布积分，得到无条件FDR的下界。
  • 必要条件：N足够大，使得经验CDF收敛到极限CDF（Glivenko-Cantelli定理）。反例的构造依赖于大样本渐近。
  • 解决的技术难点：需要将无限多个Z的取值和连续c域上的不等式验证，转化为有限个有理区间上的严格数值验证。这通过单调性（Q函数在c和|m|上的单调性）和区间算术实现。

证明路线与技术技巧

整体路线（3-5步逻辑主干）：

  1. 条件极限CDF：条件于Z=z，三个块内独立同分布。由Glivenko-Cantelli定理，p值的经验CDF \( \hat{G}_N(t) \) 和真零p值的经验CDF \( \hat{F}_{0,N}(t) \) 分别一致收敛到极限CDF \( G_z(t) \) 和 \( F_{0,z}(t) \)。

  2. BH阈值的bracketing：引理2和3将BH阈值τ_N和FDP的下界与极限CDF的符号条件联系起来。具体地，如果存在v和w使得 \( G_z(v) > v/\alpha \) 且对所有t∈[w,α]有 \( G_z(t) < t/\alpha \)，则 \( v \leq \liminf \tau_N \leq \limsup \tau_N \leq w \)，且 \( \liminf \text{FDP}_N \geq \alpha \pi_0 F_{0,z}(v)/w \)。

  3. 构造满足符号条件的v和w：对于每个z的区间B_k = [k/100, (k+1)/100]，利用Q函数的单调性，将连续c域上的不等式验证转化为有限个网格点上的严格不等式。具体地：

  4. 找到第一个网格点a_k，使得对所有c∈[c_α, a_k]和z∈B_k，有 \( h_z(c) = G_z(u(c)) - 100u(c) < 0 \)（即G_z(t) < t/α）。这通过上界U_{j,k} < 0来验证。
  5. 找到第一个网格点b_k ≥ a_k，使得对所有z∈B_k，有 \( h_z(b_k) > 0 \)（即G_z(u(b_k)) > 100u(b_k)）。这通过下界L_{ℓ_k,k} > 0来验证。

  6. 令v = u(b_k), w = u(a_k)，则满足引理2和3的条件。

  7. 条件FDP下界：由引理3，对每个z∈B_k，有 \( \liminf \text{FDP}_N \geq d_k \)，其中 \( d_k = 0.0096 Q(b_k; m_{0,k}^-, s_0) / u(a_k) \)。

  8. 积分得到无条件FDR下界：利用Fatou引理，将条件下界对Z的分布积分：

    \[\liminf_{N\to\infty} \text{FDR}_N \geq \sum_{k=-500}^{499} d_k \mathbb{P}(Z \in B_k) > 0.0104.\]
    最后一个不等式由区间算术证书严格验证。

关键跳跃点： - 从连续到离散的转化：最吃功夫的部分是将无限多个(z,c)对上的不等式验证，转化为有限个有理网格点上的严格数值验证。这依赖于Q函数在c和|m|上的单调性，以及u(c)的单调性。单调性保证了：如果上界U_{j,k}在网格点(c_j, c_{j+1})×B_k上为负，则整个矩形区域上的h_z(c)都为负。 - 区间算术证书：所有数值计算（高斯尾部概率、平方根等）都用Arb库的球算术（ball arithmetic）进行，保证每个计算结果都包含真实值。严格比较只在整个球位于零的某一侧时才接受。这提供了数学上严格的证明，而非数值近似。

技术技巧点名： - 经验过程/Glivenko-Cantelli定理：用于建立条件经验CDF的一致收敛性。 - 单调性分析：Q函数在c和|m|上的单调性，用于将连续域上的不等式验证简化为有限网格点。 - 区间算术（Arb库）：用于提供严格的数值证书，避免浮点误差。 - Fatou引理：用于将逐点（条件）下界积分得到无条件下界。 - 分层/分箱：将Z的连续分布离散化为1000个小区间，每个区间内用极值（m⁻, m⁺）来bound条件CDF。

真实例子与应用

本文包含一个蒙特卡洛实验，用于支持理论结果。

  • 数据/场景：按照定理1的因子模型生成数据，取N=50, 100, 200（对应m=5000, 10000, 20000）。使用分层抽样（对Z进行1000个等概率分层）来降低方差。
  • 方法应用：对每个生成的完整数据集，计算两尾p值，应用普通BH程序（α=0.01），记录FDP。通过100次宏复制（每次包含1000个分层样本）估计FDR。
  • 结果：
  • N=50: 估计FDR=0.009936，95% MC区间[0.009711, 0.010161]，与α无显著差异。
  • N=100: 估计FDR=0.010129，区间[0.009939, 0.010320]，仍与α重叠。
  • N=200: 估计FDR=0.010359，区间[0.010155, 0.010563]，完全位于α之上。单侧p值=3.56×10⁻⁴。
  • 说明的问题：实验验证了理论结果，即对于足够大的N，FDR确实超过α。同时显示，在较小样本量下，偏离可能不显著，这与理论要求的“N足够大”一致。实验也展示了分层抽样在降低蒙特卡洛方差方面的有效性。

🔎 结论是否比证明窄

  • 窄的方面：定理1只证明了α=0.01下的反例，且FDR超过量很小（约3.59%）。作者在讨论中明确承认“violates the nominal level only slightly”，并问“是否存在通用上界”。这意味着结论可能不适用于其他α水平（如0.05或0.1），或者偏离可能被一个与α相关的上界所限制。
  • 泛化的claim：作者在摘要和引言中声称“BH程序可以失败”，但证明仅针对一个特定构造。虽然这足以反驳“任意相关下都控制”的猜想，但并未说明失败在何种条件下普遍发生。作者在讨论中提出“是否对较小数量的检验也失败”作为开放问题，暗示了有限样本下的结论可能不同。
  • 具体语句：定理1的陈述是“for all sufficiently large N”，而模拟显示N=50和100时未检测到显著偏离。这暗示“sufficiently large”可能相当大（N≥200），且有限样本行为可能不同。

四、开放问题

  1. FDR偏离的通用上界：本文的反例显示FDR超过α约3.59%，但这是否是可能的最大偏离？是否存在一个与α和依赖结构有关的通用上界？作者在讨论中明确提出了这个问题：“whether a universal bound exists on the possible inflation of the FDR above its nominal level”。（扎根于Section 8第一段）

  2. 有限样本下的行为：本文的模拟显示，在N=50和100时，FDR估计值与α无显著差异。这是否意味着BH程序在有限样本下仍然控制FDR？或者只是统计功效不足？作者提出“determine whether the BH procedure is guaranteed to control the FDR for smaller numbers of tests”作为未来方向。（扎根于Section 8第二段）

  3. 其他α水平下的反例：本文仅针对α=0.01构造了反例。对于其他常见的α水平（如0.05, 0.1），是否存在类似的反例？构造方法是否可推广？作者未讨论这一点，但这是自然的延伸。

  4. 更一般的依赖结构：本文的因子模型是精心构造的。是否存在更简单或更自然的依赖结构（如等相关、AR(1)）也能导致BH失败？或者，本文的反例是否揭示了某种必要条件？这需要进一步的理论分析。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论