張家睿，吳耀華

(1 中國(guó)科學(xué)技術(shù)大學(xué)管理學(xué)院， 合肥 230026; 2 香港大學(xué)浙江科學(xué)技術(shù)研究院， 杭州 310000)

在過(guò)去10年里分子生物學(xué)試驗(yàn)技術(shù)的進(jìn)展給我們帶來(lái)了豐富的生物醫(yī)學(xué)數(shù)據(jù)，舉例來(lái)說(shuō)，DNA顯微序列可以用來(lái)測(cè)量一個(gè)細(xì)胞中成千上萬(wàn)的基因。這種類型的數(shù)據(jù)中樣本維度p比樣本量n要大得多，對(duì)于傳統(tǒng)的統(tǒng)計(jì)推斷方法來(lái)說(shuō)是一個(gè)巨大的挑戰(zhàn)，有很多經(jīng)典的推斷方法在這種情況下變得不適用。這種情形下有效的變量選擇方法就變得尤為重要。比較著名的高維數(shù)據(jù)變量選擇方法有Lasso[1]，SCAD[2]和MCP[3]等。

當(dāng)研究關(guān)于患者生存狀態(tài)的醫(yī)療數(shù)據(jù)時(shí)，將高維的生物醫(yī)療數(shù)據(jù)和患者的生存狀態(tài)數(shù)據(jù)結(jié)合起來(lái)分析是一個(gè)很有效的方法。因此近些年來(lái)也有很多關(guān)于高維生存分析模型的變量選擇方法，比如Bradic等[4]關(guān)于高維Cox模型的正則化方法，Gorst-Rasmussen和Scheike[5]關(guān)于高維單指數(shù)模型的篩選方法，Lin和Lyu[6]關(guān)于高維可加模型的正則化方法等等。高維生存分析模型還廣泛地應(yīng)用到信用風(fēng)險(xiǎn)分析，比如Fan等[7]。

由于在實(shí)際生活中，我們經(jīng)常會(huì)遇到帶有測(cè)量誤差的數(shù)據(jù)，所以對(duì)于帶有測(cè)量誤差數(shù)據(jù)的分析方法也是一個(gè)重要的研究方向，對(duì)于高維線性模型有Loh和Wainwright[8]以及Datta和Zou[9]的相關(guān)工作；對(duì)于變系數(shù)模型，有劉智凡等[10]的工作。對(duì)于帶有測(cè)量誤差的生存分析數(shù)據(jù)的變量選擇方法，代表文章有Song和Wang[11]關(guān)于工具變量的工作，Chen和Yi[12]關(guān)于Cox模型左截?cái)嘤覄h失數(shù)據(jù)的工作。高維生存分析模型由于其計(jì)算復(fù)雜度較高以及理論性質(zhì)較為復(fù)雜，所以對(duì)于帶有測(cè)量誤差的高維生存分析數(shù)據(jù)的工作隨著近些年大數(shù)據(jù)的迅速發(fā)展才逐步出現(xiàn)在視野之中。具有代表性的文章有Chen和Yi[13]關(guān)于高維生存分析圖模型的工作以及Chen等[14]關(guān)于高維Cox模型利用糾正似然函數(shù)的工作。本文選擇同樣具有重要應(yīng)用的可加風(fēng)險(xiǎn)模型作為基礎(chǔ)，結(jié)合處理高維線性模型的正則化方法對(duì)帶有測(cè)量誤差的生存分析數(shù)據(jù)進(jìn)行分析。

1 研究背景

本文所采用的模型為高維可加風(fēng)險(xiǎn)模型，結(jié)合高維線性模型測(cè)量誤差處理辦法對(duì)帶有測(cè)量誤差的生存分析數(shù)據(jù)進(jìn)行分析。下面對(duì)高維可加風(fēng)險(xiǎn)模型和高維線性模型測(cè)量誤差處理方法分別進(jìn)行介紹。

1.1 高維可加風(fēng)險(xiǎn)模型

對(duì)于生存分析數(shù)據(jù)的變量選擇技術(shù)的發(fā)展已經(jīng)不拘泥于Cox模型，可加風(fēng)險(xiǎn)模型便是除Cox模型以外的一種重要替代方式?？杉语L(fēng)險(xiǎn)模型假設(shè)失效時(shí)間為T(mén)的風(fēng)險(xiǎn)函數(shù)和p維的協(xié)變量X(·)有如下形式的關(guān)系

(1)

