超高維生存數(shù)據(jù)中基于相關性秩排序的變量篩選法和FDR控制

潘瑩麗，趙曉洛，張淑瑩，劉 展

（1.湖北大學a.數(shù)學與統(tǒng)計學學院；b.應用數(shù)學湖北省重點實驗室，武漢 430062；2.中南財經(jīng)政法大學 統(tǒng)計與數(shù)學學院，武漢 430073）

0 引言

隨著科學技術的快速發(fā)展，超高維數(shù)據(jù)的降維問題逐漸成為機器學習、計算機科學、統(tǒng)計學等領域的熱點話題。超高維數(shù)據(jù)中往往包含許多冗余信息，在推斷他們與響應變量的統(tǒng)計關系時，通常很難保證統(tǒng)計分析的準確性、模型的可解釋性以及算法的穩(wěn)定性。針對超高維數(shù)據(jù)帶來的挑戰(zhàn)，有效的解決方法之一是先過濾掉一些不相關的協(xié)變量，使維數(shù)從超高維降至相對高維，然后再對剩余的協(xié)變量進行正則化處理來選取重要變量。在生物醫(yī)學和臨床試驗中，經(jīng)常會收集到超高維生存數(shù)據(jù)，但生存數(shù)據(jù)通常存在刪失，傳統(tǒng)的處理超高維完全數(shù)據(jù)的變量篩選法大多不再適用，這就使得研究超高維生存數(shù)據(jù)的變量篩選法成為迫切需求。

已有文獻中關于變量篩選法的研究成果十分豐富。例如，Zhao 和Li（2012）[1]提出了一種基于標準化邊際極大似然估計的Cox 模型變量篩選法，簡稱P-SIS。為了放松Cox 模型的假設條件，Gorst-Rasmussen 和Scheike（2013）[2]基于FAST統(tǒng)計量提出了適用于所有單指標危險率模型的變量篩選法，簡稱FAST-SIS。通過對Kendall’s τ 相關系數(shù)作逆概率加權，Song 等（2014）[3]提出了一種對異常值點穩(wěn)健且不依賴于模型的變量篩選法，簡稱CRIS。張婧（2018）[4]提出了一種基于相關性秩排序的獨立篩選法CR-SIS，該方法以協(xié)方差為切入點來計算每個協(xié)變量對應的邊際效應大小，并以此來衡量協(xié)變量的重要程度，保留與響應變量相關程度高的協(xié)變量。

大多數(shù)變量篩選法往往依賴于一些閾值參數(shù)來控制重要變量和不重要變量之間的界限，并選取與響應變量相關程度較高的協(xié)變量作為重要變量，以達到降維的目的[5]。這些降維方法能以較高的概率保留所有重要變量，即具有確定篩選性能，提高了模型的可解釋性，卻未能很好地控制錯誤發(fā)現(xiàn)率（FDR），降低了穩(wěn)定性，所以若能尋找一種可以平衡模型的可解釋性和穩(wěn)定性的降維方法，則對超高維生存數(shù)據(jù)的處理有著極大的幫助，能科學地利用所獲取到的數(shù)據(jù)，更好地讓試驗數(shù)據(jù)發(fā)揮統(tǒng)計預測的功能。Fixed-X Knockoff 過濾器最早是由Barber 和Candès（2015）[6]提出的，它是一種新的能控制FDR的變量選擇過程，其主要思想是構造原始協(xié)變量的仿制變量（即Knockoff協(xié)變量）。在有條件地基于原始協(xié)變量的前提下，Knockoff協(xié)變量與響應變量無關，并模仿原始變量中的相關結(jié)構，從而實現(xiàn)FDR控制。Fixed-X Knockoff過濾器實現(xiàn)了在數(shù)據(jù)維度小于或等于樣本量的低維情況下對有限樣本的精確FDR控制，保證了有限樣本容量下變量選擇的準確性。

