跳转至

Empirical partially Bayes multiple testing and compound χ2 decisions

作者: Nikolaos Ignatiadis, Bodhisattva Sen
来源: Annals of Statistics
主题: 数理统计 / 假设检验
相关性: 7/10
链接: 期刊页 · arXiv

一、领域脉络与小综述

  • 这个方向是什么 这个子方向是多重假设检验中的经验贝叶斯（Empirical Bayes, EB）方法，特别是针对高斯数据下均值检验时，如何处理未知方差的“高维”/“多单位”问题。其根本统计问题是：在并行检验数千个基因/蛋白单元时，每个单元的观测数据（样本均值、样本方差）都来自一个未知的均值-方差对，如何利用所有单元的信息来改进每个单元的检验（即借力，borrowing strength），且能经得起错误发现率（FDR）的控制。当前成熟度较高，但有一个理论缺口——对均值和方差的贝叶斯处理不对称。

  • 发展脉络（history）

    • 奠基工作（复合决策理论, 1950s-1960s）：Robbins (1951, 1956) 奠定了复合决策理论 (compound decision theory)，提出“同时借用多个独立但相似决策问题的信息来改进每个决策”。这是EB的思想源头。
    • 主要进展一：参数经验贝叶斯与“调节方差” (1990s-2000s)：在微阵列数据爆发后，Cui, Hwang,等人（Cui-Hwang 1997, Hwang-Peddeti 2009）以及Limma包（Smyth 2004）将EB具体化为：对每个基因的样本方差施加先验（通常是逆伽马分布），然后通过“调节”（moderate）或“收缩”（shrink）方差，计算出t-统计量或修正后的F统计量。原文提到：“It is routine to 'moderate' (i.e., to shrink) the sample variances through parametric empirical Bayes methods before computing p-values for the means.” 这种做法的特点是：对方差（nuisance参数）应用贝叶斯，但对均值（primary参数）不用——是一种非对称的部分贝叶斯。
    • 主要进展二：部分贝叶斯推断的理论 (Cox, 1975)：David Cox 在 1975 年提出，若对nuisance参数有先验，则可以基于其充分统计量做条件推断（conditional inference），以消除nuisance参数影响。原文将Cox (1975) 视为“partially Bayes inference”的倡议，并作为其方法论的直接前身。
    • 当前 frontier：如何将这种不对称的、对nuisance参数的经验贝叶斯处理严格形式化，并证明其与主流FDR控制工具（Benjamini-Hochberg, BH）的兼容性。作者将其称为“empirical partially Bayes multiple testing”（经验部分贝叶斯多重检验），并声称是“initiates the formal study of this paradigm”。
    • 本文位置：作者将之前的“调节方差”做法（如limma）视为一种启发式方法，并试图用非参数最大似然（NPMLE）等现代工具，为这种部分贝叶斯策略提供严格的渐近理论，包括计算条件p值的Eddington/Tweedie公式、NPMLE的收敛速度，以及渐近FDR控制。

  • 子线索聚类

    1. 参数经验贝叶斯收缩与t-统计量：Cui-Hwang (1997), Limma (Smyth 2004)。这簇方法假设方差先验有参数形式（常为逆伽马），并据此调整检验统计量。优点是计算简单、直观；缺点是先验设定可能错误，且缺少对均值先验的对称处理。
    2. 非参数经验贝叶斯（NPEB）：这是作者借鉴的核心工具。Brown & Greenshtein (2009), Koenker & Mizera (2014), Efron (2016), Polyanskiy & Wu (2021) 等文献在均值的估计/预测问题中，用NPMLE来估计先验，并证明其良好性质（如收敛速率）。作者将其推广到对方的先验估计。
    3. 复合决策与FDR控制：Efron (2004), Jiang & Zhang (2009), Zhang (2017), 以及经典的 BH (1995) 程序。这簇关注在“大量同时决策”场景下如何控制错误率。作者的核心贡献之一是证明其条件p值可与BH联用，这直接将部分贝叶斯方法与现成的多重比较工具箱对接。

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

    • 核心问题1：在多重检验中，如何利用所有单元的信息来改进每个单元的检验，同时控制全局FDR？
    • 核心问题2：对于先验的不对称（对nuisance参数用、对primary参数不用），我们应该进行怎样的理论理解？它相比完全的联合EB（均值和方差都用先验）有何优点（无需对均值假设/计算后验分布）或代价（效率损失）？
    • 核心问题3：NPMLE作为先验估计器，在这种部分贝叶斯设定下，其收敛速度如何？以及这种速度如何影响下游的FDR控制？ 本文的核心技术贡献是回答了这个问题：NPMLE对方差先验的估计达到“近乎参数”（near-parametric）速率（定理3）。

  • ⚠️ 作者的 framing（必须明确标注成“这是作者的说法”）

    • 作者把缺口 frame 成什么？ 作者声称，尽管“调节方差”在领域内已是一种“routine”（常规操作），但它缺乏一个严格的、正式的理论框架（“Our work initiates the formal study of this paradigm”）。他们将这个框架命名为“empirical partially Bayes multiple testing”。这使其论文成为“显然的下一步”——从启发式操作上升到严格理论。
    • 哪些竞争路线被他淡化或回避了？ 作者明确提到“In our work, we restrict ourselves to ... the Gaussian sampling model where both means and variances are unknown”以及“Conditional on the sample variances”，这回避了同时包含均值和方差的完全经验贝叶斯（joint EB）。一个明显被淡化的路线是：对均值和方差都用NPMLE。为什么不用？作者在引言中暗示，部分贝叶斯方法只“修饰”了nuisance参数，而对primary参数（均值）使用了基于条件似然的频率学推理。这可能是因为在均值检验中，直接对均值施加先验可能引入偏差，且计算条件p值在技术上有优势。
    • 什么明显该被引 / 该存在、却没出现在intro里？ 引言中提到了Robbins（1951,1964）、Cox（1975）、Cui-Hwang（1997）、Smyth（2004）、Efron（2004）、Hwang-Peddeti（2009）、Brown-Greenshtein（2009）、Koenker-Mizera（2014）、Zhang（2017）、Polyanskiy-Wu（2021）等。没有提到的是：有关“高维稀疏信号”的检验理论（如Spike-and-Slab先验，或者Donoho-Jin的Higher Criticism）。本文处理的是“所有单元都有非零方差”的复合决策场景，而非稀疏场景。这可能是一个重要的缺口（值得研究者去查：如果信号是稀疏的，这个部分贝叶斯框架还能成立吗？）。
    • 张力：未见明显对立引用。所有被引工作基本都支撑或铺垫了作者的框架。

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

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

  • 符号

    • \( i = 1, \dots, m \)： 单元（如基因）的索引。
    • \( n \)： 每个单元内的样本量。注意：本文假设每个单元的样本量相同（balanced design）。这是关键假设。
    • \( X_i \)： 第 \( i \) 个单元的样本均值（scalar）。
    • \( V_i \)： 第 \( i \) 个单元的样本方差（scalar）。它自由度为 \( n-1 \) 的卡方分布。
    • \( \mu_i \in \mathbb{R} \)： 第 \( i \) 个单元的未知总体均值（primary参数）。
    • \( \sigma_i^2 > 0 \)： 第 \( i \) 个单元的未知总体方差（nuisance参数）。
    • \( G \)： \( \sigma_i^2 \) 的未知先验分布。这是本文的核心被估对象。这是一个 CDF (Cumulative Distribution Function)，支撑在 \( (0, \infty) \)。
    • \( P_{val,i} \)： 对第 \( i \) 个单元检验 \( H_{0,i}: \mu_i = 0 \) 的p值。
    • BH： Benjamini-Hochberg 程序，用于在多重比较中控制FDR。

  • 模型

    • 数据生成机制：对于每个单元 \( i \)，其数据点来自 \( \mathcal{N}( \mu_i, \sigma_i^2 ) \)。根据标准理论：
      • 样本均值和样本方差联合分布：\( X_i \mid \mu_i, \sigma_i^2 \sim \mathcal{N}( \mu_i, \sigma_i^2 / n ) \)， \( (n-1) V_i / \sigma_i^2 \sim \chi^2_{n-1} \)。
      • 且 \( X_i \) 与 \( V_i \) 在给定 \( (\mu_i, \sigma_i^2) \) 下条件独立。这是本文推导p值的关键。
    • 模型结构：总体 \( \sigma_i^2 \) 独立同分布地来自未知分布 \( G \)，即 \( \sigma_i^2 \overset{i.i.d.}{\sim} G \)。而 \( \mu_i \) 没有被假设服从任何先验；它们是固定但未知的（fixed, unknown）。本文不估计/推断 \( \mu_i \) 的先验，而是专注于在给定 \( G \) 的条件下，基于观测到的 \( (X_i, V_i) \) 对 \( \mu_i \) 做条件推断。
    • 已知：样本量 \( n \) 已知；自由度 \( n-1 \) 已知；检验问题 \( H_{0,i}: \mu_i = 0 \)；BH程序门限 \( q \) 已知。
    • 要估的对象：先验 \( G \)。

  • 可观测数据

    • 我们能观测到：每个单元i的样本均值 \( x_i \) 和样本方差 \( v_i \)。 这就是全部数据。我们有一张 \( m \times 2 \) 的数据表。
    • 想要但观测不到的：总体均值 \( \mu_i \)（我们想要检验它）；总体方差 \( \sigma_i^2 \)（它被 \( G \) 估计）；以及 \( G \) 本身。
    • 关键识别条件：我们可以通过 \( m \) 个单位的样本方差 \( V_1, \dots, V_m \) 来估计 \( G \)。这是因为 \( (n-1)V_i \) 在给定 \( \sigma_i^2 \) 下是 \( \sigma_i^2 \chi^2_{n-1} \)。NPMLE 是用于从这种带噪声的观测中恢复 \( G \) 的标准工具。本文依赖NPMLE对 \( G \) 的可行性（这在理论中被隐含假设）。

