基于距離相關(guān)的超高維競爭風(fēng)險數(shù)據(jù)特征篩選*

高 靜, 陳曉靜, 陳曉林

(曲阜師范大學(xué)統(tǒng)計與數(shù)據(jù)科學(xué)學(xué)院,273165,山東省曲阜市)

0 引 言

隨著社會經(jīng)濟(jì)和科學(xué)技術(shù)的發(fā)展與進(jìn)步,社會衛(wèi)生經(jīng)濟(jì)學(xué)、生物信息學(xué)、金融學(xué)等各個領(lǐng)域都出現(xiàn)了復(fù)雜的超高維數(shù)據(jù).在超高維數(shù)據(jù)中,協(xié)變量維數(shù)p遠(yuǎn)遠(yuǎn)大于樣本量n.特征篩選作為一種處理超高維數(shù)據(jù)的重要統(tǒng)計方法,近年來引起了研究學(xué)者們的廣泛關(guān)注.Fan和Lv(2008)[1]在超高維線性模型下基于Pearson相關(guān)系數(shù)提出了一種確定性獨立篩選方法(SIS),該方法能將協(xié)變量維數(shù)從超高維快速地降低到一個較小的范圍.這樣就可以方便地使用基于懲罰的變量選擇方法對數(shù)據(jù)進(jìn)行進(jìn)一步分析.Fan和Lv(2008)[1]首次為這種降維方法建立了理論基礎(chǔ).受該工作的啟發(fā),許多統(tǒng)計學(xué)者對這類方法進(jìn)行了研究.Fan和Song(2010)[2]在廣義線性模型下提出了基于最大邊際似然估計的確定性獨立篩選方法,這實際上是SIS方法的一個拓展.Liu等(2014)[3]在變系數(shù)模型下提出了基于條件相關(guān)系數(shù)的特征篩選方法.Zhu等(2011)[4]提出了一種無模型的特征篩選方法,適用于許多常用的參數(shù)和半?yún)?shù)模型.在度量變量之間相關(guān)性的指標(biāo)中,距離相關(guān)具有一些優(yōu)良的性質(zhì),例如兩個隨機(jī)變量獨立等價于它們的距離相關(guān)為零,它能度量變量之間非線性相關(guān)等.因此,它特別適合用來進(jìn)行超高維數(shù)據(jù)的降維.Li等(2012)[5]首先在一般的完全數(shù)據(jù)下研究了基于距離相關(guān)的特征篩選.數(shù)值模擬顯示,該方法明顯改善了SIS方法的不足之處.Huang等(2014)[6]基于Pearson卡方檢驗提出了一種特征篩選方法.Pan等(2016)[7]基于期望差提出了一種成對準(zhǔn)確獨立篩選方法.Kong等(2019)[8]基于復(fù)合決定系數(shù)提出了一種篩選方法.Zhou等(2020)[9]基于累積差分提出了一種無模型向前篩選方法.

超高維生存數(shù)據(jù)是一類涉及生存時間的特殊類型的超高維數(shù)據(jù).在這一類數(shù)據(jù)中,響應(yīng)變量是研究個體的生存時間,例如病人從確診某種疾病到死亡的時間,某人從失業(yè)到再就業(yè)的時間等等.在生物醫(yī)學(xué)研究中,研究個體可能不僅會經(jīng)歷感興趣疾病的死亡,還會經(jīng)歷其它不感興趣疾病的死亡.例如,某研究關(guān)心病人從確診到死于癌癥的時間,但是病人還可能會死于心臟病、糖尿病等其它疾病.這種情況下,病人從確診到死于癌癥的時間是不一定能觀測到的.在文獻(xiàn)中,把關(guān)心的疾病導(dǎo)致的死亡稱為感興趣事件,其它疾病或原因造成的死亡稱為競爭事件.這類涉及競爭事件的生存數(shù)據(jù)被稱為競爭風(fēng)險數(shù)據(jù).不同于標(biāo)準(zhǔn)的生存數(shù)據(jù),競爭風(fēng)險數(shù)據(jù)中不僅會出現(xiàn)由于研究終止、研究個體退出研究等原因?qū)е碌母信d趣事件觀測不到(即右刪失),還會出現(xiàn)由于競爭事件的發(fā)生導(dǎo)致感興趣事件觀測不到.在分析競爭風(fēng)險數(shù)據(jù)時,既不能刪除右刪失的數(shù)據(jù),也不能簡單地把競爭事件的發(fā)生認(rèn)為是右刪失.因此,超高維標(biāo)準(zhǔn)生存數(shù)據(jù)的特征篩選方法不能簡單地用于超高維競爭風(fēng)險數(shù)據(jù).對于超高維競爭風(fēng)險數(shù)據(jù)的特征篩選問題,文獻(xiàn)中的研究還較少.Li等(2018)[10]在比例子分布風(fēng)險(proportional subdistribution hazards,PSH)模型下提出了兩階段的方法,并證明了相應(yīng)方法的確定性篩選性質(zhì).Chen等(2020)[11]把Zhu等(2011)[4]的方法推廣到了超高維競爭風(fēng)險數(shù)據(jù),提出了一種無模型的特征篩選方法(CRSIRS).Li(2020)[12]基于相關(guān)秩提出了一種特征篩選方法(CRCRS).盡管這些方法具有很多優(yōu)點,但需要注意的是Li等(2018)[10]是基于PSH模型這一特定模型的,在實際應(yīng)用中,此類模型不一定成立.Chen等(2020)[11]和Li(2020)[12]提出的兩種特征篩選方法所選取的篩選指標(biāo)有一定的局限性,此時篩選指標(biāo)等于零并不能說明協(xié)變量與響應(yīng)變量的獨立性,故不能準(zhǔn)確探測兩者之間的非線性關(guān)系.為了改進(jìn)這些方法存在的不足,提出更加適合超高維競爭風(fēng)險數(shù)據(jù)的特征篩選方法就顯得尤為重要.

