山西醫(yī)科大學(xué)衛(wèi)生統(tǒng)計(jì)教研室（030001） 許樹紅 王 慧 孫紅衛(wèi) 王 彤

基于LASSO類方法的Ⅰ類錯(cuò)誤的控制*

山西醫(yī)科大學(xué)衛(wèi)生統(tǒng)計(jì)教研室（030001） 許樹紅 王 慧 孫紅衛(wèi) 王 彤△

全基因組關(guān)聯(lián)研究（genome-w ide association studies，GWAS）是在全基因組范圍內(nèi)同時(shí)研究上百萬個(gè)單核苷酸多態(tài)性（single nucleotide polymorphism，SNP）位點(diǎn)與疾病或某些性狀之間的關(guān)聯(lián)，從而篩選出可能的致病SNP位點(diǎn)，進(jìn)而對(duì)這些位點(diǎn)進(jìn)行人群驗(yàn)證和實(shí)驗(yàn)分析。在GWAS研究中比較傳統(tǒng)的分析方法是針對(duì)每個(gè)SNP和結(jié)局變量間關(guān)聯(lián)進(jìn)行單因素分析的假設(shè)檢驗(yàn)，而待分析的SNP數(shù)量有幾十萬甚至上百萬個(gè)，使得檢驗(yàn)次數(shù)十分巨大，如果不采用合適的方法進(jìn)行多重性校正來妥善控制Ⅰ類錯(cuò)誤，會(huì)產(chǎn)生許多假陽性結(jié)果，對(duì)這些結(jié)果進(jìn)行驗(yàn)證將耗費(fèi)很多時(shí)間和財(cái)力，造成不必要的損失。針對(duì)全基因組測(cè)序數(shù)據(jù)進(jìn)行多因素建模常采用的分析策略為降維和變量選擇。變量選擇的方法包括經(jīng)典方法和懲罰類方法，懲罰類方法在GWAS研究、測(cè)序數(shù)據(jù)分析中應(yīng)用廣泛且發(fā)展迅速，它假定模型具有稀疏性即真實(shí)情況下許多未知候選變量的回歸系數(shù)為0或者接近0，因此自變量的選擇問題轉(zhuǎn)化為確定正確的子模型問題，且在選擇變量的同時(shí)給出模型參數(shù)的估計(jì)。LASSO（least absolute shrinkage and selection operator）作為一種常用的懲罰類方法，部分研究者指出其選出的變量存在較高的假陽性問題［1-3］。針對(duì)該類問題，在多重校正以及變量選擇方法的基礎(chǔ)上，發(fā)展了一些控制Ⅰ類錯(cuò)誤的同時(shí)篩選出正確變量的方法。本文主要對(duì)多重檢驗(yàn)中控制Ⅰ類錯(cuò)誤的方法進(jìn)行總結(jié)并且對(duì)基于LASSO的Ⅰ類錯(cuò)誤的控制方法進(jìn)行綜述。

多重檢驗(yàn)中控制Ⅰ類錯(cuò)誤的策略和方法

經(jīng)典的多重檢驗(yàn)校正是建立在多個(gè)檢驗(yàn)之間或P值之間相互獨(dú)立的假設(shè)上［4］，然而在針對(duì)基因數(shù)據(jù)的多個(gè)假設(shè)檢驗(yàn)中得到的P值之間往往是不獨(dú)立的，目前處理這種不獨(dú)立情況的方法大致分為三類：第一種是不做任何假設(shè)；第二種是假設(shè)P值之間存在PDS假設(shè)（positive dependence through stochastic ordering，PDS）［5-7］，PDS假設(shè)有弱PDS假設(shè)和強(qiáng)PDS假設(shè)兩種，前者為：對(duì)每一個(gè)成立的H0對(duì)應(yīng)的P值qi，及所有qi的不減函數(shù)f，在qi＝u的條件下，E [f（q1，…，qm0）｜qi＝u]是不減的，后者為：對(duì)每一個(gè)H0對(duì)應(yīng)的P值pi及其不減函數(shù)f，在qi＝u的條件下，E [f（p1，…，pm）｜qi＝u]是不減的；第三種是基于排列抽樣的方法。目前常用的控制Ⅰ類錯(cuò)誤的方法有兩大類，分別是控制FWER（family-w ise error rate，F(xiàn)WER）和控制FDR（false discovery rate，F(xiàn)DR）。

設(shè)進(jìn)行m次假設(shè)檢驗(yàn)，原假設(shè)為H0，i，i＝1，…，m，對(duì)應(yīng)的P值表示為pi，i＝1，…，m。將成立的H0對(duì)應(yīng)的P值定義為qi，i＝1，…，m0。表1是各種可能的檢驗(yàn)結(jié)果，m0和m1分別是檢驗(yàn)中成立的H0及不成立的H0的個(gè)數(shù)。V是檢驗(yàn)中成立的H0被錯(cuò)誤拒絕的個(gè)數(shù)。π0＝m0／m，是m次假設(shè)檢驗(yàn)中成立的H0所占比例。

表1 多重檢驗(yàn)的結(jié)果

1．控制FWER

FWER定義為至少犯一次Ⅰ類錯(cuò)誤的概率，F(xiàn)WER＝P（V≥1），是傳統(tǒng)的多重檢驗(yàn)中廣泛應(yīng)用的控制Ⅰ類錯(cuò)誤的指標(biāo)?？刂艶WER的方法有經(jīng)典的Bonferroni方法，當(dāng)pi＜α／m時(shí)，拒絕相應(yīng)的假設(shè)，公式1證明了Bonferroni方法可將FWER控制在π0α水平。

