Cox比例風(fēng)險模型的彈性SCAD變量選擇方法及應(yīng)用*

孫志華，謝向東，焦東方

(1.中國海洋大學(xué)數(shù)學(xué)科學(xué)學(xué)院，山東 青島 266100；2.東北師范大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，吉林 長春 130024)

變量選擇問題是統(tǒng)計領(lǐng)域中的研究方向之一。基于懲罰思想的變量選擇方法(也稱正則化回歸方法)，在選出變量的同時也對變量參數(shù)作出了估計，其計算量相對較小，相比其他變量選擇方法呈現(xiàn)出諸多優(yōu)越性。這使得，以Lasso為代表的基于懲罰的變量選擇方法得到廣泛研究，出現(xiàn)了一系列基于懲罰的變量選擇方法：Bridge、Lasso、SCAD、Elastic Net、Adaptive Lasso、Dantzig Selector、MCP等。雖然Lasso具有較好的預(yù)測結(jié)果，但一般情況下Lasso結(jié)果是有偏的，在嚴(yán)格的條件假設(shè)下才具有相合性，并且Lasso不具有oracle性質(zhì)[1]。而SCAD、Adaptive Lasso、Dantzig Selector以及MCP具有oracle性質(zhì)。West等[2]提出當(dāng)“l(fā)arge p, small n”時，要特別重視成組變量(Grouped variables)的問題。針對此類問題，Hastie等[3]、Díaz-Uriarte[4]提出了主成分分析法，Hastie等[5]提出了監(jiān)督tree harvesting方法，Dettling和Bühlmann[6]將聚類和有監(jiān)督學(xué)習(xí)結(jié)合在了一起，Segal等[7]論證了正則化回歸方法處理成組變量的優(yōu)勢，Zou和Hastie[8]提出了變量選擇方法的組效應(yīng)(Grouping effect)。Lasso處理成組變量的效果非常差，而SCAD、Elastic Net、Adaptive Lasso、Dantizig Selector、MCP中，僅Elastic Net方法具有組效應(yīng)。

由于生存數(shù)據(jù)的刪失性，完全數(shù)據(jù)的變量選擇方法不能直接應(yīng)用于生存數(shù)據(jù)，因而一些學(xué)者研究了Cox比例風(fēng)險模型(簡稱之為Cox模型)的變量選擇問題：Tibshirani[9]、Fan和Li[10]分別將Lasso、SCAD應(yīng)用于Cox模型，Li和Luan[11]提出了Cox核轉(zhuǎn)換方法，Zhang和Lu[12], Zou[13]將Adaptive Lasso應(yīng)用到Cox模型，閆麗娜[14]將Elastic Net與Cox模型結(jié)合，侯文[15]研究了Cox模型的橋估計(Bridge)，鄧秋玲[16]研究了SCAD和Adaptive Dantizig Selector在Cox模型中的應(yīng)用，劉丹等[17]研究了Adaptive Elastic Net在Cox模型中的應(yīng)用。

為了使高維生存數(shù)據(jù)的Cox模型的變量選擇方法既有oracle性質(zhì)，又具有組效應(yīng)，本文提出了Cox模型的彈性SCAD方法，并證明了彈性SCAD方法的統(tǒng)計性質(zhì)，通過數(shù)值模擬，比較了在Cox模型下，彈性SCAD與Lasso、Adaptive Lasso、SCAD、Elastic Net、Adaptive Elastic Net等方法的變量選擇結(jié)果，得到了彈性SCAD在某些情況下的優(yōu)越性，最后再結(jié)合實(shí)例，探討了彈性SCAD及其他變量選擇方法應(yīng)用于Cox模型處理生存數(shù)據(jù)的表現(xiàn)優(yōu)劣。

1 基于Cox模型的彈性SCAD方法及性質(zhì)

首先，我們給出Cox模型的彈性SCAD方法的估計：

(1)

R(ti)表示ti時刻的危險集，δi是表示是否刪失的示性函數(shù)。γ>2,λ1>0,λ2>0為調(diào)整參數(shù)

xi=(xi1，xi2，…，xip)，表示矩陣X的第i行。

故

(2)

(3)

記Q(β)=-ln(β)+Pλ1,λ2(β)。

(4)

