鄒麗英，劉 祎

(中國海洋大學(xué)數(shù)學(xué)科學(xué)學(xué)院，山東 青島 266100)

隨著科技水平的飛速發(fā)展和數(shù)據(jù)收集能力的大幅提升，超高維數(shù)據(jù)已經(jīng)越來越頻繁地出現(xiàn)在包括金融學(xué)、基因?qū)W、醫(yī)學(xué)等各領(lǐng)域中。在高維數(shù)據(jù)中，為了選出對響應(yīng)變量有重要影響的少數(shù)預(yù)測變量，許多懲罰方法被提出，例如LASSO[1]、SCAD[2]、自適應(yīng)LASSO[3]和Dantzing[4]等。但在超高維數(shù)據(jù)中，即當(dāng)變量個數(shù)相對于樣本量呈指數(shù)增長時，這些方法面臨著計算復(fù)雜性、統(tǒng)計準(zhǔn)確性和算法穩(wěn)定性的共同挑戰(zhàn)[5]。為此許多學(xué)者提出超高維數(shù)據(jù)的特征篩選方法，用于解決面臨的困難。Fan和Lv[6]基于線性模型提出了一種依賴于邊際皮爾遜相關(guān)系數(shù)的確定獨(dú)立篩選(Sure independent screening, SIS)，它可以將超高維數(shù)據(jù)降維至合適的大小，并且通過對相關(guān)系數(shù)絕對值大小進(jìn)行排序，將對響應(yīng)變量有重要影響的預(yù)測因子選中的概率趨于1。Zhu等[7]用一種針對一般多指標(biāo)模型確定獨(dú)立排序和篩選技術(shù)(Sure independent ranking and screening,SIRS)來排序重要變量。Li等[8]基于距離相關(guān)系數(shù)發(fā)展了SIS程序，但它對重尾數(shù)據(jù)不具有穩(wěn)健性。Zhong等[9]在穩(wěn)健距離相關(guān)的基礎(chǔ)上提出了一種應(yīng)用于單指標(biāo)模型的特征篩選方法。然而上述特征篩選方法都是針對完全可觀測數(shù)據(jù)的。

在實(shí)際案例中，例如成本效益分析、教育調(diào)查、全基因組關(guān)聯(lián)研究以及基因表達(dá)研究等領(lǐng)域，由于部分受試者不愿意回答敏感問題，或者不可控因素導(dǎo)致的信息丟失，出現(xiàn)響應(yīng)變量隨機(jī)缺失(Missing at random，MAR)的情況是較為常見的。有關(guān)于缺失數(shù)據(jù)的統(tǒng)計分析，已有許多文獻(xiàn)進(jìn)行研究，例如：Wang和Rao[10]提出了處理缺失響應(yīng)問題的經(jīng)驗(yàn)似然方法。Qin、Shao和Zhang[11]將逆概率加權(quán)方法用于處理協(xié)變量相關(guān)的缺失響應(yīng)數(shù)據(jù)。Hu、Follmann和Qin[12]研究了缺失數(shù)據(jù)平均響應(yīng)的半?yún)?shù)降維估計。近年來，一系列研究開始集中于處理缺失響應(yīng)變量的超高維數(shù)據(jù)，并提出高效、準(zhǔn)確的特征篩選方法。Lai等[13]調(diào)整了Zhu等[7]的SIRS方法并結(jié)合逆概率加權(quán)技術(shù)，提出了一種無模型特征篩選方法。逆概率加權(quán)方法對缺失概率較敏感，因此非參數(shù)插補(bǔ)作為一種處理缺失數(shù)據(jù)的方法得到廣泛應(yīng)用。Fang[14]基于非參數(shù)插補(bǔ)技術(shù)提出了無模型的特征篩選方法，并說明其比逆概率加權(quán)方法具有更好的篩選效果。

