18.7 特征评估和多重检验问题

本章前面的部分先介绍了 $p\gg N$ 场景下的预测问题。本节将讨论一个更基本的问题，如何评估 $p$ 个特征变量中每一个特征的显著程度。回到第 18.4.1 节中蛋白质质谱的例子。在那个问题中，科学家们可能关注的并不是预测某个患者是否患有前列腺癌症。他们可能想要做的是识别出那些在正常和癌症细胞中差异比较大的蛋白质，从而可以加强对这个疾病的认知，并获得药物研发的方向。所以我们的目的就是评估各个特征变量的显著性。这个特征评估通常并不是通过像本章之前那些多元预测模型来完成的。特征评估问题将我们的目光从预测问题转移到了 多重假设检验（multiple hypothesis testing） 这个经典的统计学课题上。在本章下文中将使用 $M$ 而不是 $p$ 来代表特征的数量，为了避免与下文中经常提到的 p 值（p-value）混淆。

例如表 18.4 中的微阵列数据，是来自一项癌症患者对电离辐射治疗（ionizing radiation treatment）敏感度的研究（Rieger et al., 2004）。每一行包含了 58 个患者样本的基因表达水平测量：其中 44 个患者样本有正常的反应，而 14 个患者样本对辐射有严重的（不良）反应。数据是对寡核苷酸（Oligonucleotide）微阵列的测量。试验的目的是找出在对辐射敏感的患者组中表达水平有差异的那些基因。一共有 $M=12,625$ 个基因；表格中为了演示展示出了一部分基因和样本的数据。

为识别出携带信息的基因，我们为每个基因构造一个双样本 t 统计量（two-sample t-statistic）：

其中的 $\bar{x}_{kj}=\sum_{i\in{C_\ell}}x_{ij}/N_\ell$.这里的 $C_\ell$ 是样本组 $\ell$ 中的 $N_\ell$ 个样本的索引下标，$\ell=1$ 代表正常组、$\ell=2$ 代表敏感组。数值 $\text{se}_j$ 是基因 $j$ 的合并组内标准误差（pooled within-group standard error）：

图 18.18 中的橙色部分是 12,625 个 t 统计量的频率直方图，它们的取值范围从 -4.7 到 5.0。如果 $t_j$ 取值是正态分布的，就会将绝对值大于 2 的统计量认为是显著大的。这应该对应着大约 5% 的显著性水平。此例中有 1189 个基因满足 $|t_j|\geq2$。但是从 12,625 个基因中在概率上来说我们本应该得出很多个显著大的 t 统计量，甚至即便是分组与任意的基因都没有相关性。例如假设所有基因都是独立的（实际上它们当然不是独立的），那么伪显著（falsely significant）的基因数量应该是一个二项分布，其均值为 $12625\cdot0.05=631.3$、标准差为 24.5；对比之下实际的数量 1189 就太小了。

那么要如何评估所有这些 12,625 个基因的结果？这就是所谓的 多重检验（multiple testing） 问题。首先为每个基因计算出对应的 p 值。若假设特征变量是正态分布的，则可通过理论的 t 概率分布来计算 p 值。另一个很好的选择是使用置换分布（permutation distribution）来计算，因而就避免了对数据概率分布的假设。（从原则上来说）需要计算对样本标签所有的 $K=\binom{58}{14}$ 个置换（排列），并计算出每个置换的 t 统计量 $t_j^k$。则基因 $j$ 的 p 值为：

不过当然 $\binom{58}{14}$ 是个很大的数字（大概 $10^13$），所以枚举出所有可能的置换是不可能的。作为替代，我们只从可能的置换中选取一个随机样本；这里使用了 $K=1000$ 个置换的随机样本。

样本中基因（特征变量）的取值都是相似的（例如，都是以同样尺度下进行测量的）；利用这一点，我们可以在计算 p 值时将所有基因的结果合并在一起（作为置换分布）：

由于在合并的零分布（null distribution）比在每个基因各自的零分布中的取值样本要多很多，所以上式可得出一个比式 18.40 更细颗粒度的 p 值。