我們可以得出基于Cox比例風(fēng)險模型的彈性SCAD方法有以下性質(zhì)(證明見附錄)：

a=max{|xij-xiτ|},i∈{1,2,…,n},j,τ∈{1,2,…,p}。

性質(zhì)3(估計一致性) 若n→+,則存在Q(β)的一個局部最優(yōu)解滿足其中

性質(zhì)4(Oracle性質(zhì)) 若n→,,則滿足

性質(zhì)5(方差估計的漸近性)

性質(zhì)1～5的證明見附錄。

2 數(shù)值模擬

彈性SCAD方法可以看作是SCAD方法和嶺回歸的結(jié)合，具有Oracle性質(zhì)，其克服了Elastic Net有偏估計的缺點(diǎn)，屬于漸近無偏估計。與SCAD相比，具有組效應(yīng)，且與同樣具有組效應(yīng)的Elastic Net(ENet)、Adaptive elastic net(AENet)相比，其在處理小樣本高維數(shù)據(jù)、變量間存在高度共線性問題與小樣本低維、變量間存在弱相關(guān)關(guān)系兩種情況下表現(xiàn)更優(yōu)。

本部分?jǐn)?shù)據(jù)模擬選用十折交叉驗(yàn)證，研究的目的是通過運(yùn)用基于Cox模型的Lasso、Adaptive lasso(Alasso)、ENet、AENet、SCAD、彈性SCAD(Escad)6種變量選擇方法，對模擬生成的高維、并具有不同共線性強(qiáng)弱的生存數(shù)據(jù)進(jìn)行變量選擇，比較其變量篩選效果以及模型誤差等指標(biāo)，進(jìn)而評價各種方法的優(yōu)劣。

數(shù)值模擬的參數(shù)設(shè)置情況如下：設(shè)樣本量為n,協(xié)變量個數(shù)為p，生成n×p的數(shù)據(jù)矩陣X，除前2列外，其余各列服從標(biāo)準(zhǔn)的多元正態(tài)分布，第一列數(shù)據(jù)服從一元標(biāo)準(zhǔn)正態(tài)分布，第二列數(shù)據(jù)與第一列數(shù)據(jù)之間的相關(guān)系數(shù)分別設(shè)為r(v1,v2)=0.8、0.5和0.2，分別代表變量和間存在強(qiáng)共線性、中等共線性和弱共線性的情況。參數(shù)設(shè)置為:

v=(v1,v2,v3,…,vp)=(0.8,0.7,1,-0.6,0,0,…,0)。

當(dāng)變量v1與v2間存在共線性且為重要變量時，若所用變量選擇方法同時將v1、v2選進(jìn)模型，則說明該方法具有變量選擇的組效應(yīng)，否則，說明其沒有變量選擇的組效應(yīng)。生存時間刪失率設(shè)為0.3。

表1 部分?jǐn)?shù)值模擬結(jié)果

從圖1可以看出，SCAD相比于Lasso、Alasso，其模型誤差較小，模型穩(wěn)定性強(qiáng)；彈性SCAD相比于Enet、AENet，由于其具有Oracle性質(zhì)，故模型誤差較小，模型穩(wěn)定性強(qiáng)，雖然AENet也具有Oracle性質(zhì)，但模型表現(xiàn)不如Escad穩(wěn)定。

圖1 6種方法在模型誤差方面的表現(xiàn)(相關(guān)系數(shù)為0.8)

圖2給出了當(dāng)n=100,p=10且相關(guān)系數(shù)為0.8時，彈性SCAD與SCAD的求解路徑的比較，其中圖2(a)為彈性SCAD，圖2(b)為SCAD。

(n=100,p=10,r(v1,v2)=0.8。曲線1～4分別為v1～v4的估計值。n=100, p=10, r(v1,v2)=0.8.line 1～4 are the estimators of variables v1～v4.)

結(jié)合圖1和2可以看出，彈性SCAD在保留了SCAD變量選擇優(yōu)點(diǎn)的同時，克服了SCAD方法在進(jìn)行變量選擇時不具有組效應(yīng)的缺點(diǎn)，在變量間存在高度相關(guān)性時，能夠把相關(guān)的變量同時選進(jìn)模型，具有變量選擇的組效應(yīng)。

