用SCAD懲罰研究帶遺傳約束的Cox模型

管 鋒, 晏海青

(1.鄭州大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,河南 鄭州 450001;2.光輝學(xué)校,河南 光山 465450;3.中國(guó)人民銀行羅山縣支行,河南 羅山 464299)

生存分析是對(duì)生存數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析的一門(mén)新興學(xué)科, 在醫(yī)學(xué)、生命科學(xué)和可靠性工程等領(lǐng)域都有廣泛應(yīng)用.Cox比例風(fēng)險(xiǎn)模型是由Cox在1972年提出, 是生存分析中應(yīng)用最多的模型[1].該模型是一個(gè)半?yún)?shù)模型, 不僅能應(yīng)用于完整的觀測(cè)數(shù)據(jù), 而且對(duì)于帶有刪失的數(shù)據(jù)同樣適用. 1975年Cox提出了偏似然的估計(jì)方法[2]來(lái)估計(jì)模型中的參數(shù). 由于這種方法估計(jì)的結(jié)果不具有稀疏性, 一些常用的變量選擇方法, 如最佳子集選擇法, 逐步選擇、計(jì)算量大且缺乏穩(wěn)定性. Zhang等[3]將自適應(yīng)LASSO估計(jì)方法應(yīng)用到Cox比例風(fēng)險(xiǎn)模型中, 得到的估計(jì)量具有良好的Oracle性質(zhì). Yang等[4]研究了超高維Cox可加模型的特征篩選, 并且通過(guò)Monte Carlo模擬和實(shí)例分析檢驗(yàn)了其可行性, 對(duì)Cox可加模型的應(yīng)用給予了理論支持.

在生物等研究領(lǐng)域中, 一些變量間的交互作用對(duì)響應(yīng)變量也有影響. 自適應(yīng)LASSO的方法研究了帶遺傳約束的Cox模型, SCAD懲罰估計(jì)可以減少估計(jì)偏差并且估計(jì)量同時(shí)滿足無(wú)偏性、稀疏性和連續(xù)性[5].因此自適應(yīng)LASSO方法研究帶遺傳約束的Cox模型[6]的基礎(chǔ)上,用SCAD懲罰估計(jì)方法研究帶遺傳約束的Cox模型也具有重要意義.

1 模型與方法

根據(jù)Fan等[7]在2001年提出來(lái)的懲罰似然的變量選擇問(wèn)題和線性回歸模型中帶交互項(xiàng)的變量選擇問(wèn)題, 給出帶遺傳約束的Cox模型, 并得到滿足遺傳約束條件的目標(biāo)函數(shù). 運(yùn)用梯度下降法的思想并結(jié)合CCCP算法[8]得到新的參數(shù)估計(jì)算法.

1.1 模型提出

假設(shè)在一次生存研究中有n個(gè)相互獨(dú)立的個(gè)體, 第i個(gè)個(gè)體觀測(cè)到的生存時(shí)間為T(mén)i,i=1,…,n, 是獨(dú)立同分布的隨機(jī)變量, 令Yi(t)=I(Ti≥t), 用來(lái)表示第i個(gè)個(gè)體在t時(shí)刻的生存情況.X(i)=(Xi1,Xi2,…,Xip)τ,i=1,…,n, 表示第i個(gè)個(gè)體對(duì)應(yīng)的p維協(xié)變量;β= (β1,β1,…,βp)τ表示p維未知參數(shù)向量,Z(i)=(Xi1,…,Xip,(Xi1Xi2),…,(Xi,p-1Xip))τ, 表示含有交互項(xiàng)的協(xié)變量,α=(α12,…,αp-1,p)τ, 表示交互項(xiàng)的未知參數(shù)向量, 則帶交互項(xiàng)的Cox模型可寫(xiě)成

(1)

在有右刪失數(shù)據(jù)的情況下可以得到偏對(duì)數(shù)似然函數(shù)

(2)

為了使估計(jì)量滿足遺傳約束條件, 重新參數(shù)化交互項(xiàng)的系數(shù), 令αjk=γjkβjβk,j

(3)

