多重检验技术在大数据分析中的应用

Application of Multiple Test Techniques in Big Data Analysis

摘要: 在对大数据进行假设检验时，为了控制假阳性，需要采用多重检验技术。多重检验技术有多种，本文通过对大数据进行实际分析，比较各种算法的优缺点，给出不同方法的适用场合，从而对数据分析人员给以理论上的指导。文章首先阐述多重检验的必要性以及多重检验的相关概念；其次分别介绍对总体错误率和错误发现率进行控制的两类方法；最后将这几种多重检验方法应用到基因大数据中对基因的表达与否进行判断。实验结果表明，控制错误发现率的方法优于控制总体错误率的方法，在控制错误发现率的方法中，q值法的结果最好。原因在于q值法考虑了原假设的先验信息，能很好地控制错误发现率的大小，因此具有较高的精确性和检验功效。

Abstract: In the hypothesis test of big data, in order to control false positives, multiple test technology needs to be used. There are many kinds of multiple test techniques. This paper makes a practical analysis of big data, compares the advantages and disadvantages of various algorithms, and gives the application occasions of different methods, so as to give theoretical guidance to data analysts. Firstly, this paper expounds the necessity and the related concepts of multiple testing; Secondly, two kinds of methods to control the family-wise error rate and false discovery rate are introduced respectively; Finally, these multiple test methods are applied to gene big data to judge whether the genes are expressed or not. The experimental results show that the method of controlling the false discovery rate is better than the method of controlling the family-wise error rate. Among the methods of controlling the false discovery rate, the q-value method has the best result. The reason is that the q-value method considers the prior information of the original hypothesis and can well control the false discovery rate, so it has high accuracy and power.

1. 引言

在当今大数据时代，对于数据的搜集和获得不再是难题，难题出现在对数据的分析上。大数据的出现不仅仅是数据量的增多，它也使得数据分析方法发生了根本的变化，传统的统计分析方法不再适用，统计从业人员需要谨慎的分析数据，并利用合适的方法去分析，以免出现决策性的失误。多重假设检验就是相应于经典的单个假设检验提出的一种适应于大数据的假设检验方法。例如，对基因的选择问题，需要根据基因表达数据对大量的基因进行检测从而筛选出需要的基因。这就需要对大量基因进行检验，即多重检验问题。

在Lobenhofer et al. (2003) [1] 中，作者利用ORIOGEN算法 [2] 对表达的基因进行了挑选。本文利用多重检验技术，对Lobenhofer et al. (2003)中的基因重新进行了挑选，结果显示，在控制错误率的前提下，我们提高了检验的功效。

2. 多重假设检验概述

在N-P假设检验中，单个假设检验问题要求犯第一类错误的概率在可接受的范围内时，使犯第二类错误的概率控制在最小。在同时对多个假设进行检验时，对每个单独的假设检验而言第一类错误都在控制范围内，但对于整体而言犯第一类错误的概率将随着检验次数的增多而增大 [3]，超出可控范围。例如，取显著性水平 $\alpha =0.05$，当进行两次独立检验时，对总体而言犯第一类错误的概率为 $1-{\left(1-\alpha \right)}^{\text{2}}$ ；进行m次检验时总体犯第一类错误的概率为 $1-{\left(1-\alpha \right)}^m$，随着检验次数的增加总体错误率增大，如图1，当假设检验的次数超过100时，总体错误率接近于1。因此，在多重假设检验问题中，不能像单个假设检验一样控制第一类错误。

多重假设检验方法有很多种，它们都需要控制总体错误率(FWER)或错误发现率(FDR) [4]，才能使得第一类错误整体控制在 $\alpha$ 水平内。考虑m个假设检验 $H_i,i=1,2,\cdots ,m$，当原假设 $H_0$ 为真时，记为 $H_i=0$，否则 $H_i=1$。设 $m_0$ 和 $m_1$ 分别表示m个假设检验中 $H_0$ 和 $H_1$ 为真的个数。对m个假设检验结果的分类如表1所示。

Figure 1. Relationship between multiple tests on times m and type I error

图1. 多重检验次数m与I型错误的关系

Table 1. Possible results of m multiple tests

表1. m次多重检验可能出现的结果