由此，當(dāng)真實(shí)的H0全部成立時(shí)控制FWER是弱控制，此時(shí)FWER被控制在α水平；當(dāng)真實(shí)的H0部分成立時(shí)控制FWER是強(qiáng)控制，且若不成立的H0的個(gè)數(shù)很多，F(xiàn)WER將被控制在遠(yuǎn)低于α的水平，此時(shí)Bonferroni方法便會(huì)很保守［7］。因此Holm在1979年提出Holm方法［8］，在Bonferroni方法基礎(chǔ)上，迭代進(jìn)行檢驗(yàn)，第一步拒絕pi＜α／h0的假設(shè)，h0＝m；第二步排除掉第一步已經(jīng)拒絕的假設(shè)，剩余沒有被拒絕的假設(shè)數(shù)為h1，拒絕pi＜α／h1的假設(shè)；以此類推，直到?jīng)]有假設(shè)被拒絕為止。Holm證明該方法能將FWER控制在α水平。Holm方法的另一種解法是首先將P值從小到大進(jìn)行排序即p（1）≤p（2）≤…≤p（m），將p（i）與α／(m-i+1)進(jìn)行比較，得到最小的j滿足p（j）＞α／(m-j+1)，然后拒絕所有j-1假設(shè)。Bonferroni和Holm方法都沒有對(duì)P值之間的關(guān)系做假設(shè)，即對(duì)P值的聯(lián)合分布沒有假設(shè)。Hochberg［9］利用成立的H0滿足弱PDS假設(shè)，在Holm方法基礎(chǔ)上進(jìn)行了修改，Hochberg方法是首先將P值從小到大進(jìn)行排序即p（1）≤p（2）≤…≤p（m），將p（i）與α／(m-i+1)進(jìn)行比較，得到最大的j滿足p（j）＜α／(m-j+1)，然后依次拒絕所有j假設(shè)。Westfall＆Young［10］利用排列抽樣方法來控制FWER，其原理是在H0假設(shè)成立條件下，樣本來自相同的總體，可以對(duì)原始樣本進(jìn)行重新抽樣，得到統(tǒng)計(jì)量的分布并進(jìn)行檢驗(yàn)。排列抽樣方法應(yīng)用很廣，不僅可以用來控制FWER，也可以控制FDR等。Westfall＆Young的方法雖然在P值之間不獨(dú)立時(shí)效能比前面提到的方法高，但是排列抽樣對(duì)原假設(shè)間的一致性有要求。

隨著基因工程和下一代測(cè)序技術(shù)的迅猛發(fā)展，產(chǎn)生龐大的生物信息數(shù)據(jù)，通過在多重檢驗(yàn)中控制FWER來控制Ⅰ類錯(cuò)誤，對(duì)于探索性研究而言結(jié)果太過保守，F(xiàn)WER將每個(gè)假設(shè)檢驗(yàn)的可信度控制在1-α水平，這只能識(shí)別出少數(shù)有意義的基因或者變量，且在嚴(yán)格控制Ⅰ類錯(cuò)誤的同時(shí)其Ⅱ類錯(cuò)誤偏大。它的優(yōu)勢(shì)在于對(duì)每一個(gè)假設(shè)都嚴(yán)格控制Ⅰ類錯(cuò)誤，在其選出的基因集中任意子集都對(duì)Ⅰ類錯(cuò)誤有很好的控制，適用于對(duì)每個(gè)假設(shè)均有要求的驗(yàn)證性研究。

2．控制FDR

FDR與FWER之間的關(guān)系為：當(dāng)m0＝m時(shí)，即全部原假設(shè)均成立，此時(shí)FDR的控制與FWER的控制等價(jià)；當(dāng)m0＜m時(shí)，F(xiàn)DR小于或者等于FWER，因此，F(xiàn)DR相當(dāng)于弱控制FWER［11-12］。

