关于基因多效性检验中互斥假设检验的总结与探索

Summary and Exploration of Exclusive Hypothesis Test of Pleiotropy

摘要: 基因多效性是指单个基因能够影响多个性状的表达水平。能够有效检测基因多效性的存在并了解其产生的原因具有重要的统计学以及生物学意义。复杂原假设条件下的互斥检验方法十分适用于该类问题且具有较好的说服力。不失一般性的情况下，本文主要针对于两个性状的汇总统计量进行研究。我们将新构建的方法与IUT (intersection-union test)、PLACO (pleiotropic analysis under composite null hypothesis)、DACT (divide-aggregate composite-null test)、MAIUP (mixture-adjusted inter-section-union test for pleiotropy test)几种方法进行了对比。模拟结果显示，新方法与MAIUP方法相比在能够在某些特定的情形下能够更好的控制一类错误同时保证检验功效。最后将上述方法应用于冠状动脉疾病以及扁桃体切除的数据来进行实证分析。

Abstract: Pleiotropy means single gene can influence multiple traits. Detection and understanding of pleiot-ropy are statistically and biologically important. Mutually exclusive test under composite null hy-pothesis is a suitable and convincing approach. Without loss of generality, we mainly focus on ana-lyzing summary statistics of two traits in this study. We compared our methods with IUT, PLACO, DACT, and MAIUP through simulations. Simulations show that the type I error can be controlled better by proposed method compared with MAIUP under specific situations with maintaining the statistical power simultaneously.

1. 引言

在遗传学发展的开始阶段，人们对于基因的了解主要借助于对性状的研究，而在其中也产生出单个基因影响单个性状的最初概念。随着时间的推移，人们逐渐认识到基因不仅能直接控制某一性状的表现，而且通过某一性状的改变，间接能够影响到许多性状。比如，翻毛鸡的翻毛基因，不仅影响鸡羽的反卷，而且还影响其体温下降，心跳增快，心、脾扩大，生育力降低等，这说明基因通过生理生化过程而影响一系列性状的表现 [1] 。

探索基因多效性对于统计学和生物学都具有重要意义。从统计学角度来看，基因多效性的引入提高了关联分析的检验功效 [2] ，使得更加准确地发现的相关表型之间的遗传关系 [3] ，揭示无法通过分析单一表型来识别的遗传结构，提高表型预测的准确性 [4] 。此外对于生物学而言，有可能提高对全基因组学等基因医学的理解，或者是定位可影响多种性状或疾病的药物遗传靶点，从而推动高效药物的研发。

因此，大量关于研究基因多效性的统计学方法应运而生。从单个遗传变异位点与多个性状的单变量边际关联的比较，到所有性状对一个遗传变异位点同时进行回归的多变量分析，再到一个遗传变异位点对所有性状进行反向回归分析的方法，各类方法仍然在持续发展。其中针对单变量分析经常使用meta分析 [5] ，此方法虽然简单可行但忽略了性状之间的相关性；而针对多变量的分析中常用‘所有性状同时对单个变异的回归’、‘小等位基因剂量对所有性状的回归’以及‘性状与小等位基因剂量矩阵的典型相关’等方法 [6] 。

尽管对于多效性研究投入了大量的关注，但上述方法所构造的原假设并不能正确的反映基因多效性的含义。这是由于上述方法中的原假设是：没有一个性状与一个变异位点相关，但是仅有一个性状与变异有关时也可以拒绝这一原假设，此时并不符合多效性的定义，因此需要事后分析来解释是否存在多效性。然而随着基因变异数量的增加，工作量也会不断加大。本文所研究的互斥检验的统计方法，能够很好的提高检验多效性方法的可解释性和效率，它将原假设拆分为一系列的子原假设，并且这些子原假设互斥。即当且仅当子原假设其一成立时，原假设才成立。

Schaid等人 [7] 基于两个性状情形下的交并检验思想出发，提出了一个基于线性模型和多元正态分布的顺序检验框架来分析多效性，并且与以往的假设检验原假设的构造不同，此时将原假设更改为：至多有一个性状与基因变异相关联。此后基于Schaid等人提出的序列检验框架汪毅等人 [8] 给其赋予了互斥检验的名称。互斥检验改变了之前原有的原假设的构成，从根本上做出了修正，此时拒绝原假设则能够表明存在多效性，因此给检验生物基因多效性提供了强有力的研究方向。

在Schaid等人提出的方法中，当检验没有一个性状与基因变异相关这一子假设时，存在第一类错误十分保守的问题。之后，汪毅等人参考Schaid等人的方法对 值的构造实施了改进。他们基于似然比以及BIC等相关知识给予原假设中拆分成的每个子假设不同的权重，由此得到原假设的加权 值，结果显示此时得到的 值保守程度减少，能够基本上达到一类错误率的合理范围。

但上述方法仍然存在局限性，其方法构建主要使用个体水平数据(即已知基因型、性状等的相关数据)，容易涉及样本人群隐私信息，所以该类方法无法得到广泛应用。

汇总统计量的出现有效解决了这一问题，其主要包含单个变异位点与单个性状之间关联性效应大小、回归估计误差以及检验 值等总括信息。在这一信息来源的支持下，又促成了诸多基于汇总统计量的基因多效性检验方法。在其中较为突出的有Debashree Ray等人 [9] 提出的PLACO (pleiotropic analysis under composite null hypothesis)、Zhonghua Liu等人 [10] 提出的DACT (divide-aggregate composite-null test )以及Ting Wang等人 [11] 提出的MAIUP (mixture-adjusted intersection-union test for pleiotropy test)法。本文主要针对这几种方法进行对比分析，并做出了一定的改进。

2. 符号定义

在不失一般性的前提下，本文主要研究各类方法针对于两个性状的情况。