其中：λ0(·)是一個(gè)不確定的基線風(fēng)險(xiǎn)函數(shù)，β0是一個(gè)p維的回歸系數(shù)。令C為刪失時(shí)間，則定義刪失失效時(shí)間為CFT=C∧T，令CFT=t1,…,tn，失效指數(shù)定義為δ=I(T≤C)，其中I(·)為指示函數(shù)，令X(t)=(X1(t),…,Xp(t))并且假設(shè)給定X觀察到的數(shù)據(jù)為(CFT,δ,X(·))，風(fēng)險(xiǎn)函數(shù)由式(1)給出。

采用常用的計(jì)數(shù)手段，定義觀察到的失效計(jì)數(shù)序列為Ni(t)=I(ti≤t,δi=1)，風(fēng)險(xiǎn)中指數(shù)為Yi(t)=I(ti≥t)，計(jì)數(shù)過(guò)程鞅為

(2)

后文也將用N(t),Y(t)和M(t)來(lái)代表這些計(jì)數(shù)過(guò)程的廣義形式。

Lin和Ying[15]采用一種有如下形式的偽得分方程來(lái)對(duì)可加風(fēng)險(xiǎn)模型進(jìn)行分析：

{dNi(t)-Yi(t)βTXi(t)dt},

(3)

其中β∈p，并且

(4)

τ是最大的跟蹤時(shí)間(生存時(shí)間和刪失時(shí)間的最大值)。這個(gè)估計(jì)函數(shù)關(guān)于回歸系數(shù)是線性的，令

(5)

和

(6)

其中v?2=vvT，通過(guò)一些代數(shù)變換，可以寫(xiě)出如下等式

U0(β)=b0-V0β.

(7)

在沒(méi)有測(cè)量誤差的情況下，V0是半正定的，式(7)兩邊關(guān)于β積分就可以得到損失函數(shù)

(8)

Leng和Ma[16]以及Martinussen和Scheike[17]都建議用上述損失函數(shù)配合正則化方法對(duì)可加風(fēng)險(xiǎn)模型(1)進(jìn)行變量選擇。本文的相關(guān)工作也是在此基礎(chǔ)上進(jìn)行。

1.2 高維線性模型測(cè)量誤差數(shù)據(jù)的處理方法

為了進(jìn)一步構(gòu)建更深層次的討論，假設(shè)觀察到的是被污染的協(xié)變量矩陣

Z(·)=(zij(·))1≤i≤n，1≤j≤p,

(9)

而不是真實(shí)的協(xié)變量矩陣X(·)。有很多種造成測(cè)量誤差的途徑，在加法測(cè)量誤差設(shè)定中，zi,j(·)=xi,j(·)+ai,j，其中A(·)=(ai,j)是加法測(cè)量誤差。在乘法測(cè)量誤差設(shè)定中，zi,j(·)=xi,j(·)mi,j，其中mi,j就是乘法測(cè)量誤差。缺失數(shù)據(jù)可以看作乘法測(cè)量誤差的一個(gè)特殊形式，mi,j=I(xi,j(·)沒(méi)缺失)。

不失一般性，用Lasso算法來(lái)舉例說(shuō)明測(cè)量誤差的影響，對(duì)于線性模型y=Xβ+來(lái)說(shuō)，Lasso算法是最小化

(10)

這等價(jià)于最小化

(11)

(12)

然后解決下面的優(yōu)化問(wèn)題來(lái)得到β的估計(jì)：

(13)

(14)

其中R是一個(gè)跟稀疏度有關(guān)的常數(shù)。Datta和Zou[9]提出一種最近鄰正定投影矩陣的算法來(lái)解決上述問(wèn)題，對(duì)于任意方陣K:

(15)

(16)

(17)

2 帶有測(cè)量誤差的高維可加風(fēng)險(xiǎn)模型的變量選擇方法

2.1 簡(jiǎn)化偽得分方程

在第1節(jié)中已經(jīng)介紹了Lin和Ying[15]的偽得分方程的具體形式，下面將在協(xié)變量X期望值為0的前提下簡(jiǎn)化該偽得分方程，提出一種全新的更加容易計(jì)算且符合實(shí)際情況的損失函數(shù)。首先定義

(18)

以及

(19)

則有

(20)

接著定義

(21)

由于X的期望為0，所以容易得到E(U(β))=0，在如上定義的基礎(chǔ)上，類似于式(7)，有

U(β)=b-Vβ,

(22)

式(22)對(duì)β積分即可得到期望為0時(shí)的損失函數(shù)

(23)

綜上所述即為簡(jiǎn)化版本的損失函數(shù)，我們將基于這個(gè)損失函數(shù)進(jìn)行變量選擇。

2.2 兩種測(cè)量誤差數(shù)據(jù)的變量選擇方法

2.2.1 加法測(cè)量誤差