B-H方法可簡(jiǎn)單表述為：先將P值從小到大進(jìn)行排序即p（1）≤p（2）≤…≤p（m），然后將p（i）與iα／m進(jìn)行比較，將FDR控制在α水平，定義j＝max {i：p（i）≤iα／m}，然后拒絕所有j假設(shè)，若沒有滿足條件的i存在，則沒有假設(shè)被拒絕。Benjamini＆Hochberg已證明當(dāng)P值間相互獨(dú)立的時(shí)候，B-H方法能將FDR控制在π0α水平而不是α水平，故其結(jié)果有些保守。大量文獻(xiàn)在B-H的方法基礎(chǔ)上進(jìn)行了改進(jìn)，如將iα／m換為iα／（π＾0m），其中π＾0是π0的估計(jì)值，并針對(duì)π0的估計(jì)提出很多方法，該類方法為“自適應(yīng)”方法［13-14］。原始的B-H方法要求P值間相互獨(dú)立，而當(dāng)P值間不獨(dú)立時(shí)，許多研究者提出了相應(yīng)的改進(jìn)方法以及一些基于重抽樣的方法，如Yekutieli提出了ssBH方法［15］，Benjam ini＆Yekutieli提出基于排列抽樣的方法但是沒有完全證明該方法能控制FDR［16］，Joseph［17］等以及Troendle［18］提出了基于bootstrap的方法也能漸近控制FDR，但bootstrap的方法要求真實(shí)H0對(duì)應(yīng)的P值的聯(lián)合分布的極限滿足可交換性［19］。除了控制FDR，還有一類方法是對(duì)選出的結(jié)果集的FDP（false discovery proportion）進(jìn)行估計(jì)。FDR控制與FDP的估計(jì)是有區(qū)別的。例如：某個(gè)研究者重復(fù)一個(gè)實(shí)驗(yàn)z次，采用控制FDR的方法進(jìn)行分析，那么可能得到z個(gè)不同的結(jié)果集Ri，i＝1，…，z，每一個(gè)結(jié)果集都有其未知的FDP，記為Qi。當(dāng)Z很大時(shí)，Qi的期望將接近FDR的控制水平。Q＾i為第i個(gè)結(jié)果集的FDP的估計(jì)值［7］。FDP的估計(jì)方法大多是從貝葉斯角度解釋，幾個(gè)代表性的方法是Storey 2003年提出的q-value法、Tusher 2001年提出的SAM（significance analysis ofm icroarray）方法、Efron 2001年提出的基于經(jīng)驗(yàn)貝葉斯的方法等，2012年劉晉等對(duì)這類方法進(jìn)行了比較詳細(xì)的介紹［12］。在P值間不獨(dú)立的情況下，F(xiàn)an等［20-21］在Storey［22］和Friguet［23］基礎(chǔ)上，提出PFA（principal factor approximation）法，該方法在進(jìn)行多重檢驗(yàn)校正時(shí)，通過因子分析將變量間相關(guān)關(guān)系考慮在內(nèi)，提出適用于任意不獨(dú)立情況的FDP的近似公式和一致估計(jì)。除此之外，一些文獻(xiàn)對(duì)FDP的置信區(qū)間進(jìn)行了估計(jì)［24-25］，如Meinshausen基于排列抽樣的方法。

目前在基因數(shù)據(jù)分析領(lǐng)域的多重比較問題中，F(xiàn)DR已經(jīng)成為一個(gè)標(biāo)準(zhǔn)有效的反映Ⅰ類錯(cuò)誤的指標(biāo)，適用于探索性研究，在整體上控制Ⅰ類錯(cuò)誤的同時(shí)能篩選出更多有意義的基因。然而FDR也存在一定的局限性，其一是FDR是隨機(jī)變量Q的期望值，當(dāng)P值間不獨(dú)立時(shí)或者一些FDR估計(jì)方法的假設(shè)不滿足時(shí)，分析結(jié)果中Q＾的變異可能會(huì)增大［7，12］；其二是控制FDR的條件下得到的分析結(jié)果是從整體上控制了Ⅰ類錯(cuò)誤，無法保證控制了每個(gè)假設(shè)或者每個(gè)結(jié)果的子集的Ⅰ類錯(cuò)誤。因此研究者在解釋其結(jié)果時(shí)需要注意保持結(jié)果集的完整性［7］。

基于LASSO類方法的控制Ⅰ類錯(cuò)誤的策略和技術(shù)

隨著生物信息技術(shù)的快速發(fā)展、大數(shù)據(jù)時(shí)代的到來，高維數(shù)據(jù)分析問題不斷涌現(xiàn)（本文中高維數(shù)據(jù)是指p＞n的情況，p和n分別表示自變量個(gè)數(shù)和樣本量），面對(duì)復(fù)雜數(shù)據(jù)建模時(shí)如何更加準(zhǔn)確的篩選變量成為很重要的問題。隨著計(jì)算能力的提高和算法的不斷改進(jìn)，懲罰類方法因其通過好的算法能執(zhí)行較大維度的變量選擇過程、得到稀疏解而逐漸受到重視，是目前比較流行的變量選擇方法。該類方法假定模型具有稀疏性，即真實(shí)模型中很多自變量的回歸系數(shù)為0，其思想是在最小化損失函數(shù)或者最大化似然函數(shù)的同時(shí)增加一個(gè)懲罰項(xiàng)來篩選變量并估計(jì)其回歸系數(shù)，將回歸系數(shù)不為零的變量直接選入結(jié)果集中。常見的懲罰類方法有：LASSO、SCAD（smoothly clipped absolute deviation）、彈性網(wǎng)（elastic net）、MCP（minimax concave penalty）等。

1．LASSO中關(guān)于調(diào)整參數(shù)選擇的傳統(tǒng)方法