為了解決無法利用Fixed-X Knockoff過濾器對超高維數(shù)據(jù)進行降維的問題，Barber 和Candès（2019）[7]開發(fā)了一個用于檢測高維線性模型中相關性的新框架。在這個框架中，完整樣本被分為兩個不相交的子樣本，樣本大小分別為n1和n2=n-n1。降維過程分兩步進行：第一步，采用變量篩選法將數(shù)據(jù)維數(shù)降至n2以下，并保留與響應變量相關程度高的協(xié)變量；第二步，利用Fixed-X Knockoff 過濾器為選中的每個協(xié)變量構造Knockoff 協(xié)變量，他們能很好地模擬原始協(xié)變量之間的相關結(jié)構，并且可以作為控制變量來衡量原始變量的重要性。綜上，將Fixed-X Knockoff 過濾器應用于數(shù)據(jù)分析的效果有一定的理論支撐。鑒于此，本文將超高維生存數(shù)據(jù)的變量篩選法和Fixed-X Knockoff 過濾器結(jié)合起來（簡記為CR-Knockoff）解決超高維生存數(shù)據(jù)的降維問題，以確保模型的可解釋性和穩(wěn)定性，拓寬降維方法的實際應用范圍。若將完全數(shù)據(jù)視為刪失率為0 的生存數(shù)據(jù)，則該過程同樣也適用于完全數(shù)據(jù)的降維。此外，本文提出的CR-Knockoff 兩步過程同時具有確定篩選和FDR控制的性能。

1 基于相關性秩排序的變量篩選法和FDR控制

1.1 CR-Screen過程

令yˉ為生存時間，C為刪失時間，y=min(,C)表示實際觀測的生存時間，Δ=I(≤C)為刪失指示變量，I(?)為示性函數(shù)。假設數(shù)據(jù)是完全隨機刪失的，即給定p維協(xié)變量x=(x1,…,xp)，刪失時間C和生存時間是相互獨立的。假定對于給定的協(xié)變量x，的條件分布函數(shù)為：

在不指定任何模型的假設前提下，定義重要協(xié)變量指標集為：

重要變量的數(shù)量為s=|Λ |，其中，| Λ |表示集合Λ 中元素的個數(shù)。用Λ 的補集Λc來表示不重要變量指標集。假設所有協(xié)變量均已中心化，即對j=1,…,p均有E(xj)=0，若協(xié)變量不滿足這個假設條件，則可以先對其進行中心化處理。令G(h)=P(≤h)為響應變量的分布函數(shù)，定義R(x,)=E{xG()}。令R(xj,)是R(x,)的第j個元素，則：

其中，xj為x的第j個元素。對j=1,…,p,定義：

其中，CR()xj,代表第j個協(xié)變量xj所對應的邊際效應。

給定來自總體{y,x,Δ} 的隨機樣本{Y,X,Δ}，其中，Y=(Y1,…,Yn)T，設計矩陣X=(X1,…,Xp)且Xj=(X1j,…,Xnj)T，j=1,…,p,Δ=(Δ1,…,Δn)T。由于響應變量為刪失數(shù)據(jù)，則可以用Kaplan-Meier 估計來得到生存函數(shù)G(h)的估計值，其具體定義如下：

其中，cn-φ為常數(shù)閾值。

1.2 基于Knockoff協(xié)變量的FDR控制過程

Fixed-X Knockoff過濾器最早是由Barber 和Candès（2015）[6]針對確定性設計矩陣提出的，它是一種新的帶有FDR 控制的變量選擇方法。令為響應變量的n次觀測向量，x=(x1,…,xp)表示p維協(xié)變量。

本文采用SDP方法[6]來構造p維協(xié)變量的Knockoff協(xié)變量，基于{x,,} 定義如下統(tǒng)計量：

參考Liu 等（2020）[8]的研究思路，定義如下閾值參數(shù)Tα：

1.3 CR-Knockoff過程

采用簡單樣本分割思想，將完整樣本隨機拆分為兩個不重疊的、樣本量分別為n1和n2=n-n1的子樣本集。不妨令是X的隨機拆分，并且對Y進行相同的拆分，則可以將X和Y分別寫成：

CR-Knockoff過程的兩個步驟如下：

（1）CR-Screen步驟

（2）Knockoff步驟

其中，a?b=max{a,b} ?；谥匾兞恐笜思?Tα)，功效（Power）可定義為：

2 數(shù)值模擬

2.1 CR-Screen篩選性能