第二步：讲最小内核

  • 最简特例：考虑一个最简单的问题：一个单元 \( i \)，其样本量 \( n \) 足够大，我们先忽略NPMLE，仅考虑已知 \( G \) 的情况。这是所有复杂性的最小内核：当先验已知时，如何基于观测到的 \( (x_i, v_i) \) 计算一个p值，使得它的类型I错误率是均匀的？
    • 在这个特例下，本文核心命题退化成：
      • 目标：检验 \( H_{0,i}: \mu_i = 0 \)。观测到 \( (x_i, v_i) \)。
      • 我们想做什么：想构造一个函数 \( p_i = f(x_i, v_i; G) \)，使得在 \( \mu_i=0 \) 下，\( p_i \) 的分布是 \( \text{Unif}(0,1) \) （至少是随机变量均匀分布，在条件概率空间中）。
      • 经典方法（无先验）：用t-统计量 \( t = \sqrt{n} x_i / \sqrt{v_i} \)，p值 = \( P( |T_{n-1}| > |t| ) \)。这完全忽略了 \( V_i \) 的分布信息（但分布已知）。
      • 部分贝叶斯方法（Cox, 1975）的思想：因为 \( X_i \) 与 \( V_i \) 在给定 \( \sigma_i^2 \) 下独立，且 \( X_i \) 的分布完全由 \( (\mu_i, \sigma_i^2) \) 决定，我们可以这样想：如果我们知道 \( V_i \)（即观测到了样本方差），那么我们实际知道了关于 \( \sigma_i^2 \) 的一切信息（通过似然）。更关键的是，我们能写出 \( X_i \) 在给定 \( V_i \) 和 \( \mu_i=0 \) 下的分布。这个分布不再依赖于 \( \sigma_i^2 \)了。为什么？因为：
        • 在 \( H_0: \mu_i=0 \)下，\( X_i \sim \mathcal{N}(0, \sigma_i^2/n) \)。
        • 且 \( (n-1)V_i / \sigma_i^2 \sim \chi^2_{n-1} \)。
        • 联合分布 \( (X_i, V_i) \)在给定 \( \sigma_i^2 \) 下是独立的。
        • 但是，当我们条件于 \( V_i = v_i \) 时，我们实际上积分掉了 \( \sigma_i^2 \)。具体来说，给定 \( V_i \)， \( \sigma_i^2 \) 的后验分布完全由 \( G \) 和似然 \( p(V_i \mid \sigma_i^2) \) 决定。然后 \( X_i \) 在 \( H_0 \) 下的分布，是 \( \mathcal{N}(0, \sigma_i^2/n) \) 对这个后验的混合。
        • 所以，条件p值 = \( P( |X_i| > |x_i| \mid V_i=v_i, \mu_i=0, G) \)，这个概率只依赖于 \( G \)、\( x_i \)、\( v_i \) 和 \( n \)。在已知 \( G \) 下，这个条件p值可以精确计算（或通过积分）。
      • 为什么这是“内核”：这个条件p值不需要对 \( \mu_i \) 的先验，只用了 \( \sigma_i^2 \) 的先验（通过G）。它完美体现了Cox的部分贝叶斯思想。本文后面做的一切（NPMLE估计G、Eddington/Tweedie公式）都是在为这个“已知G下的条件p值”提供一个可操作的、且具有良好收敛性的替代。 当 \( G \) 未知时，他们用NPMLE从全部观测 \( v_1,\dots,v_m \) 中估计 \( G \)，然后代入计算；然后证明这个估计的p值渐近地与原条件p值一致。