從數(shù)據(jù)模擬結(jié)果可以看出：不具有變量選擇組效應(yīng)的方法是基于Cox模型的Lasso、Alasso、SCAD；具有變量選擇組效應(yīng)的方法是ENet、AENet、彈性SCAD。其中三者在對于有組效應(yīng)的變量選擇方法，在n=100,p=50,變量間存在強(qiáng)相關(guān)關(guān)系與n=100,p=10,變量間存在弱相關(guān)關(guān)系兩種情況下，彈性SCAD模型誤差最小，系數(shù)估計方面也表現(xiàn)最佳；在n=500,p=10情況、變量間存在非強(qiáng)(中等強(qiáng)度及較弱強(qiáng)度)相關(guān)關(guān)系時，彈性SCAD在系數(shù)估計方面表現(xiàn)最佳；在n=100,p=50情況，變量間存在弱相關(guān)關(guān)系時，彈性SCAD的模型誤差最小?；贑ox模型的Lasso和Alasso均傾向于多選變量，而SCAD與二者相比，雖然也存在假陽性，但除小樣本低維、變量間存在弱相關(guān)性的情況外，其假陽性的個數(shù)均小于二者，故變量選擇的一致性方面，SCAD優(yōu)勢明顯；基于Cox模型的ENet、AENet和彈性SCAD均傾向于多選變量，而彈性SCAD與二者相比，雖然也存在假陽性，但在n=100,p=50情況下假陽性的個數(shù)均小于二者，故變量選擇的一致性方面，此情況下基于Cox模型的彈性SCAD方法最優(yōu)。當(dāng)n=100,p=200時，不具有變量選擇組效應(yīng)的方法是基于Cox模型的Lasso、Alasso、SCAD；具有變量選擇組效應(yīng)的方法仍是基于Cox模型的ENet、AENet、彈性SCAD；在變量間3種不同程度的相關(guān)關(guān)系情況下，基于Cox模型的彈性SCAD相比于ENet 和AENet方法，均具有較少的噪聲系數(shù)、較低的模型誤差。

3 實(shí)例分析

本實(shí)例數(shù)據(jù)來源于Kalbfleish 和Prentice的一組肺癌治療方案的臨床試驗(yàn)數(shù)據(jù)[18]，這是一個標(biāo)準(zhǔn)的生存分析數(shù)據(jù)集。我們對這組數(shù)據(jù)利用前文提到的6種方法進(jìn)行變量選擇，結(jié)果見表2。

由于Diagnosis time, age和prior的值均為0，可以判定其對研究對象的生存時間沒有影響；結(jié)合Escad、Enet、AENet的選擇結(jié)果，可知small為噪聲系數(shù)，對結(jié)果沒有明顯影響；若假設(shè)變量間存在相關(guān)性，結(jié)合數(shù)值模擬結(jié)果可知，在小樣本低維數(shù)情況下，AENet表現(xiàn)最優(yōu)，故理論上其選擇的非零個數(shù)要比SCAD、Lasso、Alasso要多，而結(jié)合表2中數(shù)據(jù)可知，這4種方法選出的非零個數(shù)并無明顯差異，甚至不具有組效應(yīng)的Lasso要比AENet選的還多，這說明這些變量間是不具有相關(guān)關(guān)系。

表2 對肺癌數(shù)據(jù)進(jìn)行變量選擇的結(jié)果

卡式評分(Karnofsky score)表現(xiàn)得分越高，健康狀況越好，越能忍受治療給身體帶來的副作用，因而也就有可能接受徹底的治療，6種方法均將其選出，比較符合實(shí)際；肺小細(xì)胞癌是肺癌中最兇惡的，壞死典型且呈廣泛性，擴(kuò)散轉(zhuǎn)移快；而與之對應(yīng)的，鱗狀細(xì)胞癌(鱗癌)、腺癌、大細(xì)胞癌均為非小細(xì)胞型肺癌，與小細(xì)胞癌相比，它們擴(kuò)散轉(zhuǎn)移相對較晚。故根據(jù)實(shí)際，可把鱗狀細(xì)胞癌、胰癌、大細(xì)胞癌看作一組變量，表2中只有Escad、ENet、Lasso 3種方法同時選進(jìn)了這3種變量，比較符合實(shí)際。但Lasso不具有組效應(yīng)，故其雖然選進(jìn)了這3種變量，應(yīng)該是作為噪聲變量選入的，系數(shù)估計會很差，故與實(shí)際相差比較大；而AENet作為在ENet的基礎(chǔ)上的一種改進(jìn)，卻沒有將這3種因素同時選入模型，可見ENet選入這3種變量是因?yàn)榇朔椒〞噙x變量。綜上，可推知影響研究對象生存時間的主要因素有：Treatment、Karnofsky、squamous、adeno和large影響效用系數(shù)分別為：0.099、-0.028、-0.683、0.296和-0.334。