本文通過數(shù)值模擬實驗評估CR-Screen 的有限樣本性能，并將其與另外三種處理超高維生存數(shù)據(jù)的變量篩選法進行比較，他們分別是Zhao和Li（2012）[1]提出的一種基于標準化邊際極大似然估計的Cox模型變量篩選法，簡稱P-SIS；Gorst-Rasmussen和Scheike（2013）[2]基于FAST統(tǒng)計量提出的適用于所有單指標危險率模型的變量篩選法，簡稱FAST-SIS；Song等（2014）[3]提出的一種對異常值點穩(wěn)健且不依賴于模型的變量篩選法，簡稱CRIS?？紤]如下三種模型：

模型1：假設生存時間yˉ服從如下線性模型：

其中，協(xié)變量x=(x1,…,xp)由均值為0、協(xié)方差矩陣為Σ=(0.5|i-j|)(i,j=1,…,p)的多元正態(tài)分布生成；回歸參數(shù)，即只有前5 個協(xié)變量為重要變量，其他p-5 個協(xié)變量均為不重要變量；誤差項ε服從標準正態(tài)分布；刪失時間C服從參數(shù)為δ的指數(shù)分布，通過調(diào)整δ使刪失率接近20%?？紤]樣本量n=100、協(xié)變量的維數(shù)p=5000 的情況。

模型2：假設生存時間yˉ服從如下Cox比例風險模型：

模型3：假設生存時間yˉ服從如下非線性模型：

其中，協(xié)變量x、誤差項和刪失時間的設置同模型1。

在數(shù)值模擬實驗中,為了比較各種方法的效果，采用三個評價指標進行評價，定義如下：

（1）最小模型大小。通過200 次模擬實驗所得到?值的5%、25%、50%、75%、95%分位數(shù)來衡量此變量篩選法所選出模型的復雜程度。

（2）Pe:給定模型大小在200 次模擬實驗中每個重要變量單獨被選出來的概率，其中[x] 代表x的整數(shù)部分。

（3）Pa:給定模型大小在200 次模擬實驗中所有重要變量同時被選出來的概率。

表1匯總了200次模擬實驗中最小模型大小? 的5%、25%、50%、75%、95%分位數(shù)。

表1 200次模擬實驗中最小模型大小? 的5%、25%、50%、75%、95%分位數(shù)

由表1可知，CR-Screen在模型1和模型3中的表現(xiàn)均優(yōu)于其他三種變量篩選法。就模型1 而言，CR-Screen、P-SIS和FAST-SIS在? 的5%、25%、50%和75%分位數(shù)處均為真實模型大?。?），但是CR-Screen 和P-SIS 在? 的95%分位數(shù)處表現(xiàn)較好。在模型2 中，變量篩選法P-SIS和FAST-SIS 的表現(xiàn)較好，CR-Screen、P-SIS 和FAST-SIS在? 的5%、25%和50%分位數(shù)處均為5，但是CR-Screen在?的75%和95%分位數(shù)處的表現(xiàn)不如P-SIS和FAST-SIS，主要原因是P-SIS 和FAST-SIS 是基于Cox 比例風險模型的某些特征和特殊結(jié)構提出的。觀察模型3的結(jié)果可知，CR-Screen的表現(xiàn)優(yōu)于另外三種方法。

表2匯總了200次模擬實驗中每個重要變量單獨被選出的概率Pe和所有重要變量同時被選出的概率Pa。

表2 200次模擬實驗中每個重要變量單獨被選出的概率Pe 和所有重要變量同時被選出的概率Pa

由表2 可知，在模型1 中，CR-Screen 的確定篩選性能與P-SIS、FAST-SIS相差不大，均能以接近于1的概率選出所有重要變量。在模型2中，CR-Screen、P-SIS和FAST-SIS均能以接近于1的概率選出所有重要變量，這三者的表現(xiàn)明顯優(yōu)于CRIS。就模型3的結(jié)果而言，與其他三種方法相比，CR-Screen 能夠以0.895 的概率保留所有重要變量，具有較強的確定篩選性能。

2.2 CR-Knockoff篩選和FDR控制性能

本文通過數(shù)值模擬實驗評估CR-Knockoff 過程的變量篩選性能和FDR 控制性能?？紤]三個不同的模型，具體為：

模型6：假設生存時間yˉ來自如下模型：

其中，所有剩余設置同模型3。此外，樣本量、參數(shù)維數(shù)和樣本的分割的設置同模型4。

