基于一種距離相關的超高維生存數據Model-Free特征篩選

潘瑩麗, 王昊宇, 喻佳麗, 劉展

(應用數學湖北省重點實驗室, 湖北大學數學與統計學學院, 湖北 武漢 430062)

0 引言

隨著時代的進步,科學技術快速發展,高維數據乃至超高維數據頻繁出現在科學研究的各個領域。超高維數據的特征是變量維度p隨樣本量n呈指數階增長,滿足p=O(nα),α>0,其攜帶的信息量一方面讓我們更加全面、客觀、深入地了解研究對象,一方面又大大增加了數據分析的難度。模型中如果涉及的變量過多,不僅會增大模型的預測誤差,還會使模型的算法變得冗余和復雜,從而導致計算時間長、運行速度慢、算法穩定性差等一系列問題,所以建模前對要納入模型的變量進行選擇十分必要。對于協變量維度略高于樣本量的高維數據,傳統的正則化變量選擇方法,如LASSO、SCAD、自適應LASSO和MCP在實際運用中能準確選出對響應變量存在顯著影響的重要變量。但超高維數據復雜程度高,上述正則化方法在處理這類數據時運算速度、統計準確性及算法穩定性均受到極大挑戰,這就迫切要求研究者們提出新的方法來處理超高維數據的降維問題。此外,在生物醫學領域和臨床試驗中,響應變量通常會因為研究對象失聯、研究截止時個體仍然存活等原因無法被完全觀測,成為右刪失的生存數據,處理超高維完全數據的特征篩選法在生存數據中將不再適用,所以尋找一種可以處理超高維右刪失數據的穩定性降維方法顯得尤為重要。

本研究提出一種能有效處理超高維右刪失數據的特征篩選方法,該方法利用距離相關系數刻畫每個協變量和響應變量間的邊際效應,并建立與該系數相關的篩選指標,最后根據確定的篩選準則進行特征篩選。一方面,在有關超高維刪失數據特征篩選方法的研究中,大多篩選方法都是在一些參數或非參數模型的基礎上提出的,由于超高維數據復雜程度高,通常情況下難以確定其真實分布,從海量數據集中指定一個具體模型十分困難,一旦模型出錯,后續的篩選和分析都會受到影響。提出的特征篩選方法不依賴模型假定,可以有效避免模型誤判帶來的不良后果。另一方面,采用的距離協方差估計量是總體距離協方差的一個無偏估計,與其他方法相比計算精度和統計準確性更高。此外,數據刪失或異常值的存在可能會導致篩選結果出現很大偏差,提出的方法可以在響應變量存在刪失和協變量被異常值污染時仍以趨近于1的概率保留所有重要協變量。

1 國內外研究現狀

特征篩選是目前最常見也是最有效處理超高維數據的降維方式。目前國內外已有大量關于特征篩選方面研究成果。Fan和Lv(2008)[1]基于線性模型假定首次提出確定性獨立特征篩選方法SIS,該方法通過每個協變量與響應變量間的皮爾遜相關系數度量協變量的重要程度,在一定的正則條件下該方法滿足確定篩選性。SIS的提出為使用邊際相關系數來區分與響應變量“相關”和“不相關”的協變量奠定了基礎,后續絕大多數篩選方法都借鑒了SIS的思想,并從不同角度對它進行了拓展和改進。如有的學者解除了SIS線性模型的限制,提出了廣義線性模型(Fan和Song,2010)[2]、可加性模型(Fan等,2011)[3]、變系數模型(Fan等,2014)[4]下的獨立特征篩選方法,有的學者則采用了其他邊際相關系數,如距離相關系數(Li等,2012)[5]和Kendall’s相關系數(Song等,2014)[6]。然而,這些特征篩選過程都事先對模型結構作出了假定。超高維數據的結構十分復雜,從這樣復雜、龐大的數據集中找到一個具體模型結構具有挑戰性,一旦模型出錯,后續的篩選結果和分析過程都會受到影響。來鵬等(2018)[7]針對多分類的響應變量,提出了一種改進的基于條件分布的無模型特征篩選方法,該方法對服從重尾分布的協變量有較好的穩健性。針對超高維數據判別分析中的特征篩選問題,何勝美等(2020)[8]在不需要特定模型假設和有限矩條件限制下,提出一種特征篩選方法,該方法基于秩能量矩,適用于多分類因變量情形,用以篩選類別間分布存在明顯差異的協變量。為減少樣本比例不平衡和數據異質性在降維過程中產生的不良影響,何勝美和王響(2021)[9]通過最大化條件分布函數的加權C-V-M距離刻畫協變量對響應變量的邊際效用,提出基于最大邊際效用的SAD-SIS篩選法。

