張家睿，吳耀華

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

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

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

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

1 研究背景

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

1.1 高維可加風險模型

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

(1)

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

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

(2)

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

Lin和Ying[15]采用一種有如下形式的偽得分方程來對可加風險模型進行分析：

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

(3)

其中β∈p，并且

(4)

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

(5)

和

(6)

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

U0(β)=b0-V0β.

(7)

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

(8)

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

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

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

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

(9)

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

不失一般性，用Lasso算法來舉例說明測量誤差的影響，對于線性模型y=Xβ+來說，Lasso算法是最小化

(10)

這等價于最小化

(11)

(12)

然后解決下面的優(yōu)化問題來得到β的估計：

(13)

(14)

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

(15)

(16)

(17)

2 帶有測量誤差的高維可加風險模型的變量選擇方法

2.1 簡化偽得分方程

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

(18)

以及

(19)

則有

(20)

接著定義

(21)

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

U(β)=b-Vβ,

(22)

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

(23)

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

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

2.2.1 加法測量誤差

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

(24)

和

(25)

(26)

2.2.2 乘法測量誤差

(27)

以及

(28)

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

(29)

3 理論性質(zhì)

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

(30)

(31)

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

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

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

(32)

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

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

(33)

其中

(34)

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

(35)

則接下來有

引理3.2說明了乘法測量誤差的計算方法滿足近鄰條件。將引理3.1，引理3.2和定理3.1結合有

推論3.1給出了加法測量誤差估計方法和乘法測量誤差估計方法的理論保證，確定了估計量l1和l2的誤差界，下面將通過隨機模擬實驗和實際數(shù)據(jù)分析來驗證我們的理論結果。

4 實驗及結果分析

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

4.1 隨機模擬

4.1.1 加法測量誤差模型

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

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

(36)

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

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

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

4.1.2 乘法測量誤差模型

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

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

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

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

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

(37)

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

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

5 結論

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

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