假设收集的两个性状的相关数据来自样本数为 $n_1$ ， $n_2$ 的个体，设 $n_1=n_2=n$ ，假设个体 $n_1$ 与 $n_2$ 之间独立，即研究之间没有重叠样本。 $Y_k\left(k=1,2\right)$ 代表第k个性状的值，X表示基因型。其中两种性状的值表示为： $Y_1=\left(y_{11},y_{12},\cdots ,y_{1n}\right)$ ， $Y_2=\left(y_{21},y_{22},\cdots ,y_{2n}\right)$ ， $y_{ij}\left(i=1,2;j=1,2,\cdots ,n\right)$ 表示来自第 个个体的第 性状的值；位点的基因型值表示为： $X=\left(x_1,x_2,\cdots ,x_n\right)$ ， $x_i\left(i=1,2,\cdots ,n\right)$ 表示基因值。若所考虑的遗传变异是双等位单核苷酸多态性(SNP)，个体的基因型可以取值为0、1或2。为简单起见，我们假设不存在协变量，通过考虑性状残差(通过对性状上的协变量进行回归得到)而不是原始的性状值，可以很容易地放宽松这一假设条件。尽管残差结果数据不是标准化的数据，但Broadaway KA [12] 等人的研究表明它不会影响基因关联测试的有效性。

对于给定基因变量，即有标准单变量回归模型表示为：

$y_1=X{\beta }_1+{\epsilon }_1$ , $y_2=X{\beta }_2+{\epsilon }_2$ , (2.1)

其中， ${\beta }_k\left(k=1,2\right)$ 是第k个性状的遗传效应， ${\epsilon }_1~N\left(0,I_n{\sigma }_1^2\right)$ 、 ${\epsilon }_2~N\left(0,I_n{\sigma }_2^2\right)$ ， $I_n$ 表示为n维单位矩阵。

Wald检验统计量定义为： $Z_k=\frac{\widehat{{\beta }_k}}{\widehat{{\sigma }_k}}$ 。其中， $\widehat{{\sigma }_k}=\widehat{se}\left(\widehat{{\beta }_k}\right)$ ， $\widehat{{\beta }_k}$ 是关于 ${\beta }_k$ 的极大似然估计， $\widehat{{\sigma }_k}$ 是其标准差的估计。

基于上述模型检验多效性的假设检验设定为：

$H_0$ ：至多只有一个 ${\beta }_k$ 不为零 vs $H_1$ ： ${\beta }_1\ne 0$ 且 ${\beta }_2\ne 0$ (2.2)

由互斥检验的定义，可以将原假设 $H_0$ 分为三个互斥的子假设 $H_{00}$ 、 $H_{01}$ 、 $H_{10}$ 组成。即：