总体错误率表示在m次检验中，至少出现一次错误拒绝 $H_0$ 的概率，即 $\text{FWER}=\text{Prob}\left(V\ge 1\right)$ 且 $\text{FWER}\le \alpha$。错误发现率表示在R次拒绝 $H_0$ 的检验中，错误拒绝所占比例的数学期望，即 $\text{FDR}=E\left(V/R\right)$，当 $R=0$ 时， $\text{FDR}=0$ 且 $\text{FDR}\le \alpha$。针对不同问题，选取的错误率控制指标也不同。在经典的多重检验中，通常m的取值较小，一般采取控制总体错误率的方法；但是在现在的大数据分析中，m的取值较大，此时采用控制总体错误率的话就过于严格，一般选择控制错误发现率的方法。

3. 多重假设检验方法

多重假设检验方法的具体算法可以分为两大类，一类为控制总体错误率的方法，另一类为控制错误发现率的方法，下面分别介绍这两类方法。

3.1. 控制总体错误率的算法

[Bonferroin算法] [5]

给定显著性水平 $\alpha$，对m个假设进行检验。采用单步法进行算法流程，每个假设各自的显著性水平选为 $\alpha /m$，如果第j个假设检验的p值 $p_j\le \alpha /m$，则拒绝 $H_{j0}$， $j=1,2,\cdots ,m$。因此调整过的p值为 ${\tilde{p}}_j=\min \left(m{p}_j,1\right)$。

[Homl算法] [6]

Homl算法与Bonferroin算法相比保守降低并且提高了功效。首先将m个假设检验的p值从小到大排序 $p_{\left(1\right)}\le p_{\left(2\right)}\le \cdots \le p_{\left(m-1\right)}\le p_{\left(m\right)}$，对应的检验为 $\left\{H_{0\left(1\right)},\cdots ,H_{0\left(m\right)}\right\}$。如果存在

$\hat{j}=\min \left\{j:p_{\left(j\right)}>\frac{\alpha }{m-j+1}\right\}$

则拒绝 $H_{0\left(j\right)},j=1,\cdots ,\hat{j}$。如果不存在 $\hat{j}$，则拒绝所有原假设。Homl过程调整过的p值为：

${\tilde{p}}_{\left(i\right)}=\underset{k=1,\cdots ,i}{\max }\left\{\min \left(\left(m-k+1\right)p_{\left(k\right)},1\right)\right\}$。

Homl算法为Bonferroin算法的改进，其他改进方法还包括Hommel、Hochberg算法 [7] [8]。针对大数据问题，我们关注的是在能够允许R次拒绝中发生少量的错误识别时，尽可能多地检验出显著的个体，即在控制错误发现率的同时，尽可能地提高检验的检验功效。

3.2. 控制错误发现率的算法

[BH算法] [9]

在显著性水平 $\alpha$ 下，控制FDR的过程如下：将原来的m个p值从小到大进行排序，即 $p_{\left(1\right)}\le p_{\left(2\right)}\le \cdots \le p_{\left(m-1\right)}\le p_{\left(m\right)}$。如果存在

$\hat{k}=\max \left\{j:p_{\left(j\right)}\le \frac{j}{m}\alpha \right\}$

则拒绝 $H_{0\left(j\right)},j=1,\cdots ,\hat{k}$ ；如果不存在 $\hat{k}$，则不拒绝任何原假设。

BH算法控制FDR时，FDR满足关系： $\text{FDR}\le \frac{m_0}{m}\alpha \le \alpha$。

[BY算法] [9]

根据BH算法，修改FDR的上界为 $m_0\alpha /\left(m{\displaystyle \underset{i=1}{\overset{m}{\sum }}\frac{1}{i}}\right)$，将原来的m个p值从小到大进行排序，即 $p_{\left(1\right)}\le p_{\left(2\right)}\le \cdots \le p_{\left(m-1\right)}\le p_{\left(m\right)}$。如果存在

$\hat{k}=\max \left\{j:m_0\alpha \le m{\displaystyle \underset{i=1}{\overset{m}{\sum }}\frac{1}{i}}\right\}$

则拒绝 $H_{0\left(j\right)},j=1,\cdots ,\hat{k}$。如果不存在 $\hat{k}$，则不拒绝任何原假设。

相对于BH算法而言，BY算法得出的上界值过于保守，所以j值减小，拒绝原假设的个数减少。在表1中表示为假设R值一定时，V值减小，S值增大。BH算法和BY算法的基本思想是给定显著性水平 $\alpha$，找到拒绝域，从而将错误水平控制在 $\alpha$ 以下。Storey (2002)提出一种控制FDR的直接方法，其基本思想是先凭借经验给出拒绝域，再估计错误率。

