王湘玉，張寶學(xué)，齊春香

（1.河北科技師范學(xué)院 工商管理學(xué)院，河北 秦皇島 066004；2.首都經(jīng)濟(jì)貿(mào)易大學(xué) 統(tǒng)計(jì)學(xué)院，北京 100070；3.東北師范大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，長春 130024）

0 引言

多重假設(shè)檢驗(yàn)是高維統(tǒng)計(jì)推斷的基本問題，也是目前統(tǒng)計(jì)研究的熱點(diǎn)問題之一，應(yīng)用的領(lǐng)域非常廣泛。比如，在金融學(xué)中，需要對數(shù)以萬計(jì)的原假設(shè)同時(shí)進(jìn)行檢驗(yàn)，用以確定哪個(gè)客戶經(jīng)理的能力更強(qiáng)。再如，在全基因組相關(guān)性研究中，希望在大量的基因數(shù)據(jù)中找到與某性狀或疾病相關(guān)的SNP位點(diǎn)，同樣需要同時(shí)對數(shù)以萬計(jì)的原假設(shè)進(jìn)行檢驗(yàn)。而在進(jìn)行多重假設(shè)檢驗(yàn)的過程中，一個(gè)非常重要的問題，就是如何更好地控制總體的錯(cuò)誤率。

Benjamini和Hochberg(1995)[1]提出了多重假設(shè)檢驗(yàn)的FDR標(biāo)準(zhǔn)(E[V/R],即在m個(gè)原假設(shè)中，錯(cuò)誤拒絕的原假設(shè)個(gè)數(shù)V占拒絕的原假設(shè)個(gè)數(shù)R的比例的期望)。之后，Storey(2002)[2]、Efron(2010)[3]等眾多學(xué)者對FDR進(jìn)行了更為深入的研究。然而大部分的研究往往都是基于原假設(shè)相互獨(dú)立為前提。但是這種假設(shè)往往較為嚴(yán)苛，很難達(dá)到。因?yàn)閷?shí)際問題中的原假設(shè)通常具有某種相關(guān)性。舉例來說，如果某個(gè)原假設(shè)是顯著的，那么在它附近的其他原假設(shè)就有更大的可能性也是顯著的。

也有一些學(xué)者針對原假設(shè)相關(guān)的情況做了一些研究，比如Benjamini和Yekutieli（2001）[4]在正回歸相依情況下研究了FDR的估計(jì)問題。再如，Storey等（2004）[5]在弱相關(guān)情況下估計(jì)FDR。但是這些相關(guān)關(guān)系都過于特殊，急需找到一種在原假設(shè)的任意相關(guān)情況下的FDR估計(jì)方法。

Fan等（2012）[6]提出了一種任意協(xié)方差結(jié)構(gòu)下的FDP(False Discovery Proportion，即V/R)估計(jì)方法。然而在估計(jì)隨機(jī)變量W時(shí)，將其當(dāng)作參數(shù)，采用L1估計(jì)方法。實(shí)際上，所得的模型是一個(gè)混合效應(yīng)模型，對此，本文基于最優(yōu)線性無偏預(yù)測方法，估計(jì)隨機(jī)變量W，進(jìn)而提出一種新的FDP估計(jì)方法，最終估計(jì)FDR。

1 最優(yōu)線性無偏預(yù)測

最優(yōu)線性無偏預(yù)測[7]是針對混合效應(yīng)模型，估計(jì)隨機(jī)變量的一種方法。對于混合效應(yīng)模型y=Xβ+Zu+ε,其中y是有N個(gè)觀測的向量；X是N×p的已知矩陣；β是p×1的未知向量（p個(gè)未知常數(shù)），看作固定效應(yīng)部分；Z是N×q的已知矩陣；u是q×1的隨機(jī)向量，看作隨機(jī)效應(yīng)部分；ε是隨機(jī)誤差項(xiàng)，是N×1的隨機(jī)向量。

顯然有 E(u)=0,E(ε)=0 ，定義：