本文提出了一種新的特征篩選程序(Imputed distance correlation，IDC)，構(gòu)造插補(bǔ)響應(yīng)變量與協(xié)變量分布函數(shù)之間的距離相關(guān)系數(shù)作為篩選指標(biāo)進(jìn)行特征篩選。所提方法不依賴于模型假設(shè)且對協(xié)變量異常值穩(wěn)健，而且可以直接處理響應(yīng)變量是多維的情形。通過數(shù)值模擬和微陣列彌漫性大B細(xì)胞淋巴瘤(Diffuse Large-B-Cell Lymphoma，DLBCL)數(shù)據(jù)分析，對方法的有限樣本性質(zhì)進(jìn)行了驗(yàn)證。

1 方法

記Y為連續(xù)的響應(yīng)變量，X=(X1,…,Xp)T為p×1維的連續(xù)預(yù)測變量。假設(shè)維數(shù)p相對于樣本容量n呈指數(shù)增長，即log(p)=Ο(nα)，常數(shù)α>0。其中X是完全可觀測的，而Y可能缺失。由于觀測數(shù)據(jù)是不完全的，記(Xi,Yi,δi)為樣本數(shù)據(jù)，i=1,…,n，其中δi為缺失響應(yīng)的指示變量，即如果Yi缺失，則δi=0，如果Yi可觀測，則δi=1。我們假設(shè)δ只依賴于X使得傾向得分函數(shù)有π(X)=P(δ=1|X)的形式。由Little和Rubin[15]可知，上述關(guān)于缺失機(jī)制的定義為隨機(jī)缺失。Y是隨機(jī)缺失的，簡單來說，就是假設(shè)

P(δ=1|X,Y)=P(δ=1|X)。

基于稀疏假設(shè)，只有少數(shù)的預(yù)測變量與Y有關(guān)。我們定義活躍預(yù)測變量的索引集：

A={k:F(y|X)依賴于Xk,對某個y∈ψY,k=1,…,p}，

其中F(y|X)=P(Y≤y|X)為給定變量X條件下變量Y的條件分布函數(shù)，ψY為Y的支撐集，且|A|為集合A的勢。在超高維數(shù)據(jù)分析中我們假設(shè)p?n且p?|A|。記Ac={1,2,…,p}A為非活躍預(yù)測變量的索引集，即Ac={k:F(y|X)不依賴于Xk,對?y∈ψY,k=1,…,p}。令XA={Xk:k∈A}且XAc={Xk:k∈Ac}分別為預(yù)測變量的活躍集和非活躍集。

首先回顧一下距離相關(guān)系數(shù)[16]的定義。假設(shè)U和V為兩個隨機(jī)向量，維度分別為dU和dV。令φU(u)和φV(v)分別為U和V的特征函數(shù)，φU,V(u,v)為U和V的聯(lián)合特征函數(shù)。距離協(xié)方差(Distance covaria-nce)定義為非負(fù)數(shù)dcov(U,V)，即

dcov2(U,V)=

U和V的距離相關(guān)系數(shù)(Distance correlation，DC)定義為

DC作為相關(guān)關(guān)系的一種度量，具有良好的特性，即dcorr(U,V)=0當(dāng)且僅當(dāng)U和V是相互獨(dú)立的。這一性質(zhì)使得距離相關(guān)尤其適用于超高維數(shù)據(jù)的變量篩選。并且對于超高維的完整數(shù)據(jù)集，Zhong等[9]在穩(wěn)健距離相關(guān)的基礎(chǔ)上提出了一種應(yīng)用于單指標(biāo)模型的特征篩選程序,這啟發(fā)我們將穩(wěn)健的距離相關(guān)應(yīng)用于具有缺失響應(yīng)的超高維數(shù)據(jù)中。

本文采用的穩(wěn)健距離相關(guān)篩選指標(biāo)是指

ωk=dcorr2(Fk(Xk),Y)=

處理缺失數(shù)據(jù)的插補(bǔ)方法是通過非參數(shù)回歸得到的，設(shè)m(x)=E(Y|X=x)為給定X時Y的回歸函數(shù)，則由核回歸[17-18]可得到m(x)的估計為

根據(jù)插補(bǔ)后的數(shù)據(jù)集，我們得到篩選指標(biāo)的估計為

2 數(shù)值模擬