[q值法]

Storey (2002) [10] [11] 将 $E\left(V/R|R>0\right)$ 称为阳性错误拒绝率(pFDR)。设 ${\Gamma }_{\alpha }$ 为事先给定的拒绝域，则pFDR可表示为：

$\text{pFDR}=E\left[V\left({\Gamma }_{\alpha }\right)/R\left({\Gamma }_{\alpha }\right)|R>0\right]$

其中 $V\left({\Gamma }_{\alpha }\right)=\text{\#}\left\{H_i=0:T_i={\Gamma }_{\alpha }\right\}$ 表示错误发现次数， $R\left({\Gamma }_{\alpha }\right)=\text{\#}\left\{T_i\in {\Gamma }_{\alpha }\right\}$ 表示所有拒绝 $H_0$ 的次数。

定理1 [12] 对m个完全相同的假设进行检验，检验统计量为 $T_1,T_2,\cdots ,T_m$，显著区域为 $\Gamma$。假设 $\left(T_i,H_i\right)$ 是独立同分布的随机变量， $T_i|H_i~\left(1-H_i\right)\cdot F_0+H_i\cdot F_1$，其中 $F_0$ 为统计量在原假设的分布， $F_{\text{1}}$ 为统计量在备择假设的分布， $H_i~\text{Bernoulli}\left({\pi }_1\right)$， $i=1,2,\cdots ,m$。则：

$\text{pFDR}=\text{Prob}\left(H_0=0|T\in {\Gamma }_{\alpha }\right)=\frac{{\pi }_0\text{Prob}\left(T\in {\Gamma }_{\alpha }|H_0=0\right)}{\text{Prob}\left(T\in {\Gamma }_{\alpha }\right)}$

$\text{Prob}\left(T\in {\Gamma }_{\alpha }\right)={\pi }_0\text{Prob}\left(T\in {\Gamma }_{\alpha }|H_0=0\right)+\left(1-{\pi }_0\right)\text{Prob}\left(T\in {\Gamma }_{\alpha }|H_0=1\right)$

其中 ${\pi }_0=\text{1}-{\pi }_{\text{1}}=\text{Prob}\left(H_0=0\right)$ 为 $H_0$ 的先验概率。 $\text{pFDR}=\text{Prob}\left(H_0=0|T\in {\Gamma }_{\alpha }\right)$ 反映了在已经拒绝 $H_0$ 的条件下，该假设为真的概率，可认为pFDR是贝叶斯后验p值。

Storey对q值的定义如下：

$q\left(t\right)=\text{pFDR}\left({\Gamma }_{\alpha }\right)=\underset{{\Gamma }_{\alpha }:T\in {\Gamma }_{\alpha }}{\inf }\text{Prob}\left(H_0=0|T\in {\Gamma }_{\alpha }\right)$

即：统计量落入拒绝域时原假设为真的最小概率。q值不过是p值定义的一个逆过程，q值是在 $T\in {\Gamma }_{\alpha }$ 的条件下 $H_0=0$ 的概率。q值法与BH算法恰好相反，即通过选定拒绝域 ${\Gamma }_{\alpha }$ 去估计对应的q值，当 $q\le \alpha$ 时可保证 $\text{FDR}\le \alpha$。在文献 [13] 中Storey给出了 ${\pi }_0$ 和q值的具体估计算法。

4. 在基因表达数据分析中的应用

在当今的大数据时代，处处需要用到多重检验技术。下面我们以基因大数据为例，介绍几种多重检验方法的具体应用并进行比较，详细列举了每种方法的优缺点，以求能对实际数据工作者以相应的指导。

4.1. 实例

Lobenhofer et al. (2003)的微阵列实验中，评价了17-β对于MCF-7胸癌细胞的影响。该实验分别在6个时间点同时对1900个基因进行观测，每个时间点上有8个观测值。我们需要根据这些观测数据，判断这1900个基因在观测时间内是否表达。

检验问题描述如下：

$H_0:\mu \in C_0,H_1:\mu \in {\displaystyle \underset{i=1}{\overset{h}{\cup }}{C}_i}$

其中 $C_0=\left\{\mu \in R^6:{\mu }_1={\mu }_2=\cdots ={\mu }_6\right\}$， $C_{\text{1}}=\left\{\mu \in R^{\text{6}}:{\mu }_1\ge \cdots \ge {\mu }_i\ge \cdots \ge {\mu }_{\text{6}}\right\}$，， $C_i=\left\{\mu \in R^6:{\mu }_1\le \cdots \le {\mu }_i\ge \cdots \ge {\mu }_6\right\}$， $i=2,3,4,5$， $C_i=\left\{\mu \in R^{\text{6}}:{\mu }_1\ge \cdots \ge {\mu }_i\le \cdots \le {\mu }_{\text{6}}\right\}$， $i=\text{6},\text{7},\text{8},\text{9}$

Peddada et al. (2003)提出了ORIOGEN算法对该实验中的表达基因进行了选择，在该算法中为了降低假阳性，把假设检验中的p值调整为 ${\tilde{p}}_j=h\cdot p_j$， $j=1,2,\cdots ,m$，h为曲线模式的个数。

本文在ORIOGEN算法的基础上，利用多重假设检验技术对基因表达数据进行了分析。具体算法如下：

下面假定对某个选定的基因g进行分析，

第1步：选取所关注的表达曲线模式，将这些曲线形式记为 $C_1,C_2,\cdots ,C_h$。

第2步：利用PAVA算法 [14] 求出该基因在每个模式 $C_i,i=1,2,\cdots ,h$ 下的均值的估计值。

第3步：在每个曲线模式 $C_i,i=1,2,\cdots ,h$ 下，分别计算 $l_{\infty }$ 模。找到r满足 $l_{\infty }^g=l_{\infty }^{g\left(r\right)}={\max }_i{l}_{\infty }^{g\left(i\right)}$。(某个模式的 $l_{\infty }$ 模为其中任意两个参数估计值之差的最大值，参见Peddada et al. (2003)。)

第4步：对该基因根据其观测表达数据抽取N个bootstrap样本。对每个抽取的样本进行步骤2和3，从而获得该基因的统计量 $l_{\infty }^g$ 在原假设下 $H_0$ 的分布，从而根据样本统计量 $l_{\infty }^{g\left(r\right)}$ 的值求出该基因的 $p_g$ 值。

对所有的基因重复以上步骤，得到1900个p值。

第5步：进行多重假设检验：

$H_0:\mu \in C_0,H_1:\mu \in {\displaystyle \underset{i=1}{\overset{h}{\cup }}{C}_0}$

给定显著性水平 $\alpha$，对所有基因进行多重检验，挑选出表达显著的基因。

4.2. 实验结果

本文利用R统计软件，对1900个基因表达数据进行了分析，结果见下表。

Table 2. Analysis results of gene expression data controlled by FWER and FDR

表2. FWER控制和FDR控制的基因表达数据分析结果

根据控制错误率的不同，我们在两类不同的多重假设检验方法下分别进行了分析。表2报告了在 $N=1000000$，显著性水平 $\alpha =0.05$ 的情况下，利用上述多种算法控制FWER和FDR，得到的1900个基因中表达的基因个数。在原始的1900个 $p_i$ 值中，有423个 $p_i$ 值小于0.05。利用ORIOGEN算法识别出170个表达的基因 [15]。对于控制总体错误率，利用多重假设检验算法进行分析，Bonferroin算法识别出了17个表达的基因，Holm、Hochberg、Hommel算法也都识别出了17个表达基因；对于控制错误发现率，BH算法识别出151个表达的基因，由于BY算法的上界值过于保守，仅识别出31个表达的基因，q值法识别出199个表达的基因。显然，在假阳性水平相同的情况下，控制FDR的多重检验算法挑选出的表达基因个数远高于控制FWER的算法，其中q值法挑选出的表达基因最多，比ORIOGEN算法挑出的基因还要多，而且q值法的检验功效比较大。因此，在实际数据分析中，推荐使用q值法进行多重假设检验分析。

5. 总结

在大数据的假设检验分析中，需要使用多重检验技术来控制错误率。从实例数据分析中可以看出，在使用多重检验算法时，控制FWER的意义并不大。研究者更关心的问题是当错误识别个数控制在可以接受的范围内时，尽可能多地识别出显著的基因。因此推荐控制FDR的多重检验算法。在控制FDR的算法中，从本文实例分析的结果可以看出q值法比BH算法的检验功效更大，因为其考虑了先验信息。在当前的大数据时代，数据量的变化也对传统的统计理论提出了挑战。在检验问题中，利用q值法解决大数据的多重检验问题，具有很强的实际意义。

基金项目

国家自然科学基金项目11401393。

辽宁省教育厅自然科学基金项目LJC201914。

NOTES

^*通讯作者。