假設(shè)觀測(cè)到的設(shè)計(jì)矩陣Z(·)被加法測(cè)量誤差污染，即zi,j(·)=xi,j(·)+ai,j，其中A(·)=(ai,j)。同時(shí)假設(shè)A的行是獨(dú)立同分布的，均值是0，協(xié)方差矩陣是ΣA，次高斯參數(shù)是τ2。假設(shè)ΣA是已知的，則V和b的無(wú)偏估計(jì)分別為

(24)

和

(25)

(26)

2.2.2 乘法測(cè)量誤差

(27)

以及

(28)

其中∥代表向量或者矩陣對(duì)應(yīng)元素相除。和加法測(cè)量誤差模型類似，乘法測(cè)量誤差下無(wú)偏估計(jì)矩陣也有可能不是正定的，所以基于Datta和Zou[9]的方法，可以得到相應(yīng)的凸損失函數(shù)：

(29)

3 理論性質(zhì)

在這一節(jié)中給出并推導(dǎo)估計(jì)量的l1和l2誤差界。記我們的估計(jì)量為CoCo估計(jì)量。首先定義近鄰條件：

(30)

(31)

對(duì)所有1≤i,j≤p成立。其中集合S={1,2,…，s}是回歸系數(shù)β的支撐集。

同樣也需要和線性模型下一樣的特征值限制條件：

條件3.2協(xié)方差陣特征值限制條件

(32)

條件3.2是一個(gè)在高維線性模型變量選擇中比較常見(jiàn)的假設(shè)。下面給出CoCo估計(jì)量的統(tǒng)計(jì)誤差界：

定理3.1在式(30)、式(31)和式(32)成立的前提下，對(duì)于λ≤min(ε0,12ε0‖βS‖∞)和ε≤min(ε0,Ω/64s)，下式至少以概率

(33)

其中

(34)

引理3.1說(shuō)明加法測(cè)量誤差的計(jì)算方法滿足近鄰條件。下面將對(duì)乘法測(cè)量誤差進(jìn)行說(shuō)明。為了保證乘法測(cè)量誤差的計(jì)算方法也滿足近鄰條件，需要添加額外的正則化條件如下：

(35)

則接下來(lái)有

引理3.2說(shuō)明了乘法測(cè)量誤差的計(jì)算方法滿足近鄰條件。將引理3.1，引理3.2和定理3.1結(jié)合有

推論3.1給出了加法測(cè)量誤差估計(jì)方法和乘法測(cè)量誤差估計(jì)方法的理論保證，確定了估計(jì)量l1和l2的誤差界，下面將通過(guò)隨機(jī)模擬實(shí)驗(yàn)和實(shí)際數(shù)據(jù)分析來(lái)驗(yàn)證我們的理論結(jié)果。

4 實(shí)驗(yàn)及結(jié)果分析

本文的方法簡(jiǎn)記為CoCo，Loh和Wainwright[8]的方法記為NCL，在隨機(jī)模擬實(shí)驗(yàn)和實(shí)際數(shù)據(jù)分析中將對(duì)兩種方法進(jìn)行比較。

4.1 隨機(jī)模擬

4.1.1 加法測(cè)量誤差模型

從可加風(fēng)險(xiǎn)模型中產(chǎn)生數(shù)據(jù)，設(shè)定λ0=5，回歸系數(shù)為

β=(3,1.5,0,0,2,…,0).

(36)

樣本量n=100，樣本維度p=200，X的行獨(dú)立同分布，均值為0，協(xié)方差矩陣為ΣX，考慮兩種情形下的ΣX：自回歸(ΣX,ij=0.5|i-j|)和復(fù)合對(duì)稱(ΣX,ij=0.5+I(i=j)*0.5)，刪失時(shí)間服從U(0,2)的均勻分布使得刪失率維持在20%左右。首先生成3n×p的X，然后從中選出n個(gè)滿足λ0+βTX>0的樣本作為實(shí)驗(yàn)數(shù)據(jù)。加法測(cè)量誤差為矩陣A，觀測(cè)數(shù)據(jù)由Z=X+A生成，A的行是服從N(0,τ2I)的獨(dú)立同分布變量，其中τ=0.25、0.5和0.75。

表1展示了CoCo和NCL兩種方法分別在自回歸和復(fù)合對(duì)稱條件下的100次重復(fù)實(shí)驗(yàn)的結(jié)果，可以看出在兩種情形下本文方法的選對(duì)數(shù)量和估計(jì)的均方誤差方面都比NCL方法要好。

表1 加法測(cè)量誤差兩種方法的結(jié)果Table 1 The results of two methods under additive error-in-variable data

4.1.2 乘法測(cè)量誤差模型