在本節(jié)中我們將本文提出的插補(bǔ)距離相關(guān)方法(Imputed distance correlation，IDC)與其他處理超高維缺失響應(yīng)方法作比較。對比方法有：Lai等[13]提出的基于逆概率加權(quán)的超高維特征篩選方法(Inverse probability weighted，IPW)，F(xiàn)ang[14]提出的基于插補(bǔ)技術(shù)的非參數(shù)特征篩選方法(Method of Imputation Technique，ITM)以及本文方法在完整數(shù)據(jù)下的情形(Distance correlation of full sample data，F(xiàn)DC)。在整個模擬研究中，我們將傾向得分函數(shù)P(δ=1|X)設(shè)置為邏輯回歸模型，即P(δ=1|X)=exp(θX)/(1+exp(θX))，通過改變θ來改變?nèi)笔?Missing rate，MR)，設(shè)置缺失率大約為0.2和0.4。對于每種情形，設(shè)置200次重復(fù)。

為了評估所提出方法的性能，我們采用Li等[8]的三個評估準(zhǔn)則。第一個評估準(zhǔn)則是包括所有活躍的預(yù)測變量的最小模型數(shù)，用S表示。我們給出了200次重復(fù)模擬中S的5%、25%、50%、75%和95%分位數(shù)。第二個評估準(zhǔn)則是每個活躍預(yù)測變量的覆蓋率，用Pi表示，我們給出了當(dāng)給定模型大小為dn時200次重復(fù)模擬中Xi被選擇的比例。第三個評估準(zhǔn)則是所有活躍預(yù)測變量的覆蓋率，用Pa表示。我們給出了在200次重復(fù)模擬中，對于給定模型大小dn，所有活躍預(yù)測變量均被選擇的比例。我們選擇dn為dn=[n/logn]，其中[x]表示x的整數(shù)部分。

實(shí)例1(線性模型)考慮線性模型如下：

Y=XTβ+ε。

表1 線性模型中單個活躍因子被選中的概率Ps及所有活躍因子被選中的概率Pa

從表1和2可以看出，在線性模型假設(shè)，以及協(xié)變量與響應(yīng)變量都存在重尾數(shù)據(jù)下，本文提出的IDC方法相比于IPW和ITM兩種方法而言，具有較為明顯的優(yōu)勢。具體表現(xiàn)在活躍預(yù)測變量入選模型的比例更高，所需要的最小模型數(shù)更小。說明IDC方法可以將重要變量排在不重要變量的前面。完整數(shù)據(jù)下應(yīng)用本文的FDC方法表現(xiàn)最好，因?yàn)槠淅昧藰颖镜乃行畔ⅰ?/p>

實(shí)例2(非線性模型1)考慮非線性模型如下：

其中β=(β1,β2,β3,0,…,0)T=(2,2.5,3,0,…,0)T∈Rp，即只有前三個變量是活躍的，X=(X1,…,Xp)T產(chǎn)生于均值為0，協(xié)方差矩陣為Σ的多元正態(tài)分布，ε服從N(0,1)。另外為進(jìn)一步獲得重尾預(yù)測變量，我們類似于Lai等[13]的設(shè)置，替換X1和X3為服從自由度為3的t分布產(chǎn)生的隨機(jī)樣本，替換X2為服從t(3)+1分布產(chǎn)生的隨機(jī)樣本，核函數(shù)與實(shí)例1相同。在此情形下，我們設(shè)置(θ1,θ2,θp)=(1,0.3,0.1)和(3,-0.5,2)得到約20%和40%的缺失率，選取(n,p)=(100,1 000)和(200,2 000)，模擬結(jié)果如表3和4所示。

從表3和4可以看出，IDC方法對非線性關(guān)系的檢測能力要明顯優(yōu)于ITM和IPW方法。這是由于不同于ITM和IPW篩選指標(biāo)的線性表達(dá)形式，IDC基于穩(wěn)健距離相關(guān)指標(biāo)，更適用于非線性模型的特征篩選。而且隨著樣本量的增大，方法的表現(xiàn)有明顯改善。