三、这篇论文做了什么

  • 三句话 ① 研究了在多重假设检验中，当每个单元的均值和方差都未知时，如何通过仅对方差施加经验贝叶斯（NPMLE）来系统化“部分贝叶斯推断”这一非对称范式。 ② 核心工具/方法是条件p值（在给定样本方差下对均值进行检验），并证明了它满足一个Eddington/Tweedie公式（定理1/2），使得p值的计算可化为对先验的估算。当先验由非参数最大似然（NPMLE）估计时，p值的收敛速度达到近乎参数（near-parametric）（定理3）。 ③ 主要结论：这些估计出的条件p值可与Benjamini-Hochberg (BH) 程序相结合，在渐近意义上控制错误发现率（FDR）（定理5）。即使在方差为固定常数的复合决策设定下，该方法仍保留渐近类型I误差保证（定理5推论）。

  • 关键设定与假设

    • 设定：\( m \) 个独立单元，每个单元有 \( n \) 个独立同分布的高斯观测 \( Y_{ij} \sim \mathcal{N}(\mu_i, \sigma_i^2) \)（其中 \( j=1,\dots,n \)）。\( n \) 对所有单元相同。这引出了 \( X_i \) 和 \( V_i \)。目标是多重检验 \( H_{0,i}: \mu_i =0 \) 对 \( H_{1,i}: \mu_i \neq 0 \)。
    • 假设：
      • 随机效应假设：\( \sigma_i^2 \) 服从一个未知的、定义在 \( (0, \infty) \) 上的分布 \( G \)。这是将单位作为来自某个总体的随机样本（这不同于传统将单位视为固定效应的模型）。
      • 平滑假设：定理3收敛速度的证明需要假设 \( g(\cdot) = \log p(\cdot) \)（其中 \( p \) 是观测到的 \( V \) 的边缘密度）是“适度平滑的”（moderately smooth），例如属于Hölder类。定理3的速率依赖于平滑参数。
      • 边界假设：定理3要求 \( V \) 的分布不集中在0或无穷远附近，即需要保证一些矩条件。这包括假设 \( G \) 的支撑不退化。
      • 独立性假设：\( m \) 个单元之间独立；且给定 \( (\mu_i, \sigma_i^2) \) 下，\( X_i \) 和 \( V_i \) 独立。这些都是标准的高斯模型假设。
      • BH兼容性假设：定理5的证明依赖于条件p值的“渐近独立性”或“渐近正相依性”（PRDS）条件（该条件在引言中被提及）。这是将BH作用于任何依赖p值的测试的基础。原文在定理5的陈述中使用了“asymptotic independence”（渐近独立性）。
    • 相比已有文献：相比参数EB（如Limma）用了非参数先验；相比Cox的原始部分贝叶斯（假设先验已知）引入了经验成分；相比Robbins的原始经验贝叶斯（通常处理单个参数的估计/检验），本文将焦点放在多重比较和FDR控制上。

  • 主要结果

    • 定理2（已知G下的条件p值均匀性）：如果 \( G \) 已知，则条件p值 \( \hat{P}_i = P( |X| > |x_i| \mid V_i=v_i, \mu_i=0, G) \) 在 \( H_0 \) 下是均匀分布的。这是条件频率学派方法的基石（Cox观点）。
    • 定理3（NPMLE的收敛速度）：设 \( \hat{G}_{\text{NPMLE}} \) 是从 \( (V_1,\dots,V_m) \) 中得到的NPMLE估计。在一定的平滑假设下，对于任何 \( k\in\mathbb{N} \)，有
      \[\sup_{t} | \hat{F}_G(t) - F_G(t) | = O_p( m^{-(k)/(2k+1)} \log m )\]
      其中 \( F_G \) 是 \( G \) 的累积分布函数，\( k \) 是平滑参数。这意味着收敛速度是 “近乎参数”（即接近 \( m^{-1/2} \)），取决于平滑度，但不依赖于 \( n \)。这解释了为什么部分贝叶斯方法非常有效——来自其它单元的数据（大的 \( m \)）能帮助精炼对 \( G \) 的估计，从而改进每个单元的条件p值。该速率的证明依赖于NPMLE的有限样本浓度界（如SRE理论, Polyanskiy & Wu 2021）以及Donsker类性质。
    • 定理5（FDR控制）：设 \( \hat{P}_{(1)} \le \dots \le \hat{P}_{(m)} \) 是使用NPMLE估计的p值。使用BH程序在名义水平 \( q \in (0,1) \) 下拒绝所有 \( H_{0,(i)} : i \le k_{\max} \)，其中 \( k_{\max} = \max\{ k: P_{(k)} \le q \cdot k/m \} \)。则在温和条件下，当 \( m\to\infty \) 且 \( n \) 固定时，
      \[\limsup_{m\to\infty} \text{FDR} \le q .\]
      这个结果不要求 \( \mu_i \) 的先验，甚至不要求 \( \mu_i \) 是随机的（它可以是固定的）。即使在复合决策设定下（方差是固定的，不是来自G），该渐近保证也成立，这被称为“FDR robustness to model misspecification of the prior”。

  • 证明路线与技术技巧（理论型必写）

    • 整体路线：
      1. 证明条件p值的核心性质（已知G）：通过推导Eddington/Tweedie公式，将条件p值的计算转化为对边际密度 \( m(v) \) 和 \( m(v, x) \) 的积分，进而证明其均匀性。这利用了 \( X_i \) 与 \( V_i \) 的条件独立性。
      2. 估计G的正面性质：使用NPMLE从观测到的 \( V_i \) 中估计 \( G \)。证明NPMLE的收敛速度。关键是利用 \( V_i \) 的边缘密度 \( p(v; G) = \int_0^\infty f(v; \sigma^2) dG(\sigma^2) \) 这一混合分布的似然。NPMLE在混合分布上的收敛性由Srebro, Rakhlin, 和 Evgeniou (2008)以及Polyanskiy & Wu (2021)的著名结果保证。
      3. p值的误差传播分析：将第2步中 \( \hat{G} \) 的误差传播到p值估计。需要证明，如果 \( \hat{G} \) 在Kolmogorov-Smirnov距离下收敛，那么用它计算的条件p值 \( \hat{p}_i \) 与真实条件p值 \( p_i \) 之间的偏差也以相应速率收敛。这一步需要积分运算的连续性，并通过Donsker类和积分逼近实现。
      4. FDR控制证明：在渐近框架下，证明用 \( \hat{p}_i \) 的BH程序的FDR不大于 \( q \)。标准BH结果（Benjamini & Hochberg 1995）需要p值在 \( H_0 \) 下是均匀的和/或正相依。第1步证明了 \( p_i \) 是均匀的。第3步证明 \( \hat{p}_i \) 接近 \( p_i \)。因此，通过一个“渐近等值”的论证，加上条件p值的“渐近独立”性质，可以将BH的有限样本FDR控制结果推广到渐近意义。
    • 关键跳跃点：
      1. 从“已知G”到“未知G”的跃迁：这是最关键的。不能简单地“plug-in” \( \hat{G} \)，因为这会引入偏差。作者如何做到？ 他们证明了NPMLE的收敛速度足够快（近乎参数），使得“plug-in”后的p值偏差可以被忽略（以 \( m^{-1/2+} \) 速率收缩）。这利用了NPMLE在混合模型中的强大性质（如收敛速度为根号n乘以log因子）。这并非显然的，因为即使 \( \hat{G} \) 收敛，p值的计算也可能对 \( G \) 的微小变化非常敏感。
      2. 条件p值的独立性（FDR控制的关键）：即使在高维设置下，条件p值之间是否独立？它们从同一个 \( \hat{G} \) 中计算。作者假设并利用了渐近独立性性质。文中可能用了“exchangeability”（可交换性）的概念。证明这个点是工程难点。
      3. 复合设定下的稳健性：定理5声称在方差是固定而非随机时也适用。这是通过证明当方差固定时，条件p值仍渐近均匀，并且其联合结构满足BH所需的依赖条件。作者似乎在证明，即使模型被错误指定（假设方差随机，但实际是固定的），FDR控制仍然成立。这体现了该方法的鲁棒性。
    • 技术技巧点名：
      • Eddington/Tweedie公式：用于计算后验均值和p值。定理2用E/T公式形式化了条件p值。
      • NPMLE（非参数最大似然）：估计未知的 \( G \)。
      • Donsker类的性质：用于证明差商 \( \hat{G} - G \) 的弱收敛以及p值的收敛性。
      • Benjamini-Hochberg程序：用于FDR控制。
      • glrVc理论 (e.g., Polyanskiy & Wu 2021)：用于证明NPMLE的收敛速度。

  • 真实例子与应用 论文包含一个真实数据例子：来自癌症基因组图谱（TCGA）的mRNA表达数据。具体来说，他们将每个基因（~20,000个）视为一个单元，将每个单元的样本均值与方差视为已观测。他们应用经验部分贝叶斯条件p值，并用于BH选择差异表达的基因。他们展示了与limma（一种常用的参数经验贝叶斯方法）的比较，表明二者结果高度一致，验证了方法的实用性。 这个例子主要是为了展示可操作性和与现有方法的兼容性，而非展示强大的性能优势。

  • 🔎 结论是否比证明窄

    • 是的。定理3的收敛速度依赖于平滑假设（即G的分布必须在某个Hölder空间内）。这意味着，如果真实的 \( G \) 是“不光滑”的（例如，集中在两个点上，或具有尖峰），NPMLE的收敛速度可能会更慢。论文的收敛速率结论是限定在平滑先验下的。但论文在应用部分（如与limma的比较）并未检验这个假设。结论的泛化（对所有 \( G \) ）比证明的要宽。
    • 定理5的FDR控制是渐近的。实际中，m可能不够大，导致FDR控制不严格。但论文没有给出有限样本的FDR界。