針對超高維競爭風(fēng)險數(shù)據(jù),在稀疏性假定下,本文提出了一種新的基于距離相關(guān)的特征篩選方法.該方法對于協(xié)變量、感興趣事件時間、其他競爭事件時間都具有很好的穩(wěn)健性.此外,該方法不需要預(yù)設(shè)某種特定的回歸模型,可應(yīng)用于協(xié)變量、感興趣事件時間、其他競爭事件時間之間的各種線性、非線性關(guān)系.

1 方 法

記個體的生存時間為T,J為失效類型,J的取值為1和2,分別代表失效原因為感興趣事件和其它競爭事件.在實際醫(yī)學(xué)臨床研究中,由于研究時間有限、受訪者失訪等原因,個體的生存時間往往被右刪失了.記刪失時間為C,則實際觀測時間為Y=T∧C,刪失指標(biāo)為δ=I(T≤C).在這里,需要注意的是當(dāng)生存時間T被C刪失時,失效類型J也就未知了.記對感興趣事件時間有潛在影響的協(xié)變量為x=(X1,X2,…,Xp)T,其中p是維數(shù),對響應(yīng)變量有重要影響的協(xié)變量的指標(biāo)集合為

M={1≤k≤p:F1(y|x)依賴于Xk},

其中F1(y|x)為給定協(xié)變量x時感興趣事件時間的條件累積發(fā)病率函數(shù).我們的目的是把它估計出來.

對于競爭風(fēng)險數(shù)據(jù),不能直接利用Székely等(2007)[13]提出的距離相關(guān)來衡量協(xié)變量與感興趣事件時間的相關(guān)程度.為了能使用距離相關(guān)度量協(xié)變量與感興趣事件時間的相關(guān)性,我們通過它們的分布函數(shù)和累積發(fā)病率函數(shù)對它們進(jìn)行變換.對每個協(xié)變量Xk,我們記其累積分布函數(shù)為Fk(xk)=Pr(Xk≤xk)；對感興趣事件時間,我們記其累積發(fā)病率函數(shù)為F1(y)=Pr(T≤y,J=1).為了度量每個協(xié)變量Xk和感興趣事件時間的相關(guān)性,考慮Fk(Xk)和F1(T)的距離相關(guān)

其中dcov(Fk(Xk),F1(T))，dcov(Fk(Xk),Fk(Xk))和dcov(F1(T),F1(T))三者的定義類似,下面以dcov(Fk(Xk),F1(T))為例說明其表達(dá)式的具體形式,其余兩者的表達(dá)式可以相應(yīng)得出.令

下面考慮特征篩選指標(biāo)的估計. 記{Yi,δi,δiJi,xi,i=1,…,n}為一個獨立同分布樣本.對協(xié)變量Xk,累積分布函數(shù)Fk(xk)的估計為

感興趣事件時間累積發(fā)病率函數(shù)的估計為

為了證明CRDC-SIS的大樣本性質(zhì),我們給出下面的假設(shè)：

(A.1) 存在正常數(shù)τ和η,使得Pr(t≤T≤C)≥η,其中t∈(0,τ],τ是最大隨訪時間；sup{t:Pr(T>t)>0}≥sup{t:Pr(C>t)>0}；函數(shù)G(y)具有一致有界的一階導(dǎo)數(shù).