為了評估CR-Knockoff的變量篩選和FDR控制性能，定義如下6個評價指標：

（1）Size:平均模型大小,即200 次模擬實驗中選出重要變量個數(shù)的平均值。

（2）Pe:通過Knockoff 過程確定閾值Tα，從而選出滿足條件的協(xié)變量，該指標表示在200次模擬實驗中每個重要變量分別被選出來的概率。

（3）Pa:通過Knockoff 過程確定閾值Tα，從而選出滿足條件的協(xié)變量，該指標表示在200次模擬實驗中所有重要變量同時被選出來的概率。

（4）α：預先指定的FDR 控制水平，考慮α=0.10，0.15，0.20，0.25，0.30的情況。

（5）FDR：CR-Knockoff 在變量篩選過程中的錯誤發(fā)現(xiàn)率。

（6）Power：功效，其用于衡量降維過程CR-Knockoff的可靠程度。

表3匯總了基于CR-Knockoff過程進行變量選擇的結(jié)果。

表3 基于CR-Knockoff過程的變量選擇結(jié)果

由表3 可知，對于模型4，每個變量被選出的概率Pe均大于或等于0.93，設定的10 個重要變量同時被選出的概率Pa也大于0.91，且隨著α的增大而增大。對于模型5，所有重要變量同時被選出的概率Pa大于或等于0.87。對于模型4，當α≥0.20 時，CR-Knockoff 過程選出所有重要變量的概率Pa達到0.8?？傊?，CR-Knockoff 過程具有較好的確定篩選性能。

表4 匯總了基于CR-Knockoff過程的FDR控制結(jié)果。

表4 基于CR-Knockoff過程的FDR控制結(jié)果

由表4可知，基于模型4、模型5和模型6這三個模型，CR-Knockoff過程均能夠?qū)㈠e誤發(fā)現(xiàn)率控制在預先設定的錯誤發(fā)現(xiàn)率水平α下且比較接近α的取值。由模型4 中的錯誤發(fā)現(xiàn)率控制結(jié)果可知，功效值均大于0.93，并隨著α的增大而增大。在模型5 和模型6 中，功效值均大于0.83，并且模型6中的功效值隨著α的增大而增大。

參見Liu 等（2020）[8]的研究中的定理5，當α≥1時，CR-Knockoff 能夠以較高概率保留所有重要變量，并實現(xiàn)FDR 控制。相反，當α<1時，CR-Knockoff 過程卻不能篩選出所有重要變量。為了進一步驗證該定理，考慮前文中相同設置的線性回歸模型、Cox 比例風險模型和非線性模型，設定真實重要變量數(shù)量為s=10，并讓FDR 水平α以0.01 為間隔從0.01 增加到0.20。對每個設定的水平α，計算200 次模擬下CR-Knockoff 的確定性篩選概率。為了更直觀地展示結(jié)果，在圖1 中繪制出確定性篩選概率與預先指定的FDR水平α的關系圖。垂直的黑色點線為截斷點α=1。從圖1 中可以清楚地看出，當α≥1時，所有重要變量都有較高的概率被選中；當α<1時，確定篩選性能并不成立，并且當α很小時確定性篩選概率接近于零。

圖1 確定性篩選概率與預先指定的FDR水平α 的關系曲線

3 實證分析

為了驗證本文所提出的方法在實際應用中的可操作性，將提出的CR-Knockoff 過程應用于套細胞淋巴瘤（MCL）數(shù)據(jù)集，該基因表達數(shù)據(jù)集包含了8810個基因，該數(shù)據(jù)可在網(wǎng)頁http://llmpp.nih.gov/MCL/上獲得。本文的主要目標是確定對患者生存風險有重大影響的基因。在101例從未有過淋巴瘤病史的未經(jīng)治療的患者中，基于形態(tài)學和免疫表型標準，有92 例患者被歸類為患有MCL。在實驗的整個隨訪過程中，有64例患者死于MCL，剩下的28 例患者在最后一次隨訪時仍然存活，該數(shù)據(jù)集的刪失率為30.4%。將患者的生存時間看作響應變量，剔除含缺失值的基因，將剩下的6312 個基因表達水平看作協(xié)變量。在進行變量篩選之前，對所有基因表達水平進行標準化處理。