基于LASSO的變量選擇及回歸系數(shù)的估計(jì)，在很大程度上受調(diào)整參數(shù)λ的影響，因?yàn)檎{(diào)整參數(shù)值的大小影響著模型的復(fù)雜程度及模型的收斂速度，因此選擇合適的調(diào)整參數(shù)十分重要［26-27］。傳統(tǒng)的選擇調(diào)整參數(shù)的方法大致分為兩類：交叉驗(yàn)證（cross validation，CV）和信息準(zhǔn)則（information criteria）。

（1）交叉驗(yàn)證

最常用的交叉驗(yàn)證是K折交叉驗(yàn)證（K-fold cross-validation），K為整數(shù)，1＜K≤n。該方法的過程是：首先將原始樣本Γ拆分成K個(gè)相同大小的子樣本，將其中一個(gè)子樣本選為驗(yàn)證數(shù)據(jù)集Γv，其余K-1個(gè)子樣本為訓(xùn)練數(shù)據(jù)集Γ-Γv，利用訓(xùn)練數(shù)據(jù)集來建立模型，然后用擬合模型來預(yù)測(cè)驗(yàn)證數(shù)據(jù)集，計(jì)算驗(yàn)證集的預(yù)測(cè)值與真實(shí)值間的偏差，重復(fù)以上過程K次，計(jì)算平均偏差，選擇平均偏差最小時(shí)對(duì)應(yīng)的調(diào)整參數(shù)值。

基于最小二乘的K折交叉驗(yàn)證的目標(biāo)函數(shù)為［28］：

其中，β＾（-v）（λ）表示去掉（yk，xk）∈Γv后的回歸系數(shù)估計(jì)值，使CV（λ）取得最小值時(shí)對(duì)應(yīng)的λ為最終模型的調(diào)整參數(shù)。

基于對(duì)數(shù)似然函數(shù)的K折交叉驗(yàn)證的目標(biāo)函數(shù)為［27］：

其中，β＾（-v）（λ）表示去掉第v份數(shù)據(jù)后的回歸系數(shù)估計(jì)值，L（v）（β＾（-v）（λ））表示第v份數(shù)據(jù)為驗(yàn)證集時(shí)求得的對(duì)數(shù)似然函數(shù)值，使CVl（λ）取得最大值時(shí)對(duì)應(yīng)的λ為最終模型的調(diào)整參數(shù)。

當(dāng)K＝n時(shí)的交叉驗(yàn)證稱為留一法（leave-oneout），該方法是對(duì)真實(shí)的預(yù)測(cè)誤差的漸近無偏估計(jì)，但是該法方差很大且計(jì)算量龐大［26］。Tibshirani1996年使用廣義交叉驗(yàn)證（generalized cross validation，GCV）選擇調(diào)整參數(shù)，GCV在CV的基礎(chǔ)上將有效參數(shù)的個(gè)數(shù)考慮在內(nèi)估計(jì)預(yù)測(cè)誤差。

（2）信息準(zhǔn)則

信息準(zhǔn)則一般包括模型擬合程度的測(cè)量（如似然函數(shù)或損失函數(shù)）以及對(duì)模型復(fù)雜程度的懲罰（懲罰乘子an）兩部分［29］。對(duì)于一個(gè)給定的統(tǒng)計(jì)模型如線性模型、生存模型，不同的信息準(zhǔn)則方法的主要區(qū)別在于如何對(duì)模型的復(fù)雜程度進(jìn)行懲罰，然而，對(duì)于選擇真實(shí)的模型而言，如何通過對(duì)模型復(fù)雜程度進(jìn)行合理的懲罰來優(yōu)化信息準(zhǔn)則十分關(guān)鍵［30］。常用的信息準(zhǔn)則有：AIC（Akaike information criterion）、BIC（Bayesian information criterion）、RIC（Risk inflation criterion）［31-33］等。an＝2和an＝log（n）分別對(duì)應(yīng)的是AIC和BIC，an＝2log（p）為RIC，其中n為樣本量，p為非零回歸系數(shù)的個(gè)數(shù)。Yang和Wang指出CV以及GCV類似于傳統(tǒng)的AIC，三種方法在高維數(shù)據(jù)情形下均易出現(xiàn)過擬合，且AIC不是一致的信息準(zhǔn)則［34-36］。Wang［36］的研究表明針對(duì)SCAD，BIC能一致的選擇真實(shí)的模型，但是隨著樣本量增大，BIC得到的結(jié)果可能會(huì)太保守，易遺漏一些重要的變量［37］。

2．線性模型中基于LASSO類方法的控制Ⅰ類錯(cuò)誤的方法

假設(shè)有n個(gè)樣本，y表示研究的反應(yīng)變量，X為n ×p維的自變量矩陣，β是回歸系數(shù)，z是隨機(jī)誤差，z～N 0，σ2

()I。線性模型為：

基于LASSO的目標(biāo)函數(shù)為：

其中b是β的估計(jì)值，λ是調(diào)整參數(shù)。

（1）平穩(wěn)選擇（stability selection）

平穩(wěn)選擇法是將變量選擇算法和再抽樣技術(shù)相結(jié)合的方法，由Meinshausen和Bühlmann［42］提出，通過計(jì)算機(jī)對(duì)樣本進(jìn)行隨機(jī)無放回抽樣得到樣本量為n／2的子樣本I，然后利用LASSO等懲罰回歸方法進(jìn)行變量選擇，選入模型的變量估計(jì)為β＾λj()I，選入模型的變量集為：

對(duì)不同的λ，每個(gè)變量j被選入模型的概率為：

平穩(wěn)選擇法將其中被選取概率大于閾值的變量選入最終模型，最終變量集為：

平穩(wěn)選擇法在R軟件中可以實(shí)現(xiàn)。一般隨機(jī)抽取子樣本數(shù)達(dá)到500到1000即可，其計(jì)算成本較交叉驗(yàn)證相差不多甚至更低［43-45］。該方法不單獨(dú)使用，而是與其他變量選擇方法結(jié)合使用，應(yīng)用于不同模型中。該法通過再抽樣技術(shù)降低了模型結(jié)果對(duì)調(diào)整參數(shù)選擇的敏感性，并達(dá)到了控制Ⅰ類錯(cuò)誤的目的。但是由于基因之間往往存在一定的關(guān)聯(lián)，其選擇變量集可交換的假設(shè)不一定成立，且有文獻(xiàn)報(bào)道在真實(shí)有意義的變量數(shù)極少的情況下，該方法的FDR值有時(shí)會(huì)很高，失去其控制Ⅰ類錯(cuò)誤的作用［37］。

（2）SLOPE

SLOPE（sorted L-one penalized estimation）方法［3，46-48］，該法是在LASSO的基礎(chǔ)上結(jié)合控制FDR的B-H方法提出的。其模型是：

其中λi表示第i個(gè)變量的調(diào)整參數(shù)值，b（i）表示第i個(gè)變量的回歸系數(shù)絕對(duì)值，λi滿足λ1≥λ2≥…≥λp≥0降序排列，且λi的順序根據(jù)λBH（i）來決定，相應(yīng)的

SLOPE方法首先在X是正交矩陣的條件下（XTX＝Ip），根據(jù)線性模型公式5構(gòu)造y～，模型為：

其中y～∈N(β，σ2Ip)，當(dāng)時(shí)拒絕原假設(shè)。其次，在處理高維數(shù)據(jù)的LASSO方法基礎(chǔ)上進(jìn)行改進(jìn)，不同效應(yīng)的變量對(duì)應(yīng)的調(diào)整參數(shù)值不同，令λBH（i）＝z（1-i·α／2p），其中z(α)為標(biāo)準(zhǔn)正態(tài)分布的分位數(shù)，λi＝λBH（i）·σ，根據(jù)λi來選擇估計(jì)變量。Bogdan等證明在X滿足正交矩陣的條件下SLOPE方法能將FDR控制在α水平內(nèi)。

當(dāng)X不滿足正交矩陣的條件時(shí)，Bogdan等對(duì)λBH（i）進(jìn)行了調(diào)整，在計(jì)算某一變量的λBH（i）時(shí)將其他變量的影響考慮在內(nèi)，提出λG(i)：

其中k*＝k( n，p，q)為最小的λG(i)的位置，λk*為λG(i)的最小值。

SLOPE借鑒多重檢驗(yàn)校正中B-H方法的思想（B-H方法對(duì)P值從小到大排序，越小的p（i）其對(duì)應(yīng)的檢驗(yàn)水準(zhǔn)iα／m值越小即檢驗(yàn)水準(zhǔn)越嚴(yán)格）改進(jìn)LASSO，對(duì)不同效應(yīng)的變量檢驗(yàn)水準(zhǔn)不同，效應(yīng)越大的變量相應(yīng)的檢驗(yàn)水準(zhǔn)越嚴(yán)格，其目的是為了控制FDR。SLOPE與LASSO的共同點(diǎn)在于兩者對(duì)β＾i的估計(jì)均為有偏估計(jì)，不同點(diǎn)在于LASSO法只選擇一個(gè)合理的λ值，根據(jù)估計(jì)β，Zi不為0，不同變量受到相同的懲罰。SLOPE根據(jù)λi來選擇估計(jì)變量，不同的變量受到的懲罰不同，值越大的變量λi越大。鑒于SLOPE是有偏估計(jì)的方法，Bogdan等在其文章中并不推薦利用SLOPE對(duì)進(jìn)行直接估計(jì)，而是建議使用兩階段法，第一階段采用SLOPE法篩選出有意義的變量，第二階段對(duì)篩選出的變量利用普通最小二乘法估計(jì)其回歸系數(shù)。SLOPE方法在R、Matlab軟件中可以實(shí)現(xiàn)。

（3）Knockoffs filter

Barber和Candès［49］提出一種新的在有限樣本情況下篩選變量的方法，對(duì)任何自變量矩陣X∈Rn×p，n≥p，構(gòu)造相應(yīng)的knockoff變量矩陣X～，并將其與原始變量一同放入線性模型中，通過構(gòu)建關(guān)于原始變量Xj與knockoff變量X～j的統(tǒng)計(jì)量Wj、確定統(tǒng)計(jì)量的界值T來控制線性模型中的FDR。

Knockoffs filter方法的過程分為三部分：

①對(duì)每一個(gè)原始變量Xj，j＝1，…，p構(gòu)造knockoff變量X～j，要求knockoff保持原始變量之間的相關(guān)結(jié)構(gòu)，即滿足XTX＝X～TX～＝∑，XTjX～k＝XTjXk，j≠k，且假設(shè)∑矩陣可逆，同時(shí)要求knockoff變量X～j與原始變量Xj間相似程度盡可能地低，即對(duì)于XTX～＝∑－diag()s盡量選擇大的s值。

②每對(duì)原始變量Xj和knockoff變量X～j構(gòu)建一個(gè)統(tǒng)計(jì)量Wj。Barber和Candès在2014年的文章中給出很多種構(gòu)建統(tǒng)計(jì)量Wj的方式，本文只詳述其中之一說明Wj的意義。在LASSO公式（6）基礎(chǔ)上，將自變量矩陣X換為增廣矩陣X X[]～，定義Zj為原始變量Xj首次進(jìn)入模型時(shí)的λ值，Zj＝supλ：β＾j()λ≠{

}0，相應(yīng)的Z～j為變量X～j首次進(jìn)入模型時(shí)的λ值，對(duì)于每對(duì)原始變量Xj和knockoff變量X～j，構(gòu)建統(tǒng)計(jì)量Wj：

當(dāng)Zj＝Z～j時(shí)，Wj＝0，若Wj為正并且值越大，說明原始變量Xj比knockoff變量X～j越早進(jìn)入模型，則越有理由拒絕原假設(shè)。

③定義統(tǒng)計(jì)量的界值T，將滿足Wj≥T的自變量選入模型。設(shè)FDR控制水平為q，在原假設(shè)βj＝0成立情況下，原始變量Xj與knockoff變量X～j兩者誰先進(jìn)入模型是隨機(jī)性的，即Wj的正負(fù)號(hào)是隨機(jī)的，那么對(duì)于任意的t，＃｛j：βj＝0 and Wj≥t｝＝d＃｛j：βj＝0 and Wj≤－t｝，其中＝d代表同分布。因此在t時(shí)，F(xiàn)DP的knockoff估計(jì)值為：

根據(jù)FDP的knockoff估計(jì)值定義T，T＝m in｛t∈W：FDP＾（t）≤q｝，最終界值T為：

Knockoffs filter方法的關(guān)鍵之一在于knockoff變量的構(gòu)造，當(dāng)樣本量n≥p時(shí)，Barber和Candès描述了如何構(gòu)造X～以及如何選擇合適的s值，但是當(dāng)樣本量n＜p時(shí)，∑矩陣可逆的假設(shè)不再滿足，該方法不能直接使用，Barber和Candès［49］建議可以先使用一些降維的方法，再使用Knockoffs filter方法。Reid和Tibshirani［50］首先使用聚類處理自變量間相關(guān)的問題實(shí)現(xiàn)降維，然后用每個(gè)類別中選出的代表性的變量建模通過Knockoffs filter方法控制FDR。Barber和Candès［51］新提出一種分析策略應(yīng)用于樣本量n＜p的情況，該策略將觀測(cè)分成兩組，第一組觀測(cè)的數(shù)據(jù)用于篩選潛在的有意義的自變量，其次在篩選出的變量基礎(chǔ)上用第二組觀測(cè)的數(shù)據(jù)，通過Knockoffs filter方法進(jìn)行統(tǒng)計(jì)推斷。目前Knockoffs filter方法可以通過R軟件實(shí)現(xiàn)。

3．Cox模型中基于LASSO的控制FDR的方法

基于Cox模型的LASSO，其最大化目標(biāo)函數(shù)為：

其中l(wèi)β，()X是Cox模型的對(duì)數(shù)似然函數(shù)，p為自變量個(gè)數(shù)。LASSO作為生存分析領(lǐng)域中常用的處理高維數(shù)據(jù)的方法之一，篩選的變量結(jié)果中存在FDR過高的問題，而調(diào)整參數(shù)的選擇很大程度上影響了LASSO的變量選擇和回歸系數(shù)的估計(jì)。Ternès，Rotolo和M ichiels［37］針對(duì)傳統(tǒng)的交叉驗(yàn)證法容易出現(xiàn)過擬合的問題，提出了一種改進(jìn)的調(diào)整參數(shù)的選擇方法，即pcvl（penalized cross-validated log-likelihood，pcvl），該方法主要在傳統(tǒng)交叉驗(yàn)證的目標(biāo)函數(shù)中加入懲罰項(xiàng)，以實(shí)現(xiàn)模型在滿足擬合優(yōu)度與模型稀疏程度兩者之間的折中。

Ternès，Rotolo和M ichiels文章中傳統(tǒng)K折交叉驗(yàn)證的目標(biāo)函數(shù)是：

其中β＾（－v）（λ）表示去掉第v份數(shù)據(jù)后的回歸系數(shù)估計(jì)值，L（β＾（－v）（λ），X）表示全部數(shù)據(jù)求得的對(duì)數(shù)似然函數(shù)值，L（－v）（β＾（－v）（λ），X（－v））表示去掉第v份數(shù)據(jù)時(shí)求得的對(duì)數(shù)似然函數(shù)值，使cvl（λ）取得最大值時(shí)對(duì)應(yīng)的λ為模型的調(diào)整參數(shù)，記為λ＾cvl。該函數(shù)參考Verweij和van Houwelingen［52］文章中對(duì)生存分析中交叉驗(yàn)證法的介紹，與公式（4）的差別在于后者是利用全部數(shù)據(jù)的對(duì)數(shù)似然函數(shù)值與去掉第v份數(shù)據(jù)時(shí)求得的對(duì)數(shù)似然函數(shù)值之差來求第v份數(shù)據(jù)的對(duì)數(shù)似然函數(shù)，而且并沒有取平均值。

令λ0表示沒有變量進(jìn)入模型時(shí)λ的最小值，即λ0＝minλ｜p{}＝0，懲罰項(xiàng)為：

Ternès，Rotolo和M ichiels將pcvl方法與傳統(tǒng)CV法、AIC、RIC、自適應(yīng)LASSO和平穩(wěn)選擇法進(jìn)行了比較，模擬數(shù)據(jù)結(jié)果顯示pcvl方法的優(yōu)勢(shì)在于控制FDR的情況下也能將假陰性率控制在較低水平。但是該方法是從適當(dāng)增加λ值以減少LASSO篩選出的變量數(shù)的角度控制FDR，無法在建模時(shí)給定FDR的理想水平實(shí)現(xiàn)精確的FDR的控制。

表2 基于LASSO類方法的控制Ⅰ類錯(cuò)誤方法的比較

展 望

本文首先回顧了多重假設(shè)檢驗(yàn)中Ⅰ類錯(cuò)誤控制方法，根據(jù)Ⅰ類錯(cuò)誤控制指標(biāo)的不同分為兩大類：控制FWER的方法和控制FDR的方法，然而基因間不獨(dú)立的問題對(duì)多重檢驗(yàn)結(jié)果產(chǎn)生了一定的影響，同時(shí)在多變量情況下如何處理高維數(shù)據(jù)的問題也亟待解決，因此一些變量選擇方法如懲罰類方法成為如今研究的熱點(diǎn)。本文針對(duì)LASSO這一變量選擇方法存在的FDR過高的問題，介紹了四種方法，分別為平穩(wěn)選擇法、SLOPE、Knockoffs filter和pcvl方法，已有一些研究分別將平穩(wěn)選擇法、SLOPE和pcvl方法與傳統(tǒng)的調(diào)整參數(shù)選擇方法進(jìn)行了比較且結(jié)果顯示平穩(wěn)選擇法、SLOPE和pcvl方法可以有效控制Ⅰ類錯(cuò)誤，但這些方法各有優(yōu)劣及適用條件，并非在任何數(shù)據(jù)情況下都有效，因此在應(yīng)用這些方法時(shí)需注意適用條件是否滿足。其中一些方法如平穩(wěn)選擇、Knockoffs filter等也可以擴(kuò)展到其他懲罰回歸模型中。LASSO懲罰估計(jì)是有偏的，不滿足oracle的三個(gè)性質(zhì)，未來如何將這些控制Ⅰ類錯(cuò)誤的方法應(yīng)用于其他懲罰類方法如SCAD、MCP等還需要進(jìn)一步研究。

［1］Benner A，Zucknick M，Hielscher T，et al．High-dimensional Cox models：the choice of penalty as part of themodel building process．Biometrical Journal，2010，52（1）：50-69．

［2］Shojaie A，M ichailidis G．Penalized likelihood methods for estimation of sparse high-dimensional directed acyclic graphs．Biometrika，2010，97（3）：519-538．

［3］Bogdan M，Berg EVD，Su W，et al．Statistical estimation and testing via the sorted L1 norm．a(chǎn)rXiv：1310．1969，2013．

［4］王彤，易東．臨床試驗(yàn)中多重性問題的統(tǒng)計(jì)學(xué)考慮．中國衛(wèi)生統(tǒng)計(jì)，2012，29（3）：445-450．

［5］Block HW，Savits TH，Shaked M．A concept of negative dependence using stochastic ordering．Statistics＆Probability Letters，1985，3（2）：81-86．

［6］Sarkar SK，Guo W．On a generalized false discovery rate．Mathematics，2009，37（3）：1545-1565．

［7］Goeman JJ，Solari A．Multiple hypothesis testing in genomics．statistics in medicine，2014，33（11）：1946-1978．

［8］Holm S．A simple sequentially rejectivemultiple test procedure Scandinavian Journal of Statistics，1979，6（2）：65-70．

［9］Hochberg Y．A sharper Bonferroniprocedure formultiple tests of significance．Biometrika，1988，75（4）：800-802．

［10］Westfall PH，Young SS．Resampling-Based Multiple Testing：Examples and Methods for P-Value Adjustment．NewYork：W iley-Interscience，1993．

［11］Benjam ini Y，Hochberg Y．Controlling the False Discovery Rate：APractical and Powerful Approach to Multiple Testing．Journal of the Royal Statistical Society，1995，57（1）：289-300．

［12］劉晉，張濤，李康．多重假設(shè)檢驗(yàn)中FDR的控制與估計(jì)方法．中國衛(wèi)生統(tǒng)計(jì)，2012，29（2）：305-308．

［13］Benjam ini Y，Hochberg Y．On the Adaptive Control of the False Discovery Rate in Multiple Testing w ith Independent Statistics．Journal of Educational And Behavioral Statistics，2000，25（1）：60-83．

［14］Benjam ini Y，Krieger AM，Yekutieli D．Adaptive linear step-up procedures that control the false discovery rate．Biometrika，2006，93（3）：491-507．

［15］Yekutieli D．False discovery rate control for non-positively regression dependent test statistics．Journalof Statistical Planning and Inference，2007，138（2）：405-415．

［16］YekutieliD，Benjam ini Y．Resampling-based false discovery rate controlling multiple test procedures for correlated test statistics．Journal of Statistical Planning and Inference，1999，82（1-2）：171-196．

［17］Romano JP，Shaikh AM，Wolf M．Control of the false discovery rate under dependence using the bootstrap and subsampling．Test，2008，17（3）：417-442．

［18］Troendle JF．Stepw ise normal theory multiple test procedures controlling the false discovery rate．Journal of Statistical Planning and Inference，2000，84（1-2）：139-158．

［19］劉莉．兩兩多重比較的FDR控制．上海：上海交通大學(xué)，2015．

［20］Xu H，Gu W，F(xiàn)an J．Control of the False Discovery Rate Under Arbitrary Covariance Dependence．a(chǎn)rXiv：1012．4397，2010．

［21］Fan J，Han X，Gu W．Estimating False Discovery Proportion Under Arbitrary Covariance Dependence．Journal of the American Statistical Association，2012，107（499）：1019-1035．

［22］Leek JT，Storey JD．A general framework formultiple testing dependence．Proceedings of the National Academy of Sciencesof the United States of America，2008，105（48）：18718-18723．

［23］Friguet C，Kloareg M，Causeur D．A Factor Model Approach to Multiple Testing Under Dependence．Journal of the American Statistical Association，2009，104（488）：1406-1415．

［24］Goeman JJ，Solari A．Multiple testing for exploratory research．Statistical Science，2011，26（4）：584-597．

［25］Meinshausen N．False discovery control formultiple tests of association under general dependence．Scandinavian Journal of Statistics，2006，33（2）：227-237．

［26］Hastie T，TibshiraniR，F(xiàn)riedman J．The Elementsof Statistical Learning-Data M ining，Inference and Prediction，2nd edn．New York：Springer，2009，241-245．

［27］王慧．生存分析中半?yún)?shù)模型的變量選擇方法及其模擬研究．太原：山西醫(yī)科大學(xué)，2013．

