右刪失數(shù)據(jù)下比例風(fēng)險混合治愈模型的貝葉斯自適應(yīng)Lasso

王純杰，朱笑瑩，劉新蕊，羅琳琳

(長春工業(yè)大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，吉林 長春 130012)

0 引言

在生物醫(yī)學(xué)研究中，一般假定研究的個體都會經(jīng)歷感興趣的事件.然而，在實際生活中，通常存在部分個體是被治愈的，即延長隨訪時間，研究個體也不會經(jīng)歷感興趣的事件.為了解決此類情況，J.Berkson和R.P.Gage[1]提出混合治愈模型，該模型中研究對象包含治愈群體及非治愈群體兩部分；V.T.Farewell[2-3]首次對治愈率建立Logistic模型，對混合治愈模型中的潛伏期建立一個參數(shù)威布爾分布.針對混合治愈模型中潛伏期部分，即協(xié)變量對易感個體生存時間的影響，已經(jīng)提出了多種半?yún)?shù)模型，其中比例風(fēng)險混合治愈模型得到了廣泛的關(guān)注.在臨床醫(yī)學(xué)領(lǐng)域中，該模型得到了廣泛的應(yīng)用，如乳腺癌、頭頸癌、白血病、前列腺癌、黑色素瘤等[4-7].

在生存分析中經(jīng)常出現(xiàn)協(xié)變量較多的情況，同時也對各類治愈模型的參數(shù)估計帶來了困難[8-10].因此對混合治愈模型進(jìn)行變量選擇成為眾多學(xué)者關(guān)注的關(guān)鍵問題之一.目前已經(jīng)有諸多變量選擇的方法應(yīng)用于比例風(fēng)險混合治愈模型中.X.Liu等[11]采用SCAD和自適應(yīng)Lasso兩種方法對半?yún)?shù)混合治愈模型進(jìn)行變量選擇；田舒[12]考慮帶有聚類時間數(shù)據(jù)的比例風(fēng)險混合治愈模型，通過SCAD方法對該模型進(jìn)行變量選擇.以上研究都是基于頻率學(xué)派進(jìn)行的，而貝葉斯也是統(tǒng)計學(xué)中另一個重要的學(xué)派，在生物統(tǒng)計中應(yīng)用廣泛.目前在貝葉斯框架下提出了多種變量選擇的方法，例如：T.Park等[13]以及C.Hans[14]提出了貝葉斯框架下的Lasso方法.此后，貝葉斯自適應(yīng)Lasso的變量選擇法被廣泛應(yīng)用于生存分析中[13-16].然而針對右刪失數(shù)據(jù)下比例風(fēng)險混合治愈模型的貝葉斯自適應(yīng)Lasso的研究較少.

本文將采用貝葉斯自適應(yīng)Lasso方法對比例風(fēng)險混合治愈模型進(jìn)行變量選擇，采用三次樣條方法擬合基準(zhǔn)風(fēng)險函數(shù).建立貝葉斯框架下的層次模型，結(jié)合Gibbs及Metropolis-Hastings 算法對參數(shù)進(jìn)行抽樣.在不同設(shè)置下的有限樣本中驗證方法的有效性，最后應(yīng)用于一個乳腺癌的實例中.

1 數(shù)據(jù)和模型

在臨床研究中，由于某些外界因素的影響，會導(dǎo)致部分樣本的生存時間不能被準(zhǔn)確觀測到，只知道感興趣的時間發(fā)生在某一時間點之后，此時數(shù)據(jù)發(fā)生右刪失.考慮n個獨立的右刪失個體，其數(shù)據(jù)結(jié)構(gòu)可表示為

為了更好地描述混合治愈模型，引入治愈指標(biāo)Yi(i=1，2，…，n)，Yi表示第i個個體是否被治愈.當(dāng)Yi=1時表示第i個個體未被治愈，否則個體被治愈.令

π(x)=P(Y=1|X=x).

表示個體易感概率.受文獻(xiàn)[2]啟發(fā)，本文對個體易感概率構(gòu)建Logistic模型，即

其中β=(β1，β2，…，βs)為s維回歸系數(shù)向量.

令S(t|Y=1)表示未被治愈患者的條件生存函數(shù)，S(t|Y=0)表示治愈個體的條件生存函數(shù)，當(dāng)研究個體被治愈時，生存概率可以表示為

S(t|Y=0)=P(T>t|Y=0)=1.

當(dāng)研究對象未被治愈時，本文將對其生存時間建立比例風(fēng)險模型[19].其形式為

λ(t|Y=1，x)=λ0(t|Y=1)exp(zTγ)，

其中λ0(t|Y=1)為條件基準(zhǔn)風(fēng)險函數(shù)，γ=(γ1，…，γp)為p維的系數(shù)向量.因此可得到易感群體的條件生存函數(shù)為

S(t|Y=1)=e-Λ0(t)exp(zTγ)，

S(t)=π(x)S(t|Y=1)+1-π(x)=π(x)e-Λ0(t)ezTγ+1-π(x).

此時似然函數(shù)為