表2 線性模型中最小模型數(shù)S的各分位數(shù)

表3 非線性模型1中單個活躍因子被選中的概率Ps及所有活躍因子被選中的概率Pa

表4 非線性模型1中最小模型數(shù)S的各分位數(shù)

實(shí)例3(非線性模型2)為了進(jìn)一步驗(yàn)證方法在非線性模型上的性能，考慮非線性模型如下：

cos(β21X1+β22X2+β23tanhX3)+ε，

其中(β11,β12,β13,β21,β22,β23)=(0.8,1,1.2,1.6,1.4,1.2)，即只有前三個變量是活躍的，其余設(shè)置與實(shí)例2中的設(shè)置相同。在此情形下，我們選取樣本量n和協(xié)變量個數(shù)p為(n,p)=(100,1 000)和(200,1 000)，模擬結(jié)果見表5和6。

表5 非線性模型2中單個活躍因子被選中的概率Ps及所有活躍因子被選中的概率Pa

表6 非線性模型2中最小模型數(shù)S的各分位數(shù)

從表5和6可以看出，在更復(fù)雜的非線性模型中，IDC對比ITM和IPW兩種方法的優(yōu)勢更為明顯。ITM和IPW兩種方法的接近于，即基本不可能選出所有的活躍預(yù)測變量，這可能是由于模型中存在示性函數(shù)及雙曲函數(shù)，使得非線性關(guān)系復(fù)雜，變量篩選也變得更加困難。而且可以看到，隨著樣本量的增大，ITM和IPW兩種方法并沒有明顯的改善。

實(shí)例4(雙響應(yīng)模型)考慮雙響應(yīng)變量模型如下：

其中β1=(β11,β12,β13,0,…,0)=(1,0.8,0.6,0,…,0)T∈Rp，β2=(β21,β22,β23,0,…,0)=(1,1.5,2,0,…,0)T∈Rp，即只有前三個變量是活躍的，其余設(shè)置與實(shí)例2中的設(shè)置相同。另外，我們對ITM方法中的篩選指標(biāo)采取類似于文獻(xiàn)[13]中公式(7)的二范數(shù)處理方式，使其適用于雙響應(yīng)模型的篩選指標(biāo)計算。在此情形下，我們選取樣本量n和協(xié)變量個數(shù)p為(n,p)=(100,1 000)和(200,1 000)，模擬結(jié)果見表7和8。

表7 雙響應(yīng)模型中單個活躍因子被選中的概率Ps及所有活躍因子被選中的概率Pa

表8 雙響應(yīng)模型中最小模型數(shù)S的各分位數(shù)

從表7和8可以看出，在多響應(yīng)變量的非線性模型中，對比ITM和IPW兩種方法，當(dāng)缺失率提高時，插補(bǔ)技術(shù)較于逆概率加權(quán)方法表現(xiàn)較好，且在一定程度上也說明了IDC方法的指標(biāo)對于復(fù)雜模型更具穩(wěn)健性。

除了上述模擬，我們還驗(yàn)證了方法在小樣本下的表現(xiàn)。在實(shí)例1的模型設(shè)置下，我們考慮了(n,p)=(50,500)和(n,p)=(50,1 000)兩種情況，模擬結(jié)果見表9和10。從中可以看出各種變量篩選方法的表現(xiàn)趨勢與大樣本結(jié)果大致相同。

表9 在小樣本的線性模型中單個活躍因子被選中的概率Ps及所有活躍因子被選中的概率Pa

另外，在實(shí)例1線性模型(n,p)=(100,1 000)的情況下，我們增加模擬了基于多重插補(bǔ)[19]的IDC方法,其中多重插補(bǔ)方法選取為10重插補(bǔ)(MI10-DC)，從結(jié)果可以看出，IDC和MI10-DC方法都優(yōu)于IPW和ITM方法，這說明穩(wěn)健距離相關(guān)的指標(biāo)優(yōu)勢。但MI10-DC方法在Pa上表現(xiàn)劣于IDC方法，結(jié)合結(jié)果來看，可能是由于X3重尾分布對多重插補(bǔ)的抽樣過程產(chǎn)生了影響。模擬結(jié)果見表11。