利用这组 p 值，我们想要对所有的 $j=1,2,\dots,M$ 进行假设检验：

如果 $p_j<\alpha$，则在置信水平 $\alpha$ 时拒绝原假设 $H_{0j}$。这个检验的第一类错误（type I error）率等于 $\alpha$；也就是说，错误地拒绝 $H_{0j}$ 的概率是 $\alpha$。

现在如果考虑到多个检验，就不太确定应该使用什么指标来衡量总体的错误（率）。令 $A_j$ 代表 $H_{0j}$ 被错误拒绝的事件；则从定义上 $\operatorname{Pr}(A_j)=\alpha$。FWER（family-wise error rate） 是至少一个错误拒绝的概率，它是一个被广泛使用的对总体错误率的衡量指标。具体来说，若 $A=\cup_{j=1}^MA_j$ 代表至少一个错误拒绝的事件，则 FWER 就是 $\operatorname{Pr}(A)$。一般来说当 $M$ 比较大时 $\operatorname{Pr}(A)\gg\alpha$，并且它取决于检验之间的相关程度。如果检验之间时独立的而且每个第一类错误率都是 $\alpha$，那么这组检验的 FWER 就是 $(1-(1-\alpha)^M)$。而另一方面，如果检验之间存在正相关性（positive dependence），即 $\operatorname{Pr}(A_j|A_k)>\operatorname{A_j}$，那么 FWER 会小于 $(1-(1-\alpha)^M)$。正相关性在实际问题中经常出现，尤其是在基因组的研究中。

多重检验的一个最简单的方法是 邦费罗尼（Bonferroni） 方法。为了能使 FWER 最大等于 $\alpha$，它让每个单独的检验变得更严格：当 $p_j<\alpha/M$ 时 拒绝 $_{0j}$。易证明这样得出的 FWER 小于等于 $\alpha$（练习 18.16）。当 $M$ 相对不太大时，邦费罗尼方法可能是有效的。但是当 $M$ 比较大时，它会过于保守，即经它认定为显著的基因过少。

在上面的例子中，假如检验的置信水平为 $\alpha=0.05$，那么对应 p 值的阈值就是 $0.05/12625=3.9\times10^{-6}$。但是在 12,625 个基因中没有一个的 p 值小于这个树值。

这个方法的一些变体调整了各自 p 值从而达到最大 $\alpha$ 的 FWER，它们使用了一些规避了独立性假设的方法。例如可参考 Dudoit et al. (2002b)。

18.7.1 假发现率

多重检验另一个不同的方法并不会试图控制住 FWER，而是关注了假显著的基因的比例。这个方法在实践中有很大的优势。

表 18.5：$M$ 个假设检验的一组可能的结果。可见其中假阳性检验的个数为 $V$；第一类错误率为 $\operatorname{E}(V)/M_0$。第二类错误率为 $\operatorname{E}(T)/M_1$，检验功效（power）为 $1-\operatorname{E}(T)/M_1$。

表 18.5 汇总了 $M$ 个假设检验的一组可能的结果。可见这里的 FWER 是 $\operatorname{Pr}(V\geq1)$。不过我们关注的却是 假发现率（false discovery rate）

在微阵列数据中，这就是被错误地认定显著的基因在所有被认定显著的基因中的期望比例。这个期望是取在生成样本数据的总体样本上的。Benjamini and Hochberg (1995) 首先提出了假发现率的概念，然后给出了一个检验过程（算法 18.2），其 FDR 不超过一个自定义的置信区间 $\alpha$。BH（Benjamini-Hochberg）方法是建立在 p 值上的；这些 p 值可以通过对测试统计量的渐进近似（例如高斯分布）、或者置换分布来计算，这里使用了置换分布。

算法 18.2：BH 方法

1. 给定一个假发现率 $\alpha$，令 $p_{(1)}\leq p_{(2)}\leq\dots\leq p_{(M)}$ 代表排序后的 p 值。
2. 定义： $$L = \max \left\{ j: p_{(j)} < \alpha \cdot \frac{j}{M} \right\} \tag{18.44}$$
3. 则 BH 拒绝阈值为 $p_{(L)}$，拒绝所有满足 $p_j\leq p_{(L)}$ 的零假设 $H_{0j}$。