［28］Arlot S，Celisse A．A survey of cross-validation procedures formodel selection．Statistics Surveys，2010，4（0）：40-79．

［29］Nishii R．Asymptotic Properties of Criteria for Selection of Variables in Multiple Regression．Annals of Statistics，1984，12（2）：758-765．

［30］Fan Y，Tang CY．Tuning parameter selection in high dimensional penalized likelihood．Journal of the Royal Statistical Society，2013，75（3）：531-552．

［31］Akaike H．A New Look at the Statistical Model Identification．Automatic Control IEEE Transactions on，1974，19（6）：716-723．

［32］Schwarz G．Estimating the Dimension of a Model．Annals of Statistics，1978，6（2）：15-18．

［33］Foster DP，George EI．The Risk Inflation Criterion for Multiple Regression．Annals of Statistics，1994，22（4）：1947-1975．

［34］Shao J．An asymptotic theory for linear model selection．Statistica Sinica，1997，7（2）：221-242．

［35］Yang Y．Can the strengths of AIC and BIC be shared？A conflictbetween model indentification and regression estimation．Biometrika，2005，92（4）：937-950．

［36］Wang H，LiR，TsaiCL．Tuning Parameter Selectors for the Smoothly Clipped Absolute Deviation Method．Biometrika，2007，94（3）：553-568．

