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

管 鋒, 晏海青

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

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

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

1 模型與方法

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

1.1 模型提出

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

(1)

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

(2)

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

(3)

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

1.2 算法實現(xiàn)

(4)

(5)

(6)

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

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

(7)

(8)

(9)

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

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

2 漸近性質(zhì)

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

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

3 數(shù)值模擬

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

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

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

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

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

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

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

4 總結(jié)

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