在生物醫學和臨床試驗中,研究對象在實驗期間失聯、實驗截止前提前退出、實驗結束后存活等行為會導致其生存時間無法被準確記錄,成為右刪失的生存數據,在有關超高維生存數據的降維方面,一些學者提出了新的篩選方法。Gorst-Rasmussen和Scheike(2013)[10]基于FAST統計量提出單指標危險率模型下的特征篩選方法,FAST統計量是一個簡單的線性統計量,由每個協變量與協變量平均值的差值乘以生存時間聚合而成,該方法借助FAST統計量刻畫了每個協變量對響應變量的邊際效應?；贔AST統計量的獨立特征篩選方法可視為SIS方法在右刪失數據中的延伸,尤其是在單指標危險率模型的背景下,基于FAST無模型統計量的獨立特征篩選方法滿足SIS的確定篩選性。與完全數據一樣,超高維生存數據中,基于特定模型提出的特征篩選方法受到模型結構影響,其在某些情形下的表現會缺乏穩健性,因此學者們針對生存數據提出一些不依賴模型的特征篩選方法。Song等(2014)[6]提出一種不依賴模型的刪失秩獨立特征篩選方法,該方法中的秩統計量可視為逆概率加權后的Kendall’s相關系數,利用該系數估計值的大小保留與生存時間邊緣秩相關性較強的協變量。在回歸分析中,我們常用相關系數衡量兩個隨機變量間相關性的大小,基于這一思想,Zhang等(2017)[11]根據相關性排序提出一種新的無模型特征篩選方法,該方法通過各協變量與響應變量分布函數的協方差刻畫單個協變量的邊際效應,其不涉及任何非參數逼近和計算,所以計算簡單便捷,運行速度快。然而,這兩種處理刪失數據的變量篩選方法均需采用Kaplan-Meier估計處理右刪失數據,對刪失率較高的情形有一定的弊端。

Székely等(2007)[12]提出距離相關系數的概念,用以度量兩個隨機變量的依賴程度,距離相關系數有兩條優良性質。一是距離相關系數是皮爾遜相關系數絕對值的嚴格增函數,線性模型背景下,若響應變量與協變量均服從正態分布,基于距離相關系數的特征篩選等價于SIS方法。二是兩個隨機變量相互獨立時,它們的距離相關系數為零。這些性質使得它可以用來進行特征篩選。因此,基于Székely等(2007)[12]的距離相關系數,Li等(2012)[5]提出一種不依賴模型結構的特征篩選方法,該方法能很好處理分組協變量或響應變量為多元情形下的特征篩選。Zhang等(2021)[13]提出一種標準化后的基于距離相關的無模型特征篩選,該方法中響應變量由標準化后的觀測時間和刪失示性變量組成,通過計算每個協變量與二元響應變量的距離相關系數,對協變量重要程度進行排序。在響應變量存在刪失或協變量存在異常值時,該方法的表現顯著優于其他方法。然而,上述基于距離相關的特征篩選,距離相關系數估計量均采用了Székely等(2007)[12]的計算方式。本研究基于Székely和Rizzo(2014)[14]提出的距離相關系數估計量提出一種不依賴模型的獨立特征篩選方法,通過非線性、Cox比例風險、轉換模型三種模型和Python中的random包隨機產生滿足設置要求的實驗數據,模擬隨機和非隨機兩種刪失機制和協變量存在異常值的各種情形,并依據模擬結果對不同篩選方法的篩選性能進行比較分析,從而驗證提出方法在變量選擇中的優勢。最后將提出的篩選方法應用于對彌漫大B細胞淋巴瘤患者生存結果有顯著影響的基因篩選之中,進一步說明提出方法在實際應用中的可行性。

2 基于距離相關的Model-Free特征篩選

I={k:S(t|X)依賴于Xk,k=1,…,p}.

Székely等(2007)[12]提出距離協方差的概念,并利用距離協方差測度兩個隨機變量間的獨立性。給定兩個隨機變量x∈Rdx和y∈Rdy,兩者之間的距離協方差dcov(x,y)的平方定義為

(1)

(2)

由式(1)和(2),Székely和Rizzo(2014)[14]給出距離協方差平方dcov2(Xk,Y)的一個估計如下:

(3)

(4)