Var(u)=D'Var(ε)=R'Cov(u'ε')=0

因此，V=Var(y)=Var(Zu+ε)=ZDZ'+R

Henderson(1950)[8]得到β和u的最優(yōu)線性無偏預(yù)測：

由公式[9]可得：

2 基于最優(yōu)線性無偏預(yù)測估計(jì)FDP

2.1 理論分析

定義由n個(gè)樣本構(gòu)成的第j個(gè)SNP位點(diǎn)的基因型數(shù)據(jù)為個(gè)樣本的表現(xiàn)型數(shù)據(jù)為Y=(Y1' …'Yn)T?，F(xiàn)考慮和的邊際線性回歸可以得到βj的最小二乘估計(jì)。

這樣可以同時(shí)檢驗(yàn)p個(gè)假設(shè)：

其中，原假設(shè)表示第j個(gè)SNP位點(diǎn)與性狀無關(guān)。

將標(biāo)準(zhǔn)化，記為 Z1'…'Zp，則有其中是 X 相關(guān)系數(shù)矩陣，∑的第(k ' l)個(gè)元素為rkl，對角元素為1。

因此，檢驗(yàn)問題（2）等價(jià)于：

通過對∑進(jìn)行特征分解，得出其中 λ1≥…≥λp是 ∑ 的p個(gè)特征值，γ1'…'γp是對應(yīng)的p個(gè)特征向量，那么這樣 Zi可被分解成：

其中k。W1'…'Wk相互獨(dú)立，且與K1'…'Kp獨(dú)立。W=(W1…,

對于混合效應(yīng)模型（4），當(dāng) y=Z'X=0'Z=L'u=W'D=Ik'R=A時(shí)，帶入式（1）可得：

針對混合效應(yīng)模型（4），取前75%p(記為m)個(gè)最小的對應(yīng)m個(gè)最小的 | μi|，并且將L的分量對應(yīng)排序后，取前m行,從而對應(yīng)有的形式，取A=A11,這樣可以將這些 ||Zi看成是在原假設(shè)式（6）下：

其中Z是有m個(gè)觀測值的m×1向量是m×k矩陣；W 是k×1隨機(jī)向量，是隨機(jī)效應(yīng)部分；K是m×1隨機(jī)向量，是隨機(jī)誤差項(xiàng)。并且有，V=Var(Z)=∑=Var(LW+K)=LL'+A。

Tipping和Bishop(1999)[10]基于高斯?jié)撛谧兞拷嵌忍岢鱿率龈怕誓Ｐ?，此模型體現(xiàn)了主成分分析的思想：

其中

在原假設(shè)下，

因此，將模型（6）替代為模型（7），有：

進(jìn)而由式（5）可得到W 的估計(jì)為：

下面估計(jì)σ2，Z的密度函數(shù)是：

似然函數(shù)為：

對數(shù)似然函數(shù)為：

由得到：

從而得到：

最終可以估計(jì)σ2，帶入式（8）中得出：

這時(shí)，可以得出任意協(xié)方差結(jié)構(gòu)下的FDP估計(jì)：

其中是標(biāo)準(zhǔn)正態(tài)分布的累積分布函數(shù)是標(biāo)準(zhǔn)正態(tài)下的分位數(shù)，ηi的估計(jì)記為。

2.2 模擬及結(jié)果

檢驗(yàn)統(tǒng)計(jì)量現(xiàn)考慮以下六種模型下的協(xié)方差陣結(jié)構(gòu)。其中，Xi產(chǎn)生過程分別如下：

（1）等相關(guān)模型其中 ∑為對角元素是1，其余元素是1/2的矩陣。

（2）Fan和Song[6]的模型令N(0'1)，且：

其中

（3）獨(dú)立柯西模型：對于是獨(dú)立同分布的柯西分布的隨機(jī)變量，其中x0=0'γ=1。

（4）三因子模型：對于 X=(X1'…,Xp)，取 Xj=ρ(1)Wj(1)其 中

（5）雙因子模型：對于 X=(X1'…'Xp)，取其中

（6）非線性因子模型：對于取

其中

令SNP位點(diǎn)個(gè)數(shù) p=1000,n=100'σ=2,錯(cuò)誤的原假設(shè)個(gè)數(shù) p1=50，原假設(shè)下 β0=0，備擇假設(shè)下 β1=1，做1000次模擬。上述六種相關(guān)模型下的FDP估計(jì)值與真值比較結(jié)果如圖1所示。

圖1

圖2為六種模型的FDP估計(jì)值與真值的相對誤差(即

圖2

從模擬結(jié)果可以看出，本文給出的FDP估計(jì)值在真值附近波動(dòng)，理論上合理。

3 結(jié)論

多重假設(shè)檢驗(yàn)是高維統(tǒng)計(jì)推斷的基本問題，也是目前統(tǒng)計(jì)研究的熱點(diǎn)問題之一，應(yīng)用領(lǐng)域十分廣泛。比如，金融領(lǐng)域內(nèi)，要得知哪個(gè)客戶經(jīng)理的能力更強(qiáng)。再如，全基因組相關(guān)性的研究中，要在大量的基因數(shù)據(jù)中尋找與性狀或疾病相關(guān)的SNP位點(diǎn)，往往都需要同時(shí)對數(shù)以萬計(jì)的假設(shè)進(jìn)行檢驗(yàn)。而在進(jìn)行多重假設(shè)檢驗(yàn)問題時(shí)，往往需要控制總體錯(cuò)誤率FDR。棘手的是，實(shí)際問題中的原假設(shè)往往具有一定相關(guān)性。對此，本文提出一種新的估計(jì)FDR的方法，即在檢驗(yàn)統(tǒng)計(jì)量的任意協(xié)方差結(jié)構(gòu)下，對混合效應(yīng)模型使用最優(yōu)線性無偏預(yù)測方法估計(jì)隨機(jī)變量W，從而估計(jì)FDP，進(jìn)而估計(jì)FDR。模擬表明，本文給出的FDP估計(jì)值在真值附近波動(dòng)，理論上比較合理。

但是由于采用最優(yōu)線性無偏預(yù)測的方法，估計(jì)FDP在運(yùn)算速度上不夠快捷，還有待在今后的研究中進(jìn)一步地完善。

[1]Benjamini Y，Hochberg Y.Controlling the False Discovery Rate:A Practical and Powerful Approach to Multiple Testing[J].Journal of the Royal Statistical Society,SeriesB,1995,(57).

[2]Storey J D.,A Direct Approach to False Discovery Rates[J].Journal of the Royal Statistical Society,2002.

[3]Efron B.Correlated Z-Values and the Accuracy of Large-Scale Statistical Estimates[J].Journal of the American Statistical Association,2010,(105).

[4]Benjamini Y，Yekutieli D.The Control of the False Discovery Rate in Multiple Testing Under Dependency[J].The Annals of Statistics,2001,(29).

[5]Storey J D，Taylor J E，Siegmund D.Strong Control,Conservative Point Estimation and Simultaneous Conservative Consistency of False Discovery Rates:A Unified Approach[J].Journal of the Royal Statistical Society,2004,(66).

[6]Fan J,Han X,Gu W.Estimating False Discovery Proportion Under Arbitrary Covariance Dependence[J].Journal of the American Statistical Association,2012.

[7]王松桂，史建紅，尹素菊，吳密霞.線性模型引論[M].北京：科學(xué)出版社,2004.

[8]Henderson C R.Estimation of Genetic Parameters[J].Ann.Math.Statist.,1950,(21).

[9]Benameur S,Mignotte M,Destrempes F,et al.Estimation of Mixtures of Probabilistic PCA with Stochastic EM for the 3D Biplanar Reconstruction of Scoliotic Rib Cage[J].Image Processing,2004,(5).

[10]Tipping M E,Bishop C M.Probabilistic Principal Component Analysis[J].Journal of the Royal Statistical Society,1999.