即癌癥類型是否為非小細(xì)胞癌，是決定生存時間的因素。更進(jìn)一步，其中是否為鱗狀細(xì)胞癌的指標(biāo)為主要決定因素。Karnofsky表現(xiàn)得分雖然對生存時間有影響，但是相對來說其影響是次要的，不起決定性作用。再對比文獻(xiàn)[12]中研究結(jié)果，可知本文上述研究結(jié)果與之相符，但治療指標(biāo)trt是否為有效影響因素，對比之下難以確定，有待進(jìn)一步研究。由于在實(shí)際應(yīng)用中，變量間的關(guān)系十分復(fù)雜，故不會完全與模擬數(shù)據(jù)中的表現(xiàn)完全一致，經(jīng)本例的分析結(jié)果，可推知基于Cox模型的彈性SCAD在實(shí)際處理生存數(shù)據(jù)時表現(xiàn)優(yōu)于文中其他變量選擇方法。

4 結(jié)論與展望

本文主要研究了基于Cox模型的彈性SCAD變量選擇方法的理論性質(zhì)，通過數(shù)值模擬比較得出了如下結(jié)論：在n=100，p=50情況下，基于Cox模型的Elastic net和Adaptive elastic net方法在變量選擇的一致性上的表現(xiàn)不如彈性SCAD方法；對于有組效應(yīng)的變量選擇方法，在n=100，p=50情況，變量間存在強(qiáng)相關(guān)關(guān)系與小樣本低維、變量間存在弱相關(guān)關(guān)系兩種情況下，基于Cox模型的彈性SCAD模型誤差最小，系數(shù)估計方面也表現(xiàn)最佳；在n=500，p=10情況，變量間存在非強(qiáng)(中等強(qiáng)度及較弱強(qiáng)度)相關(guān)關(guān)系時，基于Cox模型的彈性SCAD在系數(shù)估計方面表現(xiàn)最佳；在n=100，p=50情況、變量間存在弱相關(guān)關(guān)系時，基于Cox模型的彈性SCAD的模型誤差最小。

當(dāng)n=100，p=200時，在變量間3種不同程度的相關(guān)情況下，基于Cox模型的彈性SCAD相比于Elastic net 和Adaptive elastic net方法，均具有較少的噪聲系數(shù)、較低的模型誤差。進(jìn)一步通過實(shí)例分析發(fā)現(xiàn)，基于Cox模型的彈性SCAD的變量選擇結(jié)果優(yōu)于文中討論的其余變量選擇方法，變量選擇結(jié)果較為合理，與實(shí)際更相符。

通過多種方法的實(shí)例比較分析，可以判斷變量間是否存在共線性，有利于在決策時選擇更適合的方法，作出更理性、正確的判斷，進(jìn)一步可知，通過比較這幾種方法在不同類型參數(shù)下的表現(xiàn)，可判斷共線性的強(qiáng)弱，在實(shí)際應(yīng)用時，提供一個較為理性的方案。

本文只選取了特定刪失比例、變量個數(shù)、數(shù)據(jù)個數(shù)進(jìn)行了數(shù)值模擬，彈性SCAD與SCAD的模擬程序也不適于處理帶有節(jié)點(diǎn)的生存數(shù)據(jù)，且由于設(shè)備的局限性，也未能研究彈性SCAD在更高維數(shù)據(jù)下的表現(xiàn)，結(jié)果難免會有一定的片面性。在進(jìn)一步的研究中，可改進(jìn)算法、程序、改變刪失比例、變量個數(shù)及數(shù)據(jù)個數(shù)進(jìn)行更全面、高效的數(shù)值模擬，發(fā)現(xiàn)其新的變量選擇特性或局限性所在，從而更好地應(yīng)用于理論研究和實(shí)踐方面。