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

(應(yīng)用數(shù)學(xué)湖北省重點實驗室, 湖北大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)學(xué)院, 湖北 武漢 430062)

0 引言

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

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

1 國內(nèi)外研究現(xiàn)狀

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

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

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

2 基于距離相關(guān)的Model-Free特征篩選

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

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

(1)

(2)

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

(3)

(4)

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

3 數(shù)值研究

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

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

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

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

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

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

顯然,Med的值越接近真實重要協(xié)變量集的元素個數(shù),IQR值越小,Pe、Pa和PSR越接近于1,FDR越接近于0,證明對應(yīng)的特征篩選方法越好。為檢驗bcDC篩選方法的穩(wěn)健性,不僅考慮隨機刪失和非隨機刪失兩種刪失機制,還考慮某些協(xié)變量存在異常值的情形。基于非線性模型、Cox比例風(fēng)險回歸模型和轉(zhuǎn)換模型3種不同模型設(shè)置產(chǎn)生數(shù)據(jù)進行模擬。

3.1 非線性模型假設(shè)生存時間T產(chǎn)生于一個非線性模型,其與協(xié)變量X和隨機誤差項ε滿足下列關(guān)系式:

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

表1 非線性模型下數(shù)值模擬結(jié)果

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

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

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

表2 Cox比例風(fēng)險模型下數(shù)值模擬結(jié)果(CR=20%)

表3 Cox比例風(fēng)險模型下數(shù)值模擬結(jié)果(CR=40%)

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

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

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

表4 轉(zhuǎn)換模型下數(shù)值模擬結(jié)果(CR=20%)

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

4 實例分析

為了對比不同篩選方法的表現(xiàn)效果,除考慮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方法中同樣也被篩選出來,由于表格過大,在此不再展示基因編號。為進一步分析數(shù)據(jù),根據(jù)bcDC篩選出的協(xié)變量擬合Cox比例風(fēng)險模型,利用正則化方法LASSO、SCAD和MCP從預(yù)先選出的43個協(xié)變量繼續(xù)篩選變量,通過十倍交叉驗證選擇調(diào)整參數(shù)?；蜻x擇結(jié)果和相應(yīng)回歸系數(shù)估計值如表6所示。

表6 基于LASSO、SCAD、MCP重要基因篩選結(jié)果

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

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

表7 DLBCL數(shù)據(jù)集生存時間統(tǒng)計表

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

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

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

5 結(jié)論

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