如果这些假设都是独立的，则 Benjamini and Hochberg (1995) 证明了无论有多少零假设是正确的以及无论当零假设不成立时 p 值的分布是什么，这个检验方法都满足：

作为演示，我们使用了 $\alpha=0.15$。图 18.19 展示了一个排序后 p 值 $p_{(j)}$ 的散点图和斜率为 $0.15/12625$ 的直线。

BH 方法自左向右地找出 p 值落在直线下方的最后一次。这时的 $j=11$，所以拒绝了 p 值最小的 11 个基因的原假设。需要注意截断点处的第 11 个最小的 p 值是 0.00012，而对应的第 11 个最大的 $|t_j|$ 取值是 4.101。所以我们拒绝原假设的条件是 $|t_j|\geq 4.101$。

从上面简要的介绍中，并没有解释为什么 BH 方法有效；或者说，为什么结果的 FDR 不会大于所使用 $\alpha$ 的值 0.15。实际上，这个结论的证明非常复杂（Benjamini and Hochberg, 1995）。

一个更直接的处理方式是 插入（plug-in） 估计方法。 它并不是从一个选定的 $\alpha$ 出发的，而是先选定一个 t 统计量的截断点，例如上文中的值 $4.101$。观测样本中 $|t_j|$ 取值大于等于 $4.101$ 的有 11 个。在所有的置换中 $|t_j^k|$ 大于等于 $4.101$ 的一共有 1518 个，也就是每个置换平均有 $1518/1000=1.518$ 个。 所以一个假发现率的直接估计就是 $\widehat{\text{FDR}}=1.518/11\approx14%$。而 $14%$ 就大约等于上面使用到的 $\alpha=0.15$（两者的差别来自于采样的离散性）。算法 18.3 概括了这个方法。简单地概括一下：

当使用了式 18.40 中的置换 p 值时，算法 18.3 中 FDR 的插入估计与算法 18.2 中的 BH 方法是等价的。

算法 18.3：FDR 的插入估计

1. 创建数据的 $K$ 个置换样本，对每个特征变量 $j=1,2,\dots,M$ 和每个置换 $k=1,2,\dots,K$ 计算 t 统计量 $t_j^k$。
2. 对一个区间内很多的截断点取值 $C$，定义： $$\begin{gather} R_\text{obs} = \sum_{j=1}^M I(|t_j| > C) \\ \widehat{\operatorname{E}(V)} = \frac{1}{K} \sum_{j=1}^M \sum_{k=1}^K I(|t_j^k|>C) \end{gather}\tag{18.46}$$
3. （对每个 $C$）根据公式 $\widehat{\text{FDR}}=\widehat{\operatorname{E}(V)}/R_\text{obs}$ 来估计 FDR。

这个 BH 方法与插入估计之间的对应关系并不是一个偶然。练习 18.17 展示了它们一般地是等价的。请注意这个方法完全没有提起到 p 值，而是直接基于检验统计量来判断。

插入估计是基于下式中的近似：

一般来说 $\widehat{FDR}$ 是 FDR 的一个一致估计（Storey, 2002、Storey et al., 2004）。由于置换分布中使用了 $M$ 而不是 $M_0$ 个零假设，所以分子 $\widehat{\operatorname{E}(V)}$ 实际上估计的是 $(M/M_0)\operatorname{E}(V)$。所以如果可获得一个对 $M_0$ 的估计，那么 $(\hat{M}_0/M)\cdot\widehat{\text{FDR}}$ 可能是对 FDR 的一个更好的估计。同样地根据不等式 18.45，一个对 $M_0$ 的估计也可被用于改进 BH 方法。

读者可能会好奇我们为什么选择了 $0.15$ 作为 FDR 的上界 $\alpha$ 取值。需要提醒一下 FDR 与第一类错误是不同的，后者习惯性的选择是 0.05。对于科学研究者来说，假发现率是统计学家认定为显著基因的列表中假阳性（显著）基因的预期比例。取值 0.15 这么高的 FDR 的微阵列试验可能也是可用的，特别是当事实上这些试验是探索性的。