這里通過(guò)對(duì)交互項(xiàng)參數(shù)的重新參數(shù)化, 使得當(dāng)βj=0或者βk=0時(shí), 必有含變量Xj或者Xk的交互項(xiàng)系數(shù)為零.應(yīng)用SCAD懲罰偏對(duì)數(shù)似然來(lái)估計(jì)參數(shù), 不僅可以選擇出主要項(xiàng)和交互項(xiàng)中對(duì)響應(yīng)變量有解釋作用的協(xié)變量, 而且還可以使得估計(jì)量具有Oracle性質(zhì), 同時(shí)滿足遺傳約束條件.

1.2 算法實(shí)現(xiàn)

(4)

(5)

(6)

其中θ*=(β*τ,α*τ)τ.通過(guò)最小化(5)式, 可以得到(4)式的最小解.

同理當(dāng)γjk固定時(shí), 目標(biāo)函數(shù)可寫(xiě)為

(7)

(8)

(9)

通過(guò)最小化 (9) 式, 可以得到 (7) 式的最小解.上述估計(jì)過(guò)程的具體算法如下:

懲罰參數(shù)λβ和λγ的選擇可參考文獻(xiàn)[9].

2 漸近性質(zhì)

其中,j,k,l=1,…,p,k

性質(zhì)1的證明類(lèi)似Fan and Li (2001)中定理1的證明可得,這里省略.

3 數(shù)值模擬

假設(shè)真實(shí)的模型參數(shù)為

θ=(1,-2,0,0,-0.2,-2,0,0,-0.2,0,0,0.4,0,0,0)τ,

其中前五項(xiàng)為主要項(xiàng)Xi1,…,Xi5, 后十項(xiàng)為交互項(xiàng)Xi1Xi2,…,Xi4Xi5.主要項(xiàng)服從均值為0向量的多元正態(tài)分布, 協(xié)方差滿足Cov(Xik,Xij)=0.5∣k-j∣,1≤k,j≤5,1≤i≤n, 交互項(xiàng)為對(duì)應(yīng)主要項(xiàng)的乘積.由于在實(shí)踐中得到的數(shù)據(jù)是刪失的, 所以根據(jù)Bender[10]提供的產(chǎn)生生存時(shí)間的方法模擬生成隨機(jī)右刪失生存時(shí)間數(shù)據(jù), 刪失時(shí)間分布漸近服從U[0,c],c由刪失率確定, 主要是研究n=100,200,400, 刪失率為 25%和 40%的有限樣本表現(xiàn).

表1是在刪失率為25%, 樣本量分別為100, 200和400的情況下重復(fù)估計(jì)500 次得到的偏差, 標(biāo)準(zhǔn)差和均方誤差.表2是在刪失率為40%情況下估計(jì)的結(jié)果, 與表1對(duì)比可以發(fā)現(xiàn), 在刪失率增加的情況下估計(jì)的精確度隨之減小, 隨著樣本量的增加, 估計(jì)的精確度也明顯增加, 與漸近無(wú)偏估計(jì)的理論結(jié)果相吻合.表3結(jié)果是在不同刪失比率和不同樣本量下對(duì)正確選擇出模型的評(píng)估, 從表中數(shù)據(jù)可以看出樣本量越大刪失率越小, 正確選出模型的概率越高.

表1 在25%刪失率下非零系數(shù)估計(jì)的bias, SD和MSE

表2 在40%刪失率下非零系數(shù)估計(jì)的bias, SD和MSE

表3 在不同刪失率下模型和零系數(shù)的正選率

4 總結(jié)

主要用SCAD懲罰估計(jì)方法研究帶遺傳約束的Cox比例風(fēng)險(xiǎn)模型的變量選擇與參數(shù)估計(jì)問(wèn)題, 利用SCAD懲罰估計(jì)方法對(duì)參數(shù)向量進(jìn)行懲罰, 在給定的正則條件下估計(jì)量具有Oracle性質(zhì).另外, 在許多實(shí)際問(wèn)題中, 主要項(xiàng)可能還不足以預(yù)測(cè)協(xié)變量和響應(yīng)變量之間的關(guān)系,而該研究方法可以解決這類(lèi)問(wèn)題.