(A.2) 存在正常數(shù)c1和κ∈(0,1/2),使得真正重要變量篩選指標(biāo)的最小值滿足

假設(shè)(A.1)在生存分析文獻(xiàn)中很常見,可以保證Kaplan-Meier估計有良好的表現(xiàn).許多文獻(xiàn)中都采用了該假設(shè),例如Zhou等(2017)[14]和Chen等(2020)[11].假設(shè)(A.2)要求篩選指標(biāo)的最小值不能太小,即真正重要變量與響應(yīng)變量之間的相關(guān)性不能太弱.

為了證明CRDC-SIS的確定性篩選性質(zhì),首先給出一個引理.記概率空間(Ω,Fz,Pr),其中Ω為樣本空間,F為σ代數(shù),Pr為概率.記

根據(jù)假設(shè)(A.1)和Zhou等(2017)[14]的引理1,可以得到

是幾乎處處成立的.記上式成立的集合為Ω0,則Pr(Ω0)=1.以下的證明都是在Ω0上進(jìn)行的.

引理1假設(shè)(A.1)下,對于任意ε∈(0,1),可得

該引理可由Li(2020)[12]的引理1得到.

在給出CRDC-SIS確定性篩選性質(zhì)之前,再引入新的記號.記

定理1假設(shè)(A.1)下,有

其中d1為一正常數(shù).進(jìn)一步地,假設(shè)(A.2)也成立,令vn=c1n-κ,可以得到

其中q是協(xié)變量中真正重要的變量的個數(shù).

?T1+T2.

對于T1,可以得到

又因為F1(Yi)-F1(Yj)|≤1,有

其中最后一個不等式由DKW不等式得到.

對于T2,可以得到

其中最后一個不等式由引理1得到. 因此

2exp(-nε2/c2)+2(n+1)exp(-nε2c3)=

其中c2=1/8,c3=1/32,d1=min{c2,c3}.

又根據(jù)Li等(2012)[5]的證明可知

最終得到

O{pexp(-d1n1-2κ)}.

因此可以得到

1-O{qexp(-d1n1-2κ)}.

定理證畢.

2 數(shù)值模擬

在本節(jié),我們通過2個例子來檢驗CRDC-SIS的有限樣本性質(zhì),并和CRSIRS[11]和CRCRS[12]2種篩選方法進(jìn)行比較.對于協(xié)變量x=(X1,X2,…,Xp)T,我們假定x～N(0p×1,Σ=(ρ|i-j|)p×p),ρ=0.9.在模擬中,樣本量n取為200,協(xié)變量維數(shù)p的取值為1000和3000,刪失時間C服從(0,6)上的均勻分布.基于500次的數(shù)據(jù)重復(fù),我們考慮下面的指標(biāo)：真正重要變量被篩選出來的頻率分別記為pk；全部重要變量都被篩選出的概率Pa.

例1假設(shè)感興趣事件發(fā)生的概率為Pr(J=1)=0.7,競爭事件發(fā)生的概率為Pr(J=2)=1-Pr(J=1)=0.3.感興趣事件時間和競爭事件時間分別為T1和T2,兩者分別滿足回歸模型

log(T1)=0.8cos(X1+0.5X2+X3)+sin(X2)+cos(1.5X3)+ε1

和

log(T2)=0.8sin(1.5X1+0.25)+sin(X2X3)+0.5tan(X3)+ε2,

其中隨機(jī)誤差ε1和ε2服從標(biāo)準(zhǔn)正態(tài)分布.從模型中,可以看出,真正重要的變量為X1,X2和X3.在本例中,刪失率大約為45%.例1的模擬結(jié)果列于表1.

表1 例1中pk和pa的模擬結(jié)果

根據(jù)表1,可以看出,在不同參數(shù)設(shè)定下,本文所提出的篩選方法將所有重要變量篩選出來的頻率pa要明顯大于CRSIRS[11]和CRCRS[12]2種方法.

例2本例中感興趣事件和其他競爭事件發(fā)生的概率同例1.T1和T2分別滿足的回歸模型如下

log(T1)=2sin(2X1+X3)+0.5tan(X1X2)+1.5cos(X2)+1.5sin(X2+X3)+ε1

和

log(T2)=0.5sin(X1+0.15X2)+0.2sin(X1)+X3+0.15tan(X3)+ε2,