其中c和κ是已知常數,稱由式(4)定義的篩選方法為基于距離相關(3)排序的獨立篩選法,簡記為bcDC。由Fan和Lv (2008)[1]可知,可以采取如下重要協變量下標估計式的定義方式進行變量篩選,即事先設定一個正整數d0=「n/logn?,

3 數值研究

本節評估基于距離相關的Model-Free篩選方法bcDC在有限樣本下的表現,并將其與現有方法進行比較。為了描述方便,將對比的方法,即Gorst-Rasmussen和Scheike(2013)[10]、Song等(2014)[6]、Zhang等(2017)[11]和Zhang等(2021)[13]提出的方法分別簡記為FAST、CRIS、CR-SIS、CDCS。采用以下5個標準評估不同方法的篩選性能。

1)最小模型大小S:所有重要協變量被包含在選出的模型中所需要保留協變量個數的最小值。

本研究采用500次數值模擬最小模型大小的中位數(Med)和四分位距(IQR)兩個指標進行描述。

2)Pe:對于給定的模型大小「n/logn?,500次數值模擬中每個重要協變量被挑選出來的概率。

3)Pa:對于給定的模型大小「n/logn?,500次數值模擬中所有重要協變量均被挑選出來的概率。

其中,#{·}表示集合中元素個數,?·」med表示序列值的中位數。

顯然,Med的值越接近真實重要協變量集的元素個數,IQR值越小,Pe、Pa和PSR越接近于1,FDR越接近于0,證明對應的特征篩選方法越好。為檢驗bcDC篩選方法的穩健性,不僅考慮隨機刪失和非隨機刪失兩種刪失機制,還考慮某些協變量存在異常值的情形?；诜蔷€性模型、Cox比例風險回歸模型和轉換模型3種不同模型設置產生數據進行模擬。

3.1 非線性模型假設生存時間T產生于一個非線性模型,其與協變量X和隨機誤差項ε滿足下列關系式:

logT=1-5(X2+X3)-3exp{1+10sin(πX1/2)+5X4}+ε,

表1 非線性模型下數值模擬結果

從表1可以看出:在非線性模型背景下,不論是響應變量隨機刪失還是非隨機刪失,bcDC的表現顯著優于FAST和CR-SIS兩種方法,其變量正確選擇率PSR為100%,500次數值模擬結果顯示,bcDC最小模型大小S的中位數為真實模型大小4,四分位距趨近于0,Pa值趨于1,說明不論在上述何種設置下bcDC篩選性都能保持穩定。此外,當協變量存在異常值點時,FAST和CR-SIS方法失效,最小模型大小S的中位數非常大,遠遠超過真實模型大小,Pe和Pa趨近于0,錯誤發現率FDR高達90%以上,而基于距離相關的CDCS和bcDC表現依舊穩健,且bcDC相較于CDCS的FDR更低。

λ(t|X)=λ0(t)exp(XTβ0),

其中基準風險函數λ0(t)=(t-0.5)2,X=(X1,…,Xp)T～Np(0,∑),∑=(0.8|i-j|)(i,j=1,…,p),真實參數β0=(0.35,0.35,0.35,0.35,0.35,0p-5)T,即只有前5個協變量是重要協變量。其他設置與非線性模型設置保持一致,500次數值模擬實驗結果如表2所示。

表2 Cox比例風險模型下數值模擬結果(CR=20%)

表3 Cox比例風險模型下數值模擬結果(CR=40%)

從表2～3可以看出:當協變量不存在異常值時,5種方法均表現良好,變量正確選擇率PSR為100%,除CR-SIS外,其他四種方法每次都能以趨近1的概率保留所有重要協變量,最小模型大小S的中位數近乎接近真實模型大小5。當協變量存在異常值時,基于Cox比例風險模型的FAST和CR-SIS方法篩選性能下降,Pa有時會降至0.3以下,當刪失率為40%時,FDR有時會趨近90%以上,說明這兩種方法在500次數值模擬中絕大多數時刻都無法完整、準確的保留所有重要協變量?；诰嚯x相關的CDCS和bcDC方法的表現則相對穩健,每次都能以95%以上的概率保留所有重要協變量,并且bcDC的錯誤發現率FDR比CDCS方法要低。

H(T)=-β0X+ε,

其中函數H(t)=log(0.5(e2t-1)),協變量X～Np(0,∑),∑=(0.8|i-j|)(i,j=1,…,p),真實參數β0=(-1,-0.9,06,0.8,1.0,0p-10)T。假定隨機誤差項ε是由3種不同分布生成,即標準正態分布,標準柯西分布和標準logistic分布。其他設置與非線性模型設置保持一致,500次數值模擬結果如表4～5所示。