表10 在小樣本的線性模型中最小模型數(shù)S的各分位數(shù)

表11 在線性模型中單個活躍因子被選中的概率Ps，所有活躍因子被選中的概率Pa及最小模型數(shù)S的各分位數(shù)

3 實(shí)際案例

我們應(yīng)用提出的方法來分析微陣列彌漫性大B細(xì)胞淋巴瘤(DLBCL)數(shù)據(jù)，見Rosenwald等[17]及Zhu等[7]。為了對比，我們也采用IPW方法、ITM方法和FDC方法進(jìn)行分析。DLBCL數(shù)據(jù)集共包含240例患者，其中隨訪期間死亡138例。因此，相應(yīng)的響應(yīng)變量由觀察到的生存時間和刪失指標(biāo)組成，從每個患者的cDNA微陣列中獲得的p=7 399個基因值是預(yù)測因子。另外，我們將數(shù)據(jù)劃分為包含n1=160名患者的訓(xùn)練集和包含n2=80名患者的測試集。由于預(yù)測因子數(shù)量p遠(yuǎn)大于樣本量n，為擬合原始模型，進(jìn)行特征篩選是有必要的。

為進(jìn)一步研究各種方法的表現(xiàn)，我們首先在訓(xùn)練集中分別利用IDC、ITM、IPW和FDC方法篩選基因并擬合Cox比例風(fēng)險模型。然后，在測試集評估各種方法的預(yù)測表現(xiàn)，計算測試集的風(fēng)險得分，將其分為低風(fēng)險組和高風(fēng)險組，其中，區(qū)分值是由訓(xùn)練集中估計得分的中值得到的。圖1展示了四種方法下兩個風(fēng)險組病人的生存曲線的Kaplan-Meier估計[21]。

從圖1我們可以看出，在MAR情況下，運(yùn)用IDC方法的曲線分離得最好，且其對數(shù)秩檢驗(yàn)產(chǎn)生的p值為0.001 2，而IPW方法和ITM方法對應(yīng)的p值分別為0.037 9和0.261 2。在完整數(shù)據(jù)集的情況下，F(xiàn)DC方法的p值為0.039 8。

(虛線代表高風(fēng)險組，實(shí)線代表低風(fēng)險組。The dotted line is high-risk group, and the solid line is low-risk group.)

4 結(jié)語

本文提出了一種新的特征篩選程序(IDC)，首先通過插補(bǔ)技術(shù)，補(bǔ)全缺失響應(yīng)變量值，再構(gòu)造插補(bǔ)響應(yīng)變量與協(xié)變量分布函數(shù)之間的距離相關(guān)系數(shù)，作為篩選指標(biāo)。與傳統(tǒng)的參數(shù)方法相比，所提出的非參數(shù)方法在穩(wěn)健性、可擴(kuò)展性和非線性重要變量檢測能力方面具有一定的優(yōu)勢。本文提出的IDC方法基于插補(bǔ)技術(shù)創(chuàng)建了一個完整的數(shù)據(jù)集，便于使用標(biāo)準(zhǔn)的完整數(shù)據(jù)算法。相比于逆概率加權(quán)方法，由于一次插補(bǔ)只為每個具有缺失值的變量創(chuàng)建一個偽觀察，因此其影響遠(yuǎn)小于基于缺失概率的逆概率加權(quán)方法的特征篩選程序，即我們的方法僅依賴于插補(bǔ)，減少了極端缺失概率的不利影響。通過數(shù)值模擬和實(shí)例可以看出，IDC方法較于對比方法在應(yīng)用上具有一定的優(yōu)勢，但在確定篩選理論性質(zhì)方面，基于穩(wěn)健距離相關(guān)發(fā)展缺失響應(yīng)變量的插補(bǔ)方法具有一定的難度，這也啟發(fā)我們進(jìn)一步改進(jìn)篩選指標(biāo)，優(yōu)化處理缺失響應(yīng)數(shù)據(jù)的超高維篩選方法。