［37］Ternes N，Rotolo F，M ichiels S．Empirical extensions of the lasso penalty to reduce the false discovery rate in high-dimensional Cox regression models．statistics in medicine，2016，35（15）：2561-2573．

［38］Bogdan M，Ghosh JK，Doerge RW．Modifying the Schwarz Bayesian information criterion to locate multiple interacting quantitative trait loci．Genetics，2004，167（2）：989-999．

［39］Bogdan g，M．？S，Ghosh JK．Selecting explanatory variables w ith themodified version of the Bayesian information criterion．Quality＆Reliability Engineering International，2008，24（6）：627-641．

［40］Chen J，Chen Z．Extended Bayesian information criteria formodel selection w ith largemodel spaces．Biometrika，2008，95（95）：759-771．

［41］Wang T，Zhu L．Consistent tuning parameter selection in high dimensional sparse linear regression．Journal of Multivariate Analysis，2011，102（7）：1141-1151．

［42］Meinshausen N，Bühlmann P．Stability selection．Journal of the Royal Statistical Society，2010，72（4）：417-473．

［43］勾建偉．懲罰回歸方法的研究及其在后全基因關(guān)聯(lián)研究中的應(yīng)用．南京醫(yī)科大學(xué)，2014．

［44］Bühlmann P，Kalisch M，Meier L．High-Dimensional Statisticswith a View Toward Applications in Biology．Annual Review of Statistics＆Its Application，2014，1（1）：255-278．