18.7.2 不对称截断点和 SAM 方法

在上述的检验方法中，我们使用的是检验统计量 $t_j$ 的绝对值，所以对统计量的正负取值使用的是同样的切割点。在一些试验中，大多或所有的区分性表达的基因可能都是在正方向上（或都是在负方向上）变化。在这种情况下，在两个方向上使用各自的截断点可能是有利的。

微阵列显著性分析（significance analysis of microarrays，SAM） 方法就是这样的一个工具。图 18.20 展示了 SAM 方法的基本原理。纵坐标上是排序后的检验统计量 $t_{(1)}\leq t_{(2)}\leq\cdots\leq t_{(M)}$，而横坐标是从对数据的置换中计算出的期望排序后统计量：$\tilde{t}_{(j)}=(1/K)\sum_{k=1}^Kt_{(j)}^k$，其中的 $t_{(1)}^k\leq t_{(2)}^k\leq\cdots\leq t_{(M)}^k$ 是置换 $k$ 的排序后的检验统计量。

图中画出了两条直线，与 45 度线平行并与之有 $\Delta$ 单位的距离。从原点出发向右侧移动，找出基因第一次冲破带状区域的位置。这就定义出了上截断点 $C_\text{hi}$，并且所有在这个点之上的基因都被认定为显著（标记为红色）。相似地，从原点向左移动，在左下角中寻找基因的下截断点 $C_\text{low}$。所以每一个调节参数 $\Delta$ 的取值都会定义上下截断点，在每个截断点的插入估计 $\widehat{\text{FDR}}$ 都可和之前一样计算。通常会在一个范围内选择多个 $\Delta$ 取值并计算对应的 $\widehat{\text{FDR}}$ 估计值，然后基于主观考虑从中选取特定的一组。

SAM 方法的优势在于其截断点可能是不对称的。在图 18.20 的例子中，当 $\Delta=0.71$ 时，得到了 11 个显著基因；而它们都位于右上角。左下角的数据点从未离开过带状区域，所以 $C_\text{low}=-\infty$。因此对这个 $\Delta$ 取值来说，在左侧（负方向）没有被认定为显著的基因。与第 18.7.1 中不同，这里对截断点没有对称性的假设，因为没有理由假设在两个方向上都有相似的性质。

这个方法与似然比检验中的非对称性 possible 有一些相似性。假设在无作用的零假设下有一个对数似然函数 $\ell_0(t_j)$，以及备择假设下一个对数似然函数 $\ell(t_j)$。那么假如对某个 $\Delta$ 下式成立，则似然比检验会拒绝零假设。

取决于似然函数、特别是它们的相对取值，这个判定条件对 $t_j$ 和对 $-t_j$ 可能会得出不同的阈值。SAM 方法拒绝零假设的条件是：

类似地，每个 $t_{(j)}$ 的阈值会取决于对应的零假设取值 $\tilde{t}_{(j)}$。

18.7.3 A Bayesian Interpretation of the FDR 😱

Storey (2002) 和 Efron and Tibshirani (2002) 从贝叶斯的角度为理解 FDR 提供了一个有趣的视角。首先将 阳性假发现率（positive false discovery rate，pFDR） 定义为

这个添加的“阳性（positive）”指的是我们只关注对存在阳性结果的情况下的错误率进行估计。这个对 FDR 微小的改动使其可以从贝叶斯的角度来理解。需要注意当 $\operatorname{Pr}(R=0)>0$ 时，常规的 FDR（式 18.43）是无法定义的。

令 $\Gamma$ 代表一个检验的拒绝域（rejection region）；在上面的例子中，拒绝域是 $\Gamma=(-\infty,-4.10)\cup(4.10,\infty)$。假设在独立同分布的统计量 $t_j,\dots,t_M$ 和拒绝域 $\Gamma$ 上进行 $M$ 次同样的简单假设检验。定义一个随机变量 $Z_j$，当零假设为真时等于 1,否则等于 0。假设每个 $(t_j,Z_j)$ 是独立同分布的随机变量对：