四、开放问题（点到为止）

  • 问题1：收敛速度的最优性。定理3的速率 \( m^{-4/5} \) 是在平滑度 \( k=2 \) 下的结果。对于NPMLE，是否能达到 \( m^{-1/2} \)（即参数速率）？这是一个理论问题。扎根于：定理3的陈述和证明中对平滑度的隐含依赖。
  • 问题2：NPMLE的计算瓶颈。NPMLE的求解是凸问题，但它需要迭代计算，对于 \( m \) 很大的场合（如百万级数据）可能计算成本高。是否有更快的替代方案（如基于核的、更直接的方法）？扎根于：论文引言和讨论中可能提到的计算挑战。
  • 问题3：非平衡设计。论文假设每个单元有相同的样本量 \( n \)。如果样本量不同，还有类似的理论吗？条件p值的计算将更复杂（因为自由度不同），NPMLE的收敛性也可能改变。这是一个明显的泛化方向，扎根于论文的假设（“balanced design”）。
  • 问题4：与稀疏信号的适配。论文处理的是“所有单元都有非零方差”的密集场景。在信号非常稀疏（例如只有5%的单元有非零均值）的情况下，FDR控制可能会变差（因为BH基于所有p值）。如何合并稀疏信号先验（如Spike-and-Slab）并仍然使用部分贝叶斯框架？扎根于：论文在对均值的处理上回避了先验，但研究稀疏信号时，不基于先验可能效率不足。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论