［45］趙俊琴．基于Lasso的高維數(shù)據(jù)線性回歸模型統(tǒng)計(jì)推斷方法比較．太原：山西醫(yī)科大學(xué)，2015．

［46］Bogdan M，Berg Evd，Su W，et al．Statistical Estimation and Testing via the Ordered L1 Norm．a(chǎn)rXiv：1310．1969，2013．

［47］Bogdan M，van den Berg E，SabattiC，et al．Slope-Adaptive Variable Selection Via Convex Optimization．Ann Appl Stat，2015，9（3）：1103-1140．

［48］Su W，Candes E．SLOPE is Adaptive to Unknown Sparsity and Asymptotically M inimax．Ann Appl Stat，2015，44（3）：1038-1068．

［49］Barber RF，Candès EJ．Controlling the false discovery rate via knockoffs．Annals of Statistics，2014，43（5）：2055-2085．

［50］Reid S，Tibshirani R．Sparse regression and marginal testing using cluster prototypes．Biostatistics，2016，17（2）：364-376．

［51］Barber RF，Candes EJ．A knockoff filter for high-dimensional selective inference．a(chǎn)rXiv：1602．03574，2016．

［52］Verweij PJ，Van Houwelingen HC．Cross-validation in survival analysis．statistics inmedicine，1993，12（24）：2305-2314．

（責(zé)任編輯：郭海強(qiáng)）

國家自然科學(xué)基金項(xiàng)目（81473073）

△通信作者：王彤，E-mail：tongwang＠sxmu．edu．cn