根據(jù)前文所提出的CR-Knockoff 過程的實現(xiàn)步驟，設定CR-Screen 篩選步驟的樣本量為42，d=[42/log(42)]=11，即先從92個個體中隨機選擇42個觀察對象，在篩選步驟中使用CR-Screen 預篩選出11 個基因；然后再在Knockoff 步中利用剩下的50 個觀察對象為預先篩選出的11 個協(xié)變量構造近似二階Knockoff 協(xié)變量。在預先指定的FDR 水平α=0.2 下，CR-Knockoff 過程選出了5 個重要變量，對所選出的5個重要變量擬合Cox比例風險模型，以獲得重要變量的系數(shù)估計值。為了進行比較，考慮Zhao和Li（2012）[1]中的篩選方法P-SIS、Gorst-Rasmussen 和Scheike（2013）[2]中的篩選方法FAST-SIS、Song 等（2014）[3]中的篩選方法CRIS、張婧（2018）[4]中的篩選方法CR-Screen。分別采用這四種變量篩選法從6312 個基因中篩選出20=[92/log(92)]個協(xié)變量。為了進一步分析這一數(shù)據(jù)，擬合基于20個所選協(xié)變量的加LASSO懲罰的Cox比例風險模型，通過10 折交叉驗證選擇優(yōu)化參數(shù)。加LASSO處理后的CRIS（簡記為CRIS+LASSO）、P-SIS（簡記為P-SIS+LASSO）、FAST-SIS（簡記為FAST-SIS+LASSO）、CR-Screen（簡記為CR-Screen+LASSO），分別選出5、5、6、6個重要變量，基因的編號（UNIQID）和選出的重要變量的系數(shù)估計值的結(jié)果總結(jié)在表5中。

表5 基于MCL數(shù)據(jù)集篩選重要基因的結(jié)果

由表5 可以看出，CR-Knockoff 方法和CR-Screen+LASSO 方法篩選出兩個共同的基因30157 和15936；CR-Knockoff 方法和FAST-SIS+LASSO 方法以及P-SIS+LASSO方法均篩選出一個共同的基因30157。此外，Zhang等（2017）[9]指出基因30157與患者的生存率密切相關。因此，基因30157是對患者生存風險有重大影響的基因。

為了評估CR-Knockoff的預測性能，將92名患者隨機分為樣本量為72 的訓練集和樣本量為20 的測試集，利用訓練集中的數(shù)據(jù)建立預測模型，并利用測試數(shù)據(jù)對模型進行評價。將CR-Knockoff 過程的FDR 控制在α=0.2 水平下，對訓練集進行篩選，設定CR-Screen 篩選步驟的樣本量為32，d=[32/log(32)]=9。CR-Knockoff 過程選出了5個重要基因，他們的UNIQIDs分別為26944、17474、16312、28640、27116。基于所選的基因擬合Cox 比例風險模型，并獲得協(xié)變量的系數(shù)估計值。計算測試集中患者的風險得分，并將他們分為低風險和高風險組，分界值由訓練集中患者風險得分的平均值確定。下頁圖2 繪制了測試集中低風險和高風險組的Kaplan-Meier 生存曲線，可以看出，對于所提出的CR-Knockoff 過程，兩條曲線可以很好地分開。接著，采用對數(shù)秩檢驗來評估這兩條生存曲線的差異，檢驗的P 值為4×10-3，表明基于CR-Knockoff 過程最終選擇的模型具有良好的預測性能。

圖2 基于測試集所得低風險和高風險組的Kaplan-Meier生存曲線

4 結(jié)論

本文將變量篩選法CR-Screen 和Fixed-X Knockoff 過濾器結(jié)合起來解決超高維生存數(shù)據(jù)的降維問題，以達到平衡模型的可解釋性和穩(wěn)定性的目的，拓寬了降維方法的實際應用范圍。本文提出使用Knockoff 協(xié)變量來指定變量篩選閾值的兩步CR-Knockoff過程，在預先設定的FDR水平α大于或等于1的情況下，提出的兩步過程同時具有確定篩選和FDR 控制的性能。模擬和實證研究結(jié)果表明，CR-Knockoff 過程可以很好地從大量超高維生存信息中找出與所感興趣的變量真正相關的協(xié)變量，并能有效控制FDR。