表4 轉換模型下數值模擬結果(CR=20%)

通過觀察表4～5響應變量在不同刪失機制下各方法的篩選結果發現:FAST、CRIS和CR-SIS這3種方法在篩選變量過程中幾乎失效,CDCS和bcDC顯著優于FAST、CRIS和CR-SIS方法。不論響應變量存在20%刪失還是40%刪失,bcDC和CDCS最小模型大小S的中位數均接近真實模型大小4,變量正確選擇率達到100%,并且在500次數值模擬研究中能以99%以上的概率保留所有重要協變量,但bcDC錯誤發現率FDR比CDCS的更低,表現更加穩健。

4 實例分析

為了對比不同篩選方法的表現效果,除考慮bcDC方法外,還考慮Gorst-Rasmussen和Scheike(2013)[10]提出的FAST方法,Song等(2014)[6]提出的CRIS、Zhang 等(2017)[11]提出的CR-SIS和Zhang等(2021)[13]提出的CDCS方法。利用這5種方法從7 399個基因中篩選出[240/log(240)]=43個重要基因。通過bcDC方法篩選出的基因與其他方法篩選出的基因重合率很高,如bcDC方法篩選出的43個重要基因中,有42個基因在不基于模型的CR-SIS和CRIS方法中同樣也被篩選出來,由于表格過大,在此不再展示基因編號。為進一步分析數據,根據bcDC篩選出的協變量擬合Cox比例風險模型,利用正則化方法LASSO、SCAD和MCP從預先選出的43個協變量繼續篩選變量,通過十倍交叉驗證選擇調整參數?；蜻x擇結果和相應回歸系數估計值如表6所示。

表6 基于LASSO、SCAD、MCP重要基因篩選結果

從表6可以看到,43個重要基因中,有6個基因31981、24376、32238、28641、28532、27592被LASSO、SCAD、MCP共同選中,說明這6個基因對病人的生存結果有重大影響。此外,觀察回歸系數估計值可以發現,31981基因對病人生存結果影響最大,這一結果與Li和Luan(2005)[16]和Zhang等(2021)[13]的結論一致,他們同樣將31981基因放在所有對病人生存結果具有重大影響基因的首位。

為進一步評估bcDC方法的預測能力,將240個樣本隨機劃分成n1=160的訓練集和n2=80的測試集,訓練集中的數據用于建立患者生存風險預測模型,測試集中的數據用于模型預測能力的評估。DLBCL生存時間摘要如表7所示。

表7 DLBCL數據集生存時間統計表

本文中首先將bcDC方法用于訓練集,篩選出[160/log(160)]=31個重要基因,然后利用選出的31個基因擬合Cox比例風險模型,計算回歸系數的估計值;在測試集中計算每個患者的風險得分,根據風險得分估計值的平均值將患者分為高、低風險組兩組。圖1描繪了bcDC方法下Kaplan-Meier生存曲線,可以觀察到:低風險組和高風險組的生存曲線可以很好地分離開來。

圖1 生存時間的Kaplan-Meier曲線

為使結果更具有說服力,進一步進行了log-rank檢驗,檢驗兩條生存曲線是否存在顯著差異,bcDC方法log-rank檢驗的p值為0.045<0.05,說明基于bcDC建立的預測模型具有較好的預測能力,可以采用提出的bcDC方法建立預測模型預測未來患者死亡的風險。

5 結論

針對超高維生存數據,基于距離相關系數,提出了一種不依賴模型的特征篩選方法bcDC,該方法不涉及生存時間的Kaplan-Meier估計或復雜的數值優化過程,實施起來較為迅速與便捷。bcDC允許協變量的維度隨著樣本量指數階增長,而且能以較高的概率在保留所有重要變量的前提下快速有效的降低數據的維度。模擬研究表明,在非線性模型、Cox比例風險回歸模型和轉換模型3種不同模型設置下,當響應變量存在刪失或協變量存在異常值時,bcDC方法均能以趨近1的概率保留所有重要協變量并且錯誤發現率也較低。尤其是在非線性模型中,與FAST、CRIS、CR-SIS相比,bcDC能更好地篩選出超高維生存數據的重要變量。實證研究表明bcDC方法具有可實施性,通過該方法篩選出的變量建立的預測模型預測效果較好。