其中的 $F_0$ 和 $F_1$ 为某两个概率分布。这也就是说每个检验统计量 $t_j$ 来自两个概率分布其中之一：如果零假设为真则是 $F_0$，否则是 $F_1$。令 $\operatorname{Pr}(Z_j=0)=\pi_0$，则边际分布为：

那么可证明（Efron et al., 2001; Storey, 2002）

因此在混合模型 18.51 下，pFDR 就是给定了检验统计量落在检验的拒绝域中条件下的零假设为真的后验概率。也就是给定了拒绝零假设的条件下，零假设为真的后验概率（练习 18.20）。

假发现率提供了对基于一个全部拒绝域，比如 $|t_j|\geq2$，的检验的准确性的一个度量。但假如这样一个检验的 FDR 是（例如）10%，那么一个（例如）$t_j=5$ 的基因会比一个 $t_j=2$ 的基因更显著。因此我们想推导出一个 FDR 的局部（具体基因）的版本。一个检验统计量 $t_j$ 的 q 值（q-value）（Storey, 2003）被定义为所有会拒绝 $t_j$ 的拒绝域中最小的 FDR。也就是说，在对称的拒绝域场景中，$t_j=2$ 的 q 值被定义为拒绝域 $\Gamma=(-\infty,-2)\cup(2,\infty)$ 的 FDR。因此 $t_j=5$ 的 q 值会比 $t_j=2$ 的小，体现了 $t_j=5$ 比 $t_j=2$ 更显著的事实。在 $t=t_0$ 处的 局部假发现率（local false discovery rate） 被定义为（Efron and Tibshirani, 2002）：

这就是在 $t_j=t_0$ 附近的一个无穷小拒绝域上的（阳性）假发现率。

本节练习

练习 18.16

多重比较的邦费罗尼方法

假设在一个多重检验问题中，零假设为 $H_{0j}$，$j=1,2,\dots,M$；对应的 p 值为 $p_j$，$J=1,2,\dots,M$。令 $A$ 为至少一个零假设被错误地拒绝的事件，令 $A_j$ 为零假设 $j$ 被错误地拒绝的事件。假设我们使用了邦费罗尼方法，即当 $p_j<\alpha/M$ 时拒绝零假设 $j$。

1. 证明 $\operatorname{Pr}(A)\leq\alpha$。
   【提示：$\operatorname{Pr}(A_j\cup A_{j’})=\operatorname{Pr}(A_j)+\operatorname{Pr}(A_{j’})+\operatorname{Pr}(A_j\cap A_{j’})$】
2. 如果零假设 $H_{0j}$，$j=1,2,\dots,M$ 是独立的，那么 $$\operatorname{Pr}(A) = 1 - \operatorname{Pr}(A^c) = 1 - \prod_{j=1}^M \operatorname(A_j^C) = 1 - (1 - \alpha/M)^M$$ 利用这个等式，证明在这个问题中 $\operatorname{Pr}(A)\approx\alpha$。

练习 18.17

BH 方法和插入方法的等价关系

练习 18.18

Use result (18.53) to show that

(Storey, 2003).

练习 18.19

Consider the data in Table 18.4 of Section (18.7), available from the book website.

1. Using a symmetric two-sided rejection region based on the t-statistic, compute the plug-in estimate of the FDR for various values of the cut-point.
2. Carry out the BH procedure for various FDR levels α and show the equivalence of your results, with those from part (a).
3. Let (q.25, q.75) be the quartiles of the t-statistics from the permuted datasets. Let ˆπ0 = {#tj ∈ (q.25, q.75)}/(.5M), and set ˆπ0 = min(ˆπ0, 1). Multiply the FDR estimates from (a) by ˆπ0 and examine the results.
4. Give a motivation for the estimate in part (c).

(Storey, 2003)

练习 18.20

Proof of result (18.53).

Write

Use the fact that given R = k, V is a binomial random variable, with k trials and probability of success $\operatorname{Pr}(H=0|T\in\Gamma$, to complete the proof.