$H_0=\{\begin{array}{ll}{H}_{00}:{\beta }_1={\beta }_2=0\hfill & ①\hfill \\ H_{01}:{\beta }_1=0,{\beta }_2\ne 0\hfill & ②\hfill \\ H_{10}:{\beta }_1\ne 0,{\beta }_2=0\hfill & ③\hfill \end{array}$ vs $H_1:{\beta }_1\ne 0且{\beta }_2\ne 0$ , (2.3)

即有 $H_0=H_{00}\cup H_{01}\cup H_{10}$ ，而 $H_1$ 代表我们最感兴趣的多效性的情况。

3. P值的构建

3.1. IUT法

上文公式(2.3)可以变换为 $H_0=H_0^{\left(1\right)}\cup H_0^{\left(2\right)}$ 。其中：

$H_0^{\left(k\right)}:{\beta }_k=0\text{vs}{H}_1^{\left(k\right)}:{\beta }_k\ne 0$ (3.1)

针对每个 $H_0^{\left(k\right)}$ 检验都会有一个对应的 $p_k$ 值，因此关于原假设 $H_0$ 的p值由交并检验知识可知为：

$p=\max \left\{p_1,p_2\right\}$ . (3.2)

此时 $p_1,p_2$ 值计算依赖于Wald检验统计量 $Z_k\left(k=1,2\right)$ 在 $H_0^{\left(k\right)}$ 下服从标准正态分布计算而得。因此IUT方法也称为最大p值法或联合显著性检验 [13] 。

然而此方法在任意有限样本量的情况下一类错误会出现保守现象，尤其是当 $H_{00}$ 成立时候，原因如下：

假设两性状独立，则有 $p_1$ 和 $p_2$ 是独立的，所以有

$\mathrm{Pr}\left(p<\alpha \right)=\mathrm{Pr}\left(p_1<\alpha \right)\mathrm{Pr}\left(p_2<\alpha \right)$ . (3.3)

在有限样本量下， $H_{01}$ 情形有 $\mathrm{Pr}\left(p_1<\alpha |{\beta }_1=0\right)=\alpha$ ，并且 $\mathrm{Pr}\left(p_2<\alpha |p_2\ne 0\right)<1$ ，所以有 $\mathrm{Pr}\left(p<\alpha \right)<\alpha$ ；同理在 $H_{10}$ 情形下亦是如此。

然而在 $H_{00}$ 会更保守的原因是：

由于 $\mathrm{Pr}\left(p_1<\alpha \right)=\mathrm{Pr}\left(p_2<\alpha \right)=\alpha$ ，即有 $\mathrm{Pr}\left(p<\alpha \right)={\alpha }^2$ ，对任意的 $\alpha \in \left(0,1\right)$ 。因此在此种情形下一类错误会更加保守。

3.2. PLACO法

为了能够改进IUT方法下一类错误的保守性，同时考虑多效性检验的复合假设的特性，Debashree Ray [9] 提出了PLACO，此方法的不同之处在于它没有依赖于p值，而依赖于Z检验统计量。通过观察假设(3.1)式中的原假设发现，原假设可以改写为：

$H_0:{\beta }_1{\beta }_2=0\text{vs}{H}_1:{\beta }_1{\beta }_2\ne 0$ (3.4)

因此需要考虑 $Z_1*Z_2$ 的分布，可以借助Craig [14] 提到的关于两变量 $a,b$ 独立且均服从于正态分布，则两变量分别标准化后得到的乘积X的密度函数为： $f\left(x\right)=\frac{{\kappa }_0\left(x\right)}{\pi }$ ，其中 ${\kappa }_0\left(x\right)$ 为第二类的阶数为零阶的修正贝塞尔函数。

换句话说，假设有如下情形： $H_{00}$ 成立时， $Z_1$ 和 $Z_2$ 独立且服从于 $N\left(0,1\right)$ ； $H_{01}$ 成立时， $Z_1$ 和 $Z_2$ 独立且分别服从于 $N\left(0,1\right)$ 和 $N\left(0,1+{\tau }_2^2\right)$ ； $H_{10}$ 成立时， $Z_1$ 和 $Z_2$ 独立且分别服从于 $N\left(0,1+{\tau }_2^2\right)$ 和 $N\left(0,1\right)$ 。假设在全局假设下 $H_{00}$ 成立的概率为 ${\pi }_{00}$ ， $H_{01}$ 以 ${\pi }_{01}$ 的概率成立， $H_{10}$ 以概率 ${\pi }_{10}$ 成立，则 $Z_2$ 服从标准正态分布 $N\left(0,1\right)$ ，而 $Z_1$ 在给定均值 ${\mu }_1$ 服从 $N\left(0,{\tau }_1^2\right)$ 前提下，服从 $N\left({\mu }_1,1\right)$ 的分布。在此假设下p值计算为：

$\begin{array}{c}{p}_{z_1{z}_2}=2*P_{H_0}\left(Z_1{Z}_2>\left|z_1{z}_2\right|\right)\\ =2*\left\{P\left(H_{00}\right)P_{H_{00}}\left(Z_1{Z}_2>\left|z_1{z}_2\right|\right)+P\left(H_{01}\right)P_{H_{01}}\left(Z_1{Z}_2>\left|z_1{z}_2\right|\right)+P\left(H_{10}\right)P_{H_{10}}\left(Z_1{Z}_2>\left|z_1{z}_2\right|\right)\right\}\\ ={\pi }_{00}F\left(z_1{z}_2\right)+{\pi }_{01}F\left(\frac{z_1{z}_2}{\sqrt{1+{\tau }_2^2}}\right)+{\pi }_{10}F\left(\frac{z_1{z}_2}{\sqrt{1+{\tau }_1^2}}\right)\end{array}$ (3.5)

其中 $z_1$ 和 $z_2$ 为观测值，且 $F\left(u\right)=2{\displaystyle {\int }_{\left|u\right|}^{\infty }f\left(x\right)\text{d}x}$ 是标准正态乘积分布的双边尾部概率。

因此，该方法主要需要知道均值 ${\mu }_1$ 和 ${\mu }_2$ 的数值，以及了解三个子原假设可能成立的概率。PLACO方法通过设置额外的模型假设(即假设两性状独立)来克服上述问题带来的困难，最终根据Huang YT [15] 中提到的方法中在不涉及估计未知参数的同时可以将p值近似的估计为：

${\hat{p}}_{z_1{z}_2}\approx F\left(z_1\frac{z_2}{\sqrt{\mathrm{var}\left(z_2\right)}}\right)+F\left(\frac{z_1}{\sqrt{\mathrm{var}\left(z_1\right)}}\right)-F\left(z_1{z}_2\right)$ . (3.6)

3.3. DACT法

与PLACO不同，DACT直接估计三个子原假设的比例，最终将单个子原假设获得的P值的加权总和作为结果。

假设取自m个样本，则需要进行m个假设检验来进行检验 $H_0:{\beta }_1{\beta }_2=0\text{vs}{H}_1:{\beta }_1{\beta }_2\ne 0$ 。在回归模型中针对于第一个性状所涉及到的参数 ${\beta }_1$ 检验有m个假设检验： $H_{0j}^{{\beta }_1}:{\beta }_1=0$ ，同样针对第二个性状所涉及到的参数 ${\beta }_2$ 检验亦有m个假设检验： $H_{0j}^{{\beta }_2}:{\beta }_2=0$ ，这里 $\left(1\le j\le m\right)$ 。同时将 $H_j^{{\beta }_1}$ (对于 $H_j^{{\beta }_2}$ 亦是如此)定义为伯努利随机变量，即：如果假设 $H_{0j}^{{\beta }_1}$ 为真，则 $H_j^{{\beta }_1}=0$ ；如果假设 $H_{0j}^{{\beta }_1}$ 为假，则 $H_j^{{\beta }_1}=1$ ， $\left(1\le j\le m\right)$ 。正如前文独立的假设条件，可以得到 $H_j^{{\beta }_1}$ 和 $H_j^{{\beta }_2}$ 是独立的。对于每一个基因，通过检验统计量可以获得检验 ${\beta }_1$ 的p值 $p_{{\beta }_{1}j}$ 以及检验 ${\beta }_2$ 的p值 $p_{{\beta }_{2}j}$ ， $\left(1\le j\le m\right)$ 。

由Efron [16] 假设 $P\left(H_j^{{\beta }_1}=0\right)={\pi }_0^{{\beta }_{{}_1}}$ ， $P\left(H_j^{{\beta }_2}=0\right)={\pi }_0^{{\beta }_2}$ ，其中，参数 ${\pi }_0^{{\beta }_1}$ 是第一个性状不受遗传基因影响在全基因组中的比例， ${\pi }_0^{{\beta }_2}$ 是第二个性状不受遗传基因影响在全基因组中的比例。即有：(2.8)式情况①有 $\mathrm{Pr}\left(H_j^{{\beta }_1}=0,H_j^{{\beta }_2}=0\right)={\pi }_0^{{\beta }_1}{\pi }_0^{{\beta }_2}$ ；情况②有 $\mathrm{Pr}\left(H_j^{{\beta }_1}=0,H_j^{{\beta }_2}=1\right)={\pi }_0^{{\beta }_1}\left(1-{\pi }_0^{{\beta }_2}\right)$ ；情况③有 $\mathrm{Pr}\left(H_j^{{\beta }_1}=1,H_j^{{\beta }_2}=0\right)=\left(1-{\pi }_0^{{\beta }_1}\right){\pi }_0^{{\beta }_2}$ ；备择假设有 $\mathrm{Pr}\left(H_j^{{\beta }_1}=1,H_j^{{\beta }_2}=1\right)=\left(1-{\pi }_0^{{\beta }_1}\right)\left(1-{\pi }_0^{{\beta }_2}\right)$ 。

由上述公式成立可知，多效性的假设检验(2.8)式中的三个子假设的相对比例则分别为： ${\omega }_1=\frac{{\pi }_0^{{\beta }_1}{\pi }_0^{{\beta }_2}}{c}$ ， ${\omega }_2=\frac{{\pi }_0^{{\beta }_1}\left(1-{\pi }_0^{{\beta }_2}\right)}{c}$ ， ${\omega }_3=\frac{\left(1-{\pi }_0^{{\beta }_1}\right){\pi }_0^{{\beta }_2}}{c}$ ， $c={\pi }_0^{{\beta }_1}{\pi }_0^{{\beta }_2}+{\pi }_0^{{\beta }_1}\left(1-{\pi }_0^{{\beta }_2}\right)+\left(1-{\pi }_0^{{\beta }_1}\right){\pi }_0^{{\beta }_2}$ 。

此时只需要需要估计 ${\pi }_0^{{\beta }_1},{\pi }_0^{{\beta }_2}$ 即可。通过Jin和Cai [17] 提出使用经验特征函数和傅里叶分析来估计真原假设比例的估计方法，并且是基于检验 ${\beta }_1=0,{\beta }_2=0$ 的Z统计量来进行估计。假设有m个Z检验统计量，即： $Z_j~N\left({\mu }_j,{\sigma }_j^2\right),\left(1\le j\le m\right)$ 。这里我们设置为 ${\mu }_j={\mu }_0=0,{\sigma }_j^2={\sigma }_0^2=1$ 。给定 $r\in \left(0,\frac{1}{2}\right)$ ，关于真原假设 的比例可以一致地估计为：

$\widehat{{\pi }_0}=\underset{\left\{0\le t\le \sqrt{2r\log \left(m\right)}\right\}}{\sup }\left\{{\displaystyle {\int }_{-1}^1\left(1-\left|\zeta \right|\right)\left(\mathrm{Re}\left({\psi }_m\left(t\zeta ,Z_1,\cdots ,Z_m\right)\exp \left(-i{\mu }_{0}t\zeta +{\sigma }_0^2{t}^2{\zeta }^2/2\right)\right)\right)\text{d}\zeta }\right\}$ (3.7)

此处 $\mathrm{Re}\left(a\right)$ 表示复数a的实部。其中 ${\psi }_m$ 为经验特征函数，为 ${\psi }_m\left(t\right)=\frac{1}{m}{\displaystyle \underset{j=1}{\overset{m}{\sum }}\exp \left(it{Z}_j\right)}$ 。因此通过此方法可以估计 ${\pi }_0^{{\beta }_1}$ 和 ${\pi }_0^{{\beta }_2}$ 以及相应的权重 ${\omega }_1$ 、 ${\omega }_2$ 、 ${\omega }_3$ 。

在有了权重后，还需要考虑如何定义三个子原假设中进行多效性的检验。则有在检验 ${\beta }_1\ne 0$ 且 ${\beta }_2=0$ 这种情形时只需要通过 $p_{{\beta }_2}$ 来检验 ${\beta }_2=0$ 即可，这是因为 ${\beta }_1\ne 0$ ；同样在检验 ${\beta }_1=0$ 且 ${\beta }_2\ne 0$ 时，由于 ${\beta }_2\ne 0$ ，所以只需要通过 $p_{{\beta }_1}$ 来检验 ${\beta }_1=0$ 即可；对于 ${\beta }_1={\beta }_2=0$ 这种情况，需要检验 ${\beta }_1$ 和 ${\beta }_2$ 是否都非零。在给定显著性水平 $\alpha$ 下，当 $MaxP=\max \left(p_{{\beta }_1},p_{{\beta }_2}\right)<\alpha$ 时，就会拒绝此种情况。给出的子假设 ${\beta }_1={\beta }_2=0$ 情况下中，由于 $MaxP$ 服从 $\beta \left(2,1\right)$ 分布，其平方服从均匀分布，因此在此情况下的p值取为 ${\left(MaxP\right)}^2$ 。

即有p值计算为：

$p=\{\begin{array}{l}{p}_{1j}={\left(Max{P}_j\right)}^2;\text{if}{\beta }_1={\beta }_2=0\\ p_{2j}=p_{{\beta }_{1j}};\text{if}{\beta }_1=0,{\beta }_2\ne 0\\ p_{3j}=p_{{\beta }_{2j}};\text{if}{\beta }_1\ne 0,{\beta }_2=0\end{array}$ (3.8)

因此原假设的p值计算为：

$p=\widehat{{\omega }_1}{p}_{1j}+\widehat{{\omega }_2}{p}_{2j}+\widehat{{\omega }_3}{p}_{3j}$ . (3.9)

3.4. MAIUP法

由于基于IUT的方法时，尤其是在检验多效性中的两个基因与某一性状都无关这一检验的时候，所获得的一类错误往往过于保守。为了修正IUT在探索多效性时一类错误的保守性，Ting Wang [11]提出了MAIUP的方法。MAIUP基于传统的交并检验，以两组独立的p值作为输入值，并基于最初在高维中介分析 [18] (mediation analysis)框架下而提出的一个新想法。MAIUP的主要修正在于它通过拟合三组子原假设的分布，将多效性检验的复合假设检验的形式的原假设考虑在内，从而最终生成有良好校准性的p值。简言之，MAIUP通过遵循IUT的原则来检测基因多效性，并且用一对表型的每个遗传位点的两个p值的最大值(用 $p_{\max }$ )表示作为显著性的度量。具体来说，从IUT方法中可以获得的一组 $p_{\max ,j}\left(j=1,2,\cdots ,S\right)$ ；(S表示样本个数)。此时，我们不依赖于其从0到1的常规的均匀分布，而是通过借用整个基因组中所有基因的有用信息来估计每个子原假设的比例，并通过拟合三个子原假设混合的分布来构建 $p_{\max ,j}$ 的新的分布。即：

$\begin{array}{l}\mathrm{Pr}\left(p_{\max ,j}\le t|H_{0j}\right)\\ =\mathrm{Pr}\left(p_{\max ,j}\le t|H_{01,j}\right)\mathrm{Pr}\left(H_{01,j}\right)+\mathrm{Pr}\left(p_{\max ,j}\le t|H_{10,j}\right)\mathrm{Pr}\left(H_{10,j}\right)\\ +\mathrm{Pr}\left(p_{\max ,j}\le t|H_{00,j}\right)\mathrm{Pr}\left(H_{00,j}\right)\\ ={\pi }_{01}\mathrm{Pr}\left(p_{1j}\le t|H_{01,j}\right)\mathrm{Pr}\left(p_{2j}\le t|H_{01,j}\right)+{\pi }_{10}\mathrm{Pr}\left(p_{1j}\le t|H_{10,j}\right)\mathrm{Pr}\left(p_{2j}\le t|H_{10,j}\right)\\ +{\pi }_{00}\mathrm{Pr}\left(p_{1j}\le t|H_{00,j}\right)\mathrm{Pr}\left(p_{2j}\le t|H_{00,j}\right)\\ ={\pi }_{01}t\mathrm{Pr}\left(p_{2j}\le t|H_{01,j}\right)+{\pi }_{10}t\mathrm{Pr}\left(p_{2j}\le t|H_{10,j}\right)+{\pi }_{00}{t}^2\end{array}$ (3.10)

上式成立前提是 $p_{1j}$ 和 $p_{2j}$ 两者独立，其中， ${\pi }_{00}$ 、 ${\pi }_{01}$ 和 ${\pi }_{10}$ 是相对应的三个子原假设所对应的比例，t是评价多效性显著性给定的一个阈值。同时， $\mathrm{Pr}\left(p_{1j}\le t|H_{10,j}\right)$ 是在假设 $H_{10,j}$ 下拒绝 ${\beta }_{1j}=0$ 的功效函数，同理 $\mathrm{Pr}\left(p_{2j}\le t|H_{01,j}\right)$ 是在假设 $H_{01,j}$ 下拒绝 ${\beta }_{2j}=0$ 的功效函数。当在样本量比较大时，两者趋于1，因此最终结果变为 ${\pi }_{01}t+{\pi }_{10}t+{\pi }_{00}{t}^2$ 。

在计算上式时还需要估计参数 ${\pi }_{01},{\pi }_{10},{\pi }_{00}$ 。依据James Y. Dai [18] 和John D. Storey [19] 提出的估计方法，可以进行上述参数的估计。以 ${\pi }_{0+}\left({\lambda }_1\right)$ 表示为 ${\beta }_1=0$ 时的假设下的原假设成立的比例， ${\pi }_{+0}\left({\lambda }_2\right)$ 表示为 ${\beta }_2=0$ 时的假设下的原假设成立的比例。因此 ${\pi }_{0+}\left({\lambda }_1\right)$ 和 ${\pi }_{+0}\left({\lambda }_2\right)$ 可以保守的估计为：

$\widehat{{\pi }_{0+}}\left({\lambda }_1\right)=\frac{1}{S\left(1-{\lambda }_1\right)}\#\left\{p_{1j}>{\lambda }_1\right\}$ ， $\widehat{{\pi }_{+0}}\left({\lambda }_2\right)=\frac{1}{S\left(1-{\lambda }_2\right)}\#\left\{p_{2j}>{\lambda }_2\right\}$ . (3.11)

其中S表示为在分析中样本的所有数目， ${\lambda }_1,{\lambda }_2$ 表示两个调谐参数，用于确定观察到的p值是否来自原假设情况。与此同时 ${\pi }_{00}$ 可以估计为：

$\widehat{{\pi }_{00}}\left({\lambda }_1,{\lambda }_2\right)=\frac{1}{S\left(1-{\lambda }_1\right)\left(1-{\lambda }_2\right)}\#\left\{p_{1j}>{\lambda }_1,p_{2j}>{\lambda }_2\right\}$ . (3.12)

所以 ${\pi }_{01},{\pi }_{10}$ 计算分别为：

$\widehat{{\pi }_{01}}\left({\lambda }_1\right)=\widehat{{\pi }_{0+}}\left({\lambda }_1\right)-\widehat{{\pi }_{00}}\left({\lambda }_1,{\lambda }_2\right)$ ; $\widehat{{\pi }_{10}}\left({\lambda }_1\right)=\widehat{{\pi }_{+0}}\left({\lambda }_1\right)-\widehat{{\pi }_{00}}\left({\lambda }_1,{\lambda }_2\right)$ . (3.13)

借助James Y. Dai [18] 试用的结果，此处采用 ${\lambda }_1={\lambda }_2=0.5$ 。将上述估计值带入式(3.10)中，即可获得最终的结果。

3.5. 新构造的P值

DACT方法中在估计三个子原假设的比例时借用了独立的严格前提条件，而在MAIUP方法中估计三个子原假设的比例时条件相对较为宽泛，所以在两种方法估计出来的比例进行比较后，发现MAIUP方法做得到的比例更为接近，因此考虑利用MAIUP方法中估计三个子原假设的比例作为DACT方法中子原假设的比例。即有 ${\pi }_{00},{\pi }_{01},{\pi }_{10}$ 估计为：

$\widehat{{\pi }_{00}}\left({\lambda }_1,{\lambda }_2\right)=\frac{1}{S\left(1-{\lambda }_1\right)\left(1-{\lambda }_2\right)}\#\left\{p_{1j}>{\lambda }_1,p_{2j}>{\lambda }_2\right\}$ ；

$\widehat{{\pi }_{01}}\left({\lambda }_1\right)=\widehat{{\pi }_{0+}}\left({\lambda }_1\right)-\widehat{{\pi }_{00}}\left({\lambda }_1,{\lambda }_2\right)$ ；

$\widehat{{\pi }_{10}}\left({\lambda }_1\right)=\widehat{{\pi }_{+0}}\left({\lambda }_1\right)-\widehat{{\pi }_{00}}\left({\lambda }_1,{\lambda }_2\right)$

其中， $\widehat{{\pi }_{0+}}\left({\lambda }_1\right)=\frac{1}{S\left(1-{\lambda }_1\right)}\#\left\{p_{1j}>{\lambda }_1\right\}$ ， $\widehat{{\pi }_{+0}}\left({\lambda }_2\right)=\frac{1}{S\left(1-{\lambda }_2\right)}\#\left\{p_{2j}>{\lambda }_2\right\}$ ，并且S表示为在分析中样本的所有数目， ${\lambda }_1$ 和 ${\lambda }_2$ 表示两个调谐参数，用于确定观察到的p值是否来自原假设情况。借助James Y. Dai [18] 比较试用的结果，此处采用 ${\lambda }_1={\lambda }_2=0.5$ 。即各个子原假设权重为 ${\omega }_1=\frac{\widehat{{\pi }_{00}}}{c}$ ， ${\omega }_2=\frac{\widehat{{\pi }_{01}}}{c}$ ， ${\omega }_3=\frac{\widehat{{\pi }_{10}}}{c}$ ， $c=\widehat{{\pi }_{00}}+\widehat{{\pi }_{01}}+\widehat{{\pi }_{10}}$ 。

对于提出的新方法里各个子假设情形下每种情况p值取值同DACT方法。即在检验 ${\beta }_1\ne 0,{\beta }_2=0$ 这种情形时只需要通过 $p_{{\beta }_2}$ 来检验 ${\beta }_2=0$ ，是因为 ${\beta }_1\ne 0$ ；同样在检验 ${\beta }_1=0,{\beta }_2\ne 0$ 时，由于 ${\beta }_2\ne 0$ ，所以只需要通过p值 $p_{{\beta }_1}$ 来检验 ${\beta }_1=0$ 即可；而对于 ${\beta }_1=0,{\beta }_2=0$ 这种情况，需要检验 ${\beta }_1$ 和 ${\beta }_2$ 是否都非零，在此种情况下取 ${\left(MaxP\right)}^2$ 。即有p值计算为：

$p=\{\begin{array}{l}{p}_{1j}={\left(Max{P}_j\right)}^2;{\beta }_1={\beta }_2=0\\ p_{2j}=p_{{\beta }_{1j}};{\beta }_1=0,{\beta }_2\ne 0\\ p_{3j}=p_{{\beta }_{2j}};{\beta }_1\ne 0,{\beta }_2=0\end{array}$ (3.14)

因此原假设的p值的最终计算为：

$p={\omega }_1{p}_{1j}+{\omega }_2{p}_{2j}+{\omega }_3{p}_{3j}$ . (3.15)

4. 模拟检验

4.1. 模拟设置

为了评估这些多效性检验的方法的一类错误以及功效的性能，本文进行广泛的模拟研究。由于所有这些方法都是使用汇总统计量类型的数据集进行的，因此我们直接从给定的二元正态分布在各种子原假设情景下生成Z统计量。具体地说，在 $H_{00}$ 情形下以 ${\pi }_{00}$ 的概率生成Z统计量为 $BVN\left(\left(0,0\right),\Sigma \right)$ ，在 $H_{01}$ 情形下以 ${\pi }_{01}$ 的概率生成Z统计量为 $BVN\left(\left(0,{\mu }_{01}\right),\Sigma \right)$ ，在 $H_{10}$ 情形下以 ${\pi }_{10}$ 的概率生成Z统计量为 $BVN\left(\left({\mu }_{10},0\right),\Sigma \right)$ ，或者在 $H_1$ 情形下以 ${\pi }_{11}$ 的概率生成Z统计量为 $BVN\left(\left({\mu }_1,{\mu }_2\right),\Sigma \right)$ 。同时设置了样本量的数量S为100、1000或2000， ${\mu }_{10}={\mu }_{01}={\mu }_1={\mu }_2$ 且取值为4、10或20，并考虑了多种概率设置，主要以设置 ${\pi }_{11}=0$ 用于检验控制一类错误的情况，以及设置 ${\pi }_{11}\ne 0$ 来检验功效水平。与此同时 ${\mu }_{10}$ (或 ${\mu }_{01}$ 、 ${\mu }_1$ 和 ${\mu }_2$ )的大小量化了基因与性状之间关联的强度大小，值越大表示关联的信号越强。同时基于正态分布的知识近似的将Z统计量转换为对应的p值。在主要的模拟设置中，将 $\Sigma$ 设置为二维的单位矩阵；同时进行了另一项模拟研究，以评估在分析的性状中存在相关关系的情形，研究去相关性后和不去相关性时的这些方法的一类错误和功效的表现情况，来检验相关性将如何影响这些基于原假设的复合多效性测试方法的性能。此时将 $\Sigma$ 的非对角线元素指定为非零值，以生成相关检验统计的数据。

对于最终一类错误的控制和功效的评估，我们在每个设置中重复模拟1000次，最终以这些重复结果得到的平均值作为最终结果。

4.2. 比例估计比较

我们首先通过设置在不同基因个数、不同均值大小等情况下来比较DACT和MAIUP两种方法对于子原假设比例的估计。具体结果见表1~10，表头下标用于区分不同方法。

从表格可以发现，对于DACT方法比例估计的精确度要逊于MAIUP方法获得的比例估计，即估计比例与真实比例之间偏差更大，同时可以看出MAIUP方法在均值大于等于6时估计能够趋于稳定，同时随着基因个数的增加，估计的精确度会更加精准。

Table 1. Estimation of ratio of sub-null hypotheses (100 genes; π 00 = 0.8   ,   π 01 = 0.1   ,   π 10 = 0.1   ,   π 1 = 0 )

表1. 子原假设的比例估计(基因个数为100； ${\pi }_{00}=0.8,{\pi }_{01}=0.1,{\pi }_{10}=0.1,{\pi }_1=0$ )

Table 2. Estimation of ratio of sub-null hypotheses (100 genes; π 00 = 0.04   ,   π 01 = 0.48   ,   π 10 = 0.48   ,   π 1 = 0 )

表2. 子原假设的比例估计(基因个数为100； ${\pi }_{00}=0.04,{\pi }_{01}=0.48,{\pi }_{10}=0.48,{\pi }_1=0$ )

Table 3. Estimation of ratio of sub-null hypotheses (100 genes; π 00 = 1   ,   π 01 = 0   ,   π 10 = 0   ,   π 1 = 0 )

表3. 子原假设的比例估计(基因个数为100； ${\pi }_{00}=1,{\pi }_{01}=0,{\pi }_{10}=0,{\pi }_1=0$ )

Table 4. Estimation of ratio of sub-null hypotheses (100 genes; π 00 = 0   ,   π 01 = 0   ,   π 10 = 0   ,   π 1 = 1 )

表4. 子原假设的比例估计(基因个数为100； ${\pi }_{00}=0,{\pi }_{01}=0,{\pi }_{10}=0,{\pi }_1=1$ )

Table 5. Estimation of ratio of sub-null hypotheses (100 genes; π 00 = 0.1   ,   π 01 = 0.05   ,   π 10 = 0.05   ,   π 1 = 0.8 )

表5. 子原假设的比例估计(基因个数为100； ${\pi }_{00}=0.1,{\pi }_{01}=0.05,{\pi }_{10}=0.05,{\pi }_1=0.8$ )

Table 6. Estimation of ratio of sub-null hypotheses (1000 genes; π 00 = 0.8   ,   π 01 = 0.1   ,   π 10 = 0.1   ,   π 1 = 0 )

表6. 子原假设的比例估计(基因个数为1000； ${\pi }_{00}=0.8,{\pi }_{01}=0.1,{\pi }_{10}=0.1,{\pi }_1=0$ )

Table 7. Estimation of ratio of sub-null hypotheses (1000 genes; π 00 = 0.04   ,   π 01 = 0.48   ,   π 10 = 0.48   ,   π 1 = 0 )

表7. 子原假设的比例估计(基因个数为1000； ${\pi }_{00}=0.04,{\pi }_{01}=0.48,{\pi }_{10}=0.48,{\pi }_1=0$ )

Table 8. Estimation of ratio of sub-null hypotheses (1000 genes; π 00 = 1   ,   π 01 = 0   ,   π 10 = 0   ,   π 1 = 0 )

表8. 子原假设的比例估计(基因个数为1000； ${\pi }_{00}=1,{\pi }_{01}=0,{\pi }_{10}=0,{\pi }_1=0$ )

Table 9. Estimation of ratio of sub-null hypotheses (1000 genes; π 00 = 0   ,   π 01 = 0   ,   π 10 = 0   ,   π 1 = 1 )

表9. 子原假设的比例估计(基因个数为1000； ${\pi }_{00}=0,{\pi }_{01}=0,{\pi }_{10}=0,{\pi }_1=1$ )

Table 10. Estimation of ratio of sub-null hypotheses (1000 genes; π 00 = 0.1   ,   π 01 = 0.05   ,   π 10 = 0.05   ,   π 1 = 0.8 )

表10. 子原假设的比例估计(基因个数为1000； ${\pi }_{00}=0.1,{\pi }_{01}=0.05,{\pi }_{10}=0.05,{\pi }_1=0.8$ )

4.3. 一类错误

在假设检验中一类错误是一个重要的衡量指标，在检验水平为0.01、0.05下，我们设置在不同的均值、基因个数以及各个子原假设在不同比例下进行了相应一类错误的模拟得到如表11、表12的结果。表中New_method表示新构造的p值方法。

Table 11. Comparison of type I error of different methods (significance level is 0.01)

表11. 不同方法下一类错误的比较(显著性水平为0.01)

Table 12. Comparison of type I error of different methods (significance level is 0.05)

表12. 不同方法下一类错误的比较(显著性水平为0.05)

我们发现在两个性状都与同一基因都无关这一情况成立比例较高时，不管均值多大IUT方法都十分的保守，同时DACT方法相对来说也比较保守；但当两个性状都与同一基因都无关这一情况成立比例较低时，IUT方法能够控制一类错误，而DACT方法出现膨胀现象。而PLACO方法只有在三个子假设比例分别为1、0、0时能够控制住一类错误，其余情况出现膨胀状况。新构建的p值、MAIUP方法在三个子假设比例分别为1、0、0时虽然能够控制一类错误，但MAIUP方法相对来说有点膨胀，但新构建的p值的方法在此种情况下表现更好

为了更加直观地进行对比，绘制了在均值为10、20，样本量为1000、2000情况下三个子假设比例分别为1、0、0；0.8、0.1、0.1；0.04、0.48、0.48这三种情形下在显著性水平为0.01的一类错误柱状图。根据图1可以更直观的理解上述的一类错误数据规律。

Figure 1. Comparison of type I error (significance level is 0.01)

图1. 一类错误比较(显著性水平为0.01)

上述结果是在两个性状独立时所得到，现在比较分析两性状存在相关关系时一类错误在去相关性 (图2中C、D)和不去相关性(图2中A、B)大小比较如图2所见，未去相关性时的一类错误随着相关性的加强而变大，此时PLACO、新构建的p值、MAIUP方法尤为明显，在去完相关性后，PLACO方法在检验公式(2.3)中三个子假设比例分别为0.8、0.1、0.1情形下去完相关性后的一类错误有些膨胀(如图2中C所见)。而新构建的p值、MAIUP方法能够较好地维持一类错误，但新构建的p值的方法在D情形下更为稳定。

4.4. 功效

功效是假设检验中另一重要衡量目标，在了解了相关一类错误后，下面我们模拟了不同功效水平不同均值以及基因数下的各个方法的功效大小，如表13、表14。

Figure 2. Comparison of type I error after decorrelation (significance level is 0.01)

图2. 一类错误去相关性对比(显著性水平为0.01)

Table 13. Comparison of power of different methods (significance level is 0.01)

表13. 不同方法下一类功效的比较(显著性水平为0.01)

Table 14. Comparison of power of different methods (significance level is 0.05)

表14. 不同方法下一类功效的比较(显著性水平为0.05)

通过表13、表14发现在备择假设成立的概率越高时，各个方法的功效水平会有所提升，这符合常理。同时发现MAIUP以及New_method所得到的功效水平更高。

同时通过设置备择假设成立的概率在0.7时，三个子假设成立的比例分别为0.1、0.1、0.1的情形，比较在水平为0.01下的不同均值以及样本量间的功效如图3，我们发现在均值较小的时候IUT以及DACT方法功效较差，同时可见New_method以及MAIUP方法能够维持一个稳定功效水平。

Figure 3. Comparison of power (significance level is 0.01)

图3. 功效对比(显著性水平为0.01)

同时将上述方法应用于去相关性的功效中，在备择假设成立的比例为0.7，三个子假设成立的比例分别为0.1、0.1、0.1，均值为4的情形下进行对比(如图4，A表示未去相关性、B表示去相关性)。两性状在有相关性时去相关性前后功效进行了对比，IUT、DACT方法会受相关性的影响，去相关性前后的差距较大。而PLACO、MAIUP方法和新构建的p值的方法则受影响较小，但后两者的功效较高，效果较好。

Figure 4. Comparison of power after decorrelation (significance level is 0.01)

图4. 去相关性功效对比(显著性水平为0.01)

4.5. 实际应用

多效性的显著例子在许多疾病和性状中广泛存在，如自身免疫性疾病 [20] 和精神疾病 [21] 。除此之外多效性也存在于看似不相关的性状中。例如哮喘和精神疾病 [22] 、前列腺癌和血脂 [23] ，以及冠状动脉疾病和扁桃体切除术 [24] 。本文将上述五种方法应用于冠状动脉疾病以及扁桃体切除的实际数据中，检验多

Figure 5. Venn diagram of detected genes of different methods

图5. 不同方法检测出基因的韦恩图

效性基因的个数。从网站(https://grasp.nhlbi.nih.gov/FullResults.aspx)上查找相关的数据并进行下载并随机抽取两千个相关的基因进行多效性检验。由于DACT方法的不及新构建的p值的方法，而MAIUP方法的结果和新构建的p值的方法结果相差无几，因此使用IUT、PLACO、新构建的p值的三种方法来检验多效性基因的个数，分别检测到158、162、420个。并绘制韦恩图如图5。

其中A代表IUT方法、B代表PLACO方法、C代表新构建的p值方法。可以看到新构建的p值的方法能够额外检测到更多多效性基因，也许可以为关研究提供一定参考价值。

5. 结论

本文在此次检验多效性的研究中，以两性状为例，从互斥检验角度出发，介绍了IUT、PLACO、DACT、MAIUP方法，同时在DACT方法的基础上借助MAIUP方法来修正DACT方法即新构造的p值。

通过模拟研究发现，MAIUP方法无论是在维持一类错误还是保持功效水平等方面，都表现更好。IUT方法在检验两性状与一个基因无关时一类错误出现保守现象；PLACO方法在基因与性状之间关联的强度越大时，即关联的信号越强时会不能很好的控制一类错误，一类错误出现膨胀现象，但对于功效能够较好的维持功效水平。因此可以在数据关联的信号强度较小时采用此方法进行检验。DACT方法在检验两性状与一个基因无关这一情形成立概率较大时一类错误会出现保守的现象，而提出的新方法能够较好地控制一类错误，同时较未修正的DACT方法能够提高检验的功效水平，所以表现良好。值得注意，上述情形是基于两性状之间是独立的情形，当两性状之间存在关联时，首先采用去相关性，然后再试用上述方法，此时新方法、MAIUP方法无论是在控制一类错误还是检验功效水平方面都有较好的效果，在实际应用中若出现两性状关联的情形时，检验多效性可以优先考虑新构造的p值的方法和MAIUP方法。