其中累積基準(zhǔn)風(fēng)險函數(shù)Λ0(t)是完全未指定的非降函數(shù).本文采用樣條的方法對基準(zhǔn)風(fēng)險函數(shù)λ0(t)建立三次樣條模型

從而得到Λ0(t)為

其中：(t-κm)+=max(0，t-κm)，節(jié)點個數(shù)m是固定的，并且κ1，κ2，…，κm是預(yù)先選定的.令h=(h0，h1，h2，h3，h31，…，h3m)，得到似然函數(shù)的變換形式為

則對數(shù)似然函數(shù)為

其中：

2 貝葉斯自適應(yīng)Lasso

本文采用J.Fan等[21]提出的最小化懲罰似然函數(shù)法，對比例風(fēng)險混合治愈模型進(jìn)行變量選擇.基于該模型，參數(shù)的貝葉斯自適應(yīng)Lasso定義為

其中參數(shù)λj≥0，j=1，2，…，p.

本文給定待估參數(shù)γ一個條件拉普拉斯先驗，形式為

令θ=(βT，γT，hT)T基于右刪失數(shù)據(jù)，由層次模型可得后驗分布為

選取合適的MCMC算法，從上述后驗分布中抽取隨機(jī)數(shù)，即可進(jìn)行后驗推斷.

3 后驗推斷

根據(jù)上述后驗分布，可求得各個參數(shù)的滿條件后驗.由于參數(shù)β，γ，h對應(yīng)的滿條件后驗分布不是常用的分布類型，因此本文使用Metropolis-Hastings算法對其進(jìn)行參數(shù)抽樣.除此之外，其余各個參數(shù)的滿條件后驗分布為常用的分布類型.因此本文采用Gibbs算法對其進(jìn)行抽樣.抽樣算法步驟如下：

⑤α-2j(m)的滿條件后驗分布為逆高斯分布.從下式中進(jìn)行抽樣：

其中j=1，2，…，q.

⑥參數(shù)τ-2j(m)的滿條件后驗分布是逆高斯分布.可以通過下式進(jìn)行抽樣：

其中j=1，2，…，p.

⑦υ-2j(m)的滿條件后驗分布為逆高斯分布.從下式中進(jìn)行抽樣：

其中j=1，2，…，s.

⑧參數(shù)β(m)的聯(lián)合后驗分布為

β(m)滿條件后驗分布函數(shù)不是常用的分布類型，所以本文采用的是Metropolis-Hastings算法對抽取β(m)，其算法步驟如下：

ⅰ.從提議分布中選出一個候選點β(prop)，本文選取的是正態(tài)分布g(β(prop)|β(m-1))=N(β(prop)|β(m-1)，Σ)，其中Σ為常數(shù)矩陣.

ⅱ.從均勻分布U(0，1)中抽取隨機(jī)數(shù)u.令接受概率r為

判斷u和r的大小，如果u

⑨參數(shù)h(m)的滿條件后驗分布為

同理，h(m)滿條件后驗分布函數(shù)不是常用的分布類型，因此仍采用Metropolis-Hastings算法對h(m)進(jìn)行抽樣估計，步驟同上.

⑩參數(shù)γ(m)的滿條件后驗分布為

同理，γ(m)滿條件后驗分布函數(shù)不是常用的分布類型，因此仍采用Metropolis-Hastings算法對γ(m)進(jìn)行抽樣估計，步驟同上.

(3)重復(fù)步驟(2)，直至所有參數(shù)的馬爾科夫鏈?zhǔn)諗?

4 數(shù)值模擬

在本節(jié)中，通過模擬研究評估所提出的貝葉斯自適應(yīng)Lasso方法在變量選擇和參數(shù)估計方面的性能.考慮如下比例風(fēng)險混合治愈模型：

S(t)=π(x)e-Λ0(t)ezTγ+1-π(x)，

其中基線危險函數(shù)設(shè)置為λ0(t|Y=1)=2，系數(shù)向量γ=(0.8，0，1，0，1，0，1，0)T，β=(1，0.9)T.對協(xié)變量設(shè)置考慮兩種情況：

設(shè)置1Z為正態(tài)分布N(0，1)中產(chǎn)生的8維協(xié)變量.

設(shè)置2Z1，Z2，…，Z4為正態(tài)分布N(0，1)中產(chǎn)生，Z5，Z6，Z7從獨立的二項分布B(1，0.5)中產(chǎn)生，Z8從泊松分布P(1)中產(chǎn)生.易感概率從Logistic模型中產(chǎn)生，即

設(shè)定協(xié)變量X=(1，X0)T，其中X0來自正態(tài)分布N(0，0.5).

刪失時間來自于的均勻分布U(0，uc)里產(chǎn)生，產(chǎn)生右刪失比例為50%，治愈比例為30%的有限樣本.選取λ先驗分布中的超參數(shù)aj0=0.1，bj0=0.05.使用具有4個等距節(jié)點的三次樣條估計λ0(t)和Λ0(t)，即q=8.