與加法測(cè)量誤差模擬類似，依舊從可加風(fēng)險(xiǎn)模型中產(chǎn)生數(shù)據(jù)，λ0=5，回歸系數(shù)，樣本量和樣本維度都保持不變，X的行獨(dú)立同分布，均值為0，協(xié)方差矩陣為ΣX，依舊考慮ΣX在自回歸和復(fù)合對(duì)稱兩種條件下的情形，并且與加法測(cè)量誤差中的設(shè)定保持一致。刪失時(shí)間服從U(0,2)的均勻分布使得刪失率維持在20%左右，首先生成3n×p的X，然后從中選出n個(gè)滿足λ0+βTX的作為實(shí)驗(yàn)數(shù)據(jù)。乘法測(cè)量誤差矩陣為M=((mi,j))，觀測(cè)數(shù)據(jù)由Z(·)=X(·)⊙M生成，log(mi,j)是服從N(0,τ2I)的獨(dú)立同分布變量，其中τ=0.25、0.5和0.75。與上一個(gè)隨機(jī)模擬實(shí)驗(yàn)一樣，依舊采用5折的交叉驗(yàn)證方法來(lái)估計(jì)CoCo估計(jì)量和NCL的參數(shù)R。同樣記錄C和IC分別代表選對(duì)的系數(shù)數(shù)量和錯(cuò)誤的數(shù)量，還記錄均方誤差(MSE)以及其標(biāo)準(zhǔn)差(se)。總共進(jìn)行100次實(shí)驗(yàn)取平均數(shù)作為最后的結(jié)果，在表2中展示。

表2展示了乘法測(cè)量誤差中，CoCo和NCL兩種方法分別在自回歸和復(fù)合對(duì)稱條件下的100次重復(fù)實(shí)驗(yàn)結(jié)果，可以看出在兩種情形下本文方法的選對(duì)數(shù)量和估計(jì)的均方誤差都比NCL方法要好。但是隨著測(cè)量誤差變大，CoCo和NCL方法的估計(jì)精確度都會(huì)有明顯下降。

表2 乘法測(cè)量誤差兩種方法的結(jié)果Table 2 The results of two methods under multiplicative error-in-variable data

4.2 實(shí)際數(shù)據(jù)分析

為了檢驗(yàn)我們方法的有效性，將295個(gè)樣本隨機(jī)分成包含235個(gè)樣本的訓(xùn)練集和60個(gè)樣本的驗(yàn)證集并重復(fù)100次，在每一次實(shí)驗(yàn)中，都采用隨機(jī)模擬實(shí)驗(yàn)中的兩種方法，即CoCo和NCL，用訓(xùn)練集訓(xùn)練模型參數(shù)并用驗(yàn)證集來(lái)篩選表現(xiàn)最好的估計(jì)量。計(jì)算

(37)

作為檢驗(yàn)兩種方法效果的指標(biāo)。具體的結(jié)果展示在表3中。從表3中可以看出我們的方法依舊有比較高的預(yù)測(cè)精確度，這也和隨機(jī)模擬實(shí)驗(yàn)的結(jié)果相符。我們方法的指標(biāo)相比NCL方法要好一些，并且變量選擇的數(shù)量上也比較相近。

表3 加法測(cè)量誤差情形下兩種方法應(yīng)用在乳腺癌數(shù)據(jù)中的結(jié)果Table 3 The results of two methods in breast cancer data under additive measurement error

5 結(jié)論

本文提出一種針對(duì)高維可加風(fēng)險(xiǎn)模型中帶有測(cè)量誤差情況下的變量選擇方法。在已知的生存分析數(shù)據(jù)相關(guān)文獻(xiàn)中，尚未有針對(duì)測(cè)量誤差數(shù)據(jù)的變量選擇方法。本文基于高維線性模型測(cè)量誤差數(shù)據(jù)的估計(jì)方法，重構(gòu)了高維可加風(fēng)險(xiǎn)模型，并給出了加法和乘法兩種測(cè)量誤差模型的變量選擇算法。簡(jiǎn)化偽得分方程的形式更加簡(jiǎn)潔且實(shí)用性強(qiáng)。隨機(jī)模擬實(shí)驗(yàn)和實(shí)際數(shù)據(jù)分析的相關(guān)結(jié)果證實(shí)了本文方法的有效性和精確性。

在未來(lái)的工作中，我們將致力于將簡(jiǎn)化偽得分方程應(yīng)用于高維可加風(fēng)險(xiǎn)模型的變量選擇中。同時(shí)也會(huì)對(duì)Cox模型，加速失效模型等其他生存分析模型中的測(cè)量誤差數(shù)據(jù)利用最近鄰半正定投影的方法進(jìn)行變量選擇方面的探索。