其中隨機(jī)誤差ε1和ε2服從標(biāo)準(zhǔn)正態(tài)分布.真正重要的變量為X1,X2和X3.在本例中,刪失率大約為46%.例2的模擬結(jié)果列于表2.

表2 例2中pk和pa的模擬結(jié)果

我們可以從表2中看出CRDC-SIS方法一致地優(yōu)于另外2種方法.此外,根據(jù)表2發(fā)現(xiàn),在非線性模型中,CRSIRS[11]方法將各個重要變量篩選出來的頻率p1,p2和p3都非常小,在此情況下,將所有真正重要變量全部篩選出來的頻率pa也很小.對于CRCRS[12]方法來說,雖然頻率較CRSIRS[11]方法大,但其表現(xiàn)稍差于本文的篩選方法.因此,本文所提出的篩選方法對于本例的非線性模型也有很好的表現(xiàn).

3 實例分析

本節(jié)將CRDC-SIS方法應(yīng)用到膀胱癌數(shù)據(jù)集[15].該數(shù)據(jù)集包含301位病人的數(shù)據(jù),其中響應(yīng)變量數(shù)據(jù)分為3類：患者死于膀胱癌,即感興趣事件；死于其他原因,即其他競爭事件；由于臨床研究時間限制或者失訪而產(chǎn)生的刪失數(shù)據(jù).除了生存時間,本數(shù)據(jù)集還包含1381個基因和5個臨床變量.據(jù)了解,人體內(nèi)某些基因會影響膀胱癌的發(fā)病幾率.因此,我們擬利用本文提出的方法篩選出對膀胱癌有重要影響的協(xié)變量.首先,我們將CRDC-SIS、CRSIRS[11]和CRCRS[12]3種篩選方法應(yīng)用到膀胱癌實例中,篩選出[n/log(n)]=52個基因. 3種方法共同篩選出了16個基因.這些基因是更有可能導(dǎo)致疾病發(fā)生的重要基因.

由于篩選出來的基因個數(shù)還較多,我們采用基于懲罰的變量選擇方法進(jìn)行進(jìn)一步的降維.在競爭風(fēng)險數(shù)據(jù)分析中,應(yīng)用最廣泛的回歸模型是Fine和Gray(1999)[16]提出的PSH模型.借助R軟件中的程序包crrp和cmprsk來實施基于PSH模型的變量選擇方法.基于各種方法的變量選擇結(jié)果列于表3,其中Est.和Se.分別表示回歸系數(shù)估計和標(biāo)準(zhǔn)誤差估計.回歸系數(shù)估計越大,標(biāo)準(zhǔn)誤差估計越小,則說明識別出的基因?qū)τ诎螂装┑挠绊懺斤@著.

從下頁表3看出,CRDC-SIS方法下,SEQ820、SEQ279、SEQ377和SEQ833在4種變量選擇方法下都被識別出來.同時,SEQ820、SEQ279在CRCRS[12]方法下也被4種方法所選擇,說明這兩個基因極有可能是導(dǎo)致膀胱癌的關(guān)鍵原因.CRSIRS[11]方法下除ALASSO外,其余3種變量選擇方法都將SEQ820選擇出來,這進(jìn)一步驗證了SEQ820對于膀胱癌的重要影響.此外,在同一懲罰函數(shù)下,CRDC-SIS方法所確定的基因個數(shù)少于其余兩種,根據(jù)回歸系數(shù)的估計可知CRDC-SIS方法選出的基因更為顯著.這進(jìn)一步說明了本文所提方法的優(yōu)越性.

表3 不同懲罰的變量選擇方法下的結(jié)果

4 結(jié) 論

基于Székely等(2007)[13]提出的距離相關(guān),本文構(gòu)建了一種適用于超高維競爭風(fēng)險數(shù)據(jù)的特征篩選方法CRDC-SIS，從理論和數(shù)值模擬兩方面對其進(jìn)行了詳細(xì)的研究.理論研究表明在很弱的假設(shè)下,CRDC-SIS具有確定性篩選性質(zhì).數(shù)值結(jié)果驗證了本文方法的有限樣本性質(zhì),并表明新提出的CRDC-SIS要優(yōu)于文獻(xiàn)中已經(jīng)有的方法CRSIRS[11]和CRCRS[12].注意到本文的實際數(shù)據(jù)中包含5個臨床變量,一般來說,這些變量都是對感興趣事件的生存時間有重要影響的變量.在變量篩選的過程中,需要把這樣的先驗信息考慮進(jìn)來.因此,超高維競爭風(fēng)險數(shù)據(jù)的條件特征篩選是進(jìn)一步要研究的方向.