在模擬設(shè)置中考慮樣本量分別為500、1 000、2 000的數(shù)據(jù).MCMC算法的鏈長設(shè)置為15 000，退火5 000.模擬過程重復(fù)了1 000次.使用偏差的平均值(BIAS)、標(biāo)準(zhǔn)誤差(SEE)、均方根誤差(RMSE)以及95%置信區(qū)間覆蓋率(CP)四種指標(biāo)對貝葉斯自適應(yīng)Lasso方法進(jìn)行評估.為了判斷系數(shù)γ是否顯著，設(shè)定臨界值為0.15.如果參數(shù)滿足|γ|<0.15，那么就認(rèn)為其對應(yīng)的系數(shù)就是不顯著的.模擬結(jié)果見表1—2.

表1 協(xié)變量來自正態(tài)分布的模擬結(jié)果

表2 協(xié)變量來自正態(tài)分布、二項分布及泊松分布的模擬結(jié)果

由表1可以看出，貝葉斯自適應(yīng)Lasso方法下的BIAS、SSE、RMS的值都較小，CP值都接近于0.95.與此同時，隨著樣本量的增大，各參數(shù)的大部分偏差都相應(yīng)的減小，標(biāo)準(zhǔn)誤差及均方誤差也越來越小.同時，所對應(yīng)的95%置信區(qū)間覆蓋率越來越接近0.95.除此之外，參數(shù)γ選擇的準(zhǔn)確性也越來越高，當(dāng)樣本量達(dá)到2 000時，參數(shù)γ選擇的準(zhǔn)確性達(dá)到100%.因此通過模擬研究得出隨著樣本量的增大，參數(shù)估計結(jié)果和變量選擇的準(zhǔn)情性都隨之變好.

通過表2中的數(shù)據(jù)可得如下結(jié)論：隨著樣本量不斷增大，參數(shù)γ及β的BIAS、SSE、RMS的值都隨之減小，CP值也越來越接近0.95.并且，參數(shù)γ選擇的準(zhǔn)確性隨著樣本量的增大也逐漸升高.綜合表1—2，本文提出的貝葉斯自適應(yīng)Lasso方法在協(xié)變量是連續(xù)分布和離散分布情況下表現(xiàn)都較好.

5 實例應(yīng)用

將貝葉斯自適應(yīng)Lasso方法應(yīng)用于一個實例中.該數(shù)據(jù)來自于686個乳腺癌患者，本文對可能造成乳腺癌疾病的因素進(jìn)行篩選，目的是選擇出影響該疾病的8個協(xié)變量.該數(shù)據(jù)的8個可能影響乳腺癌疾病的變量為：horTh(是否進(jìn)行激素治療)、age(患者的年齡/年)、menostat(更年期狀態(tài)，其中pre表示更年期之前，post表示更年期之后)、tsize(腫瘤大小/mm)、tgrade(腫瘤水平因子，水平1<水平2<水平3)、pnodes：(正節(jié)點數(shù)/個)、progrec(孕酮受體的個數(shù)/個)、ester(雌激素受體的個數(shù)/個).采用貝葉斯自適應(yīng)Lasso方法對影響乳腺癌發(fā)生的主要因素進(jìn)行篩選，抽取長度為15 000的馬爾科夫鏈，退火5 000.其結(jié)果見表3.

表3 乳腺癌數(shù)據(jù)變量選擇的結(jié)果

設(shè)置閾值為0.05，從表3可以得出，變量7及變量8對應(yīng)的估計值小于閾值，所以被剔除.變量2、變量4及變量6與乳腺癌的發(fā)生成正相關(guān).變量1及變量3與乳腺癌的發(fā)生成正相關(guān).即患者年齡越大患乳腺癌的風(fēng)險就越大，患者體內(nèi)的腫瘤越大就越可能惡化，同時，腫瘤水平越高，最后患乳腺癌的可能性越大，正節(jié)點數(shù)越多，越容易發(fā)病.患者可以通過激素治療來降低患病的風(fēng)險，沒發(fā)生更年期的群體不易患乳腺癌，與實際情況相符.

6 結(jié)語

本文利用自適應(yīng)Lasso的方法，對比例風(fēng)險混合治愈模型進(jìn)行參數(shù)估計及變量選擇.在貝葉斯的框架下，通過給各個參數(shù)不同的先驗構(gòu)造層次模型，進(jìn)而求出聯(lián)合后驗分布，利用Gibbs及Metropolis-Hastings 算法對各個參數(shù)進(jìn)行抽樣，得到了平穩(wěn)的馬爾科夫鏈.模擬研究表明，在不同樣本及協(xié)變量分布不同的情況下，針對比例風(fēng)險混合治愈模型的參數(shù)估計及變量選擇問題，所采用的貝葉斯自適應(yīng)Lasso效果較好.最后將該方法應(yīng)用到實例數(shù)據(jù)中，以解決實際問題.

然而本文僅討論了右刪失數(shù)據(jù)下比例風(fēng)險混合治愈模型的貝葉斯自適應(yīng)Lasso方法的可行性，在其他復(fù)雜數(shù)據(jù)類型下需進(jìn)一步討論.