高 靜, 陳曉靜, 陳曉林

(曲阜師范大學(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ù)提出了一種確定性獨(dú)立篩選方法(SIS),該方法能將協(xié)變量維數(shù)從超高維快速地降低到一個較小的范圍.這樣就可以方便地使用基于懲罰的變量選擇方法對數(shù)據(jù)進(jìn)行進(jìn)一步分析.Fan和Lv(2008)[1]首次為這種降維方法建立了理論基礎(chǔ).受該工作的啟發(fā),許多統(tǒng)計學(xué)者對這類方法進(jìn)行了研究.Fan和Song(2010)[2]在廣義線性模型下提出了基于最大邊際似然估計的確定性獨(dú)立篩選方法,這實(shí)際上是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ī)變量獨(dú)立等價于它們的距離相關(guān)為零,它能度量變量之間非線性相關(guān)等.因此,它特別適合用來進(jìn)行超高維數(shù)據(jù)的降維.Li等(2012)[5]首先在一般的完全數(shù)據(jù)下研究了基于距離相關(guān)的特征篩選.數(shù)值模擬顯示,該方法明顯改善了SIS方法的不足之處.Huang等(2014)[6]基于Pearson卡方檢驗(yàn)提出了一種特征篩選方法.Pan等(2016)[7]基于期望差提出了一種成對準(zhǔn)確獨(dú)立篩選方法.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)點(diǎn),但需要注意的是Li等(2018)[10]是基于PSH模型這一特定模型的,在實(shí)際應(yīng)用中,此類模型不一定成立.Chen等(2020)[11]和Li(2020)[12]提出的兩種特征篩選方法所選取的篩選指標(biāo)有一定的局限性,此時篩選指標(biāo)等于零并不能說明協(xié)變量與響應(yīng)變量的獨(dú)立性,故不能準(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àn)楦信d趣事件和其它競爭事件.在實(shí)際醫(yī)學(xué)臨床研究中,由于研究時間有限、受訪者失訪等原因,個體的生存時間往往被右刪失了.記刪失時間為C,則實(shí)際觀測時間為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}為一個獨(dú)立同分布樣本.對協(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,可以得到

又因?yàn)镕1(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個例子來檢驗(yàn)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 實(shí)例分析

本節(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)用到膀胱癌實(shí)例中,篩選出[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軟件中的程序包c(diǎn)rrp和cmprsk來實(shí)施基于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)一步驗(yà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é)果驗(yàn)證了本文方法的有限樣本性質(zhì),并表明新提出的CRDC-SIS要優(yōu)于文獻(xiàn)中已經(jīng)有的方法CRSIRS[11]和CRCRS[12].注意到本文的實(shí)際數(shù)據(jù)中包含5個臨床變量,一般來說,這些變量都是對感興趣事件的生存時間有重要影響的變量.在變量篩選的過程中,需要把這樣的先驗(yàn)信息考慮進(jìn)來.因此,超高維競爭風(fēng)險數(shù)據(jù)的條件特征篩選是進(jìn)一步要研究的方向.