刪失數(shù)據(jù)下部分線性模型的貝葉斯P-樣條估計

王純杰，羅琳琳，李純凈，袁曉惠

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

0 引言

近年來，卵巢癌的發(fā)病率逐年上升.由于早期卵巢癌不易察覺，故其對女性的身體健康產(chǎn)生了巨大的潛在威脅.由于其發(fā)病原因并不明確，因此，對卵巢癌數(shù)據(jù)的研究具有重大的實用價值.2009年，Cho和Shih[1]對卵巢癌進(jìn)行了細(xì)致描述并給出了針對卵巢癌的一些可用模型；2012年，李陽敏等[2]就早期卵巢癌的診斷數(shù)據(jù)建立多元線性回歸模型，以期找到早期卵巢癌的病理特征；2018年，赫艷玲等[3]對卵巢癌各階段的治療方案進(jìn)行了闡述.

線性回歸模型假定響應(yīng)變量與協(xié)變量的關(guān)系是線性的，在實際研究中，如果數(shù)據(jù)的協(xié)變量與響應(yīng)變量的關(guān)系為非線性時，線性回歸模型的擬合效果較差，于是便產(chǎn)生了易解釋且具有更好擬合效果的部分線性模型.部分線性模型最初由Engle等[4]提出，用于研究天氣對電力的影響.隨后很多學(xué)者對此類模型進(jìn)行了拓展研究.1995年，黃四民和梁華[5]運用部分線性模型分析居民消費結(jié)構(gòu)；2015年，Müller和Geer[6]將部分線性模型拓展到高維情形；2018年，Wu[7]在部分線性模型的DBE(Difference-Based Estimator)估計的基礎(chǔ)上進(jìn)行改進(jìn)，得到了估計效果更好的DBRE(Difference-Based Ridge Estimator)估計.

在醫(yī)療領(lǐng)域，學(xué)者們將部分線性模型與生存分析中的刪失數(shù)據(jù)相結(jié)合，并對此進(jìn)行了大量研究.2002年，陳敏和朱力行[8]給出了隨機(jī)刪失數(shù)據(jù)下的部分線性模型的擬合優(yōu)度檢驗；2018年，Wang[9]研究了右刪失數(shù)據(jù)下部分線性模型的經(jīng)驗似然推斷.

相比于頻率學(xué)派，貝葉斯充分利用了先驗的信息，在小樣本情況下得到了更好的估計效果.2004年，Lang和Brezger[10]首次針對AM(Additive Model)模型提出貝葉斯P-樣條，隨后此估計方法得到了廣泛的運用；2006年，Brezger和Lang[11]將貝葉斯P-樣條應(yīng)用于GAM(Generalized Additive Model)模型；2016年，Bruno[12]等人將貝葉斯P-樣條應(yīng)用于完整數(shù)據(jù)的部分線性模型；2018年，William和Gholamreza[13]將貝葉斯P-樣條應(yīng)用于RSM(Risk-Sharing Model)模型.

本文基于貝葉斯P-樣條，針對右刪失數(shù)據(jù)，運用部分線性模型分析一些因素對卵巢癌患者生存時間的影響，得到了有效的參數(shù)估計，為卵巢癌的防治提供了有效建議.

1 部分線性模型及貝葉斯P-樣條

部分線性模型最初由Engle等[4]提出，由于其可以處理非線性關(guān)系，因此具有很高的靈活性，其模型如下：

log(T)=g(Z)+αTX+ε.

其中：T為響應(yīng)變量；X為p維變量；Z為一維變量；g(·)是未知函數(shù)；α為p維變量X的系數(shù)；ε為隨機(jī)誤差，滿足Eε=0，Eε2=σ2，其分布可取標(biāo)準(zhǔn)正態(tài)分布(Normal)、標(biāo)準(zhǔn)極值分布(Extreme)和標(biāo)準(zhǔn)邏輯斯諦分布(Logistic).

在生存分析中很難得到完整的數(shù)據(jù)，所以需要處理不完整數(shù)據(jù)，其中最為常見的為右刪失數(shù)據(jù).令事件發(fā)生的確切時間為T，但一些病人在實驗結(jié)束后要觀察的事件仍然沒有發(fā)生，只知道事件發(fā)生的時間大于刪失點C.用δ=I(T≤C)表示數(shù)據(jù)是否刪失，故得到Y(jié)=min(T，C)的觀測數(shù)據(jù)，其數(shù)據(jù)結(jié)構(gòu)為D={Y，δ，X，Z}.

部分線性模型中較為重要的部分在于如何估計非參數(shù)部分g(Z).自從部分線性模型提出后，出現(xiàn)了許多關(guān)于g(Z)的估計方法，如樣條估計、核估計、M估計等.在這些方法中較為成熟的估計方法是B-樣條估計.在B-樣條中假定協(xié)變量的定義域為[a，b]，將協(xié)變量分成m段，a=k0

其中B(Z)為樣條基函數(shù)，其定義為

在B-樣條中，當(dāng)m很大時，容易出現(xiàn)過擬合現(xiàn)象.為了確保其有效性，1996年Eilers和Marx[14]提出了P-樣條，即在每個樣條基函數(shù)的系數(shù)上定義一個粗略的懲罰項以保證擬合出的曲線充分光滑.根據(jù)Eilers和Marx[14]的文章得到的懲罰似然函數(shù)如下：

其中：Δtβ為β的t階差分；f(yi|Xi，Z，β，α)為在給定Xi，Z，β的條件下Yi的密度函數(shù)；F(Yi|Xi，Z，β，α) 為給定Xi，Z，β的條件下Yi的分布函數(shù).在標(biāo)準(zhǔn)正態(tài)分布、標(biāo)準(zhǔn)極值分布和標(biāo)準(zhǔn)邏輯斯諦分布下的似然函數(shù)如下：

P-樣條在處理大樣本問題中具有很高的有效性，且相比于B-樣條而言，節(jié)點的選擇對估計的影響較小.但樣本量不夠多時，P-樣條的估計效果并不理想，于是Lang和Brezger[10]在P-樣條的基礎(chǔ)上加入了貝葉斯，得到的貝葉斯P-樣條方法在樣本量較小時仍有很好的估計效果.根據(jù)Lang和Brezger[10]的文章，對系數(shù)β取如下二階差分：

βj=2βj-1-βj-2+ωj.

其中ωj～N[0，τ/ψ]，τ用于保證函數(shù)的平滑度，ψ用于保證每一段上的曲率是不同的.根據(jù)Song等[17]的研究結(jié)果，非線性部分有如下約束：

故在貝葉斯部分，對β取有約束的高斯先驗，各參數(shù)的先驗分布如下：

在模擬與實例中假定樣條的段數(shù)m=18，差分階數(shù)t=2，參考文獻(xiàn)[10]，取α1=1，α2=0.005，v=1，α0=0.5，?=1.由此得出參數(shù)的后驗如下：

τ-1～Gamma(α1+(m-t)/2，α2+βTM(ψ)β/2)，

ψ～Gamma((v+m-t)/2，v/2+(βj-2βj-1+βj-2)/2τ)，

2 MCMC算法

(1) 給定初值α0，β0.

(3) 從U(0，1)中抽取u1t，若

則αt=α(new)，否則αt=αt-1.

β(new)=β(e)-B(BTB)-1BTβ(e).

(5) 從U(0，1)中抽取u2t，若

則βt+1=β(new)，否則βt+1=βt.

重復(fù)循環(huán)10 000次，為去除初始值的影響去掉開始的5 000次抽樣結(jié)果.剩余5 000次循環(huán)結(jié)果的均值即為β和α的估計值.

3 模擬研究

通過模擬檢驗貝葉斯P-樣條在刪失數(shù)據(jù)部分線性模型下的估計效果.部分線性模型的具體形式如下：

log(T)=g(Z)+αX+ε.

其中：g(z)=arctan(z)，Z～U(0，1)，α=2，X為成功概率為0.5的伯努利分布隨機(jī)數(shù).對誤差ε取3種不同的分布(見表1).

表1 誤差分布

令刪失比為40%，算得3種情況下的刪失點分別為C=αX+g(Z)+U[0，0.15]，C=αX+g(Z)+U[0，0.02]，C=αX+g(Z)+U[0，0.1].考慮不同樣本量的右刪失數(shù)據(jù).用來衡量α的估計好壞的指標(biāo)為偏和均方誤差.具體計算結(jié)果見表2.

表2 兩種不同方法的估計效果比較

由表2的結(jié)果可以看到，貝葉斯P-樣條方法在不同的誤差分布下α均擁有較小的bias與MSE，且樣本量越大，bias與MSE的值越小.說明在不同的誤差分布下部分線性模型的參數(shù)部分?jǐn)M合效果很好.而用B-樣條方法時，當(dāng)樣本量較大時的估計效果與貝葉斯P-樣條結(jié)果類似，但在樣本量較小時貝葉斯的估計效果要優(yōu)于B-樣條的方法，說明貝葉斯方法要優(yōu)于普通的B-樣條方法.取不同誤差分布的貝葉斯P-樣條非參數(shù)的g(Z)部分的擬合結(jié)果如圖1所示.圖2為B-樣條方法得出的曲線圖.對比圖1—2可以看出，貝葉斯P-樣條的擬合圖像具有較好的擬合效果，而B-樣條在前端與真實曲線相差較大，且貝葉斯P-樣條估計出的曲線更為平滑.由此可以看出貝葉斯P-樣條得到的曲線更光滑且樣條節(jié)點的選擇對其影響更小.

圖1 貝葉斯P-樣條下g(z)的擬合圖

4 實證研究

本文將部分線性模型的貝葉斯P-樣條估計應(yīng)用于一個醫(yī)學(xué)問題中.數(shù)據(jù)來自于Edmunson等[16]的研究，該研究記錄了對患卵巢癌的女性進(jìn)行兩種不同的治療(化療與放療加化療)的療效，結(jié)合患者的自身因素，分析了不同的治療方案以及患者自身因素對卵巢癌患者生存時間的影響.該研究共收集了26例卵巢癌病人的有效數(shù)據(jù)，刪失比為46%.選取4個影響因素：患者的年齡X1，患者對治療的耐受能力X2，治療方案X3，患者是否有其他疾病X4.對生存時間的對數(shù)log(T)采用部分線性模型建模，非參數(shù)部分運用貝葉斯P-樣條的方法進(jìn)行擬合.

圖2 B-樣條下g(z)的擬合圖

先對4個變量進(jìn)行篩選，用生存時間的對數(shù)對4個協(xié)變量建立加速失效模型，結(jié)果見表3.由于發(fā)現(xiàn)有些變量的P值大于0.05，故進(jìn)行逐步回歸篩選變量.逐步回歸結(jié)果見表4，逐步回歸過程剔除了患者對治療的耐受能力X2與患者是否有其他疾病X4，保留的患者的年齡與治療方案均有較小的P值，故接下來用這兩個變量建模.

表3 AFT模型估計結(jié)果

表4 逐步回歸后估計結(jié)果

患者年齡與生存時間對數(shù)圖像見圖3.由圖3可知，患者的年齡與生存時間的對數(shù)之間存在非線性趨勢，故令Z為患者的年齡，X為治療方案，建立如下部分線性模型：

log(T)=g(Z)+αX+ε.

由于誤差ε的分布未知，故分別選取標(biāo)準(zhǔn)正態(tài)分布，標(biāo)準(zhǔn)極值分布與標(biāo)準(zhǔn)邏輯斯諦分布進(jìn)行擬合.線性部分治療方案的系數(shù)估計結(jié)果見表5.

表5 線性部分治療方案系數(shù)的估計

非線性部分的擬合結(jié)果見圖4.

圖4 g(·)的擬合圖

在有可能影響卵巢癌患者生存時間的4個因素中，年齡與治療時采取的治療方案對卵巢癌患者的生存時間影響顯著.對誤差分布取3種不同的形式建立部分線性模型.其中治療方案對卵巢癌患者生存時間的影響是線性的，且呈正相關(guān)，即化療的治療方法比化療加放療的治療方法好，可更好地延長患者的壽命.年齡對卵巢癌患者生存時間的影響呈現(xiàn)非線性關(guān)系，這也符合客觀規(guī)律.由圖4可知60歲以下卵巢癌患者的生存時間與年齡成正比，而60歲以上的卵巢癌患者生存時間與年齡成反比.

5 總結(jié)

本文運用部分線性模型的貝葉斯P-樣條估計方法分析了卵巢癌數(shù)據(jù)，在模擬中貝葉斯P-樣條方法擬合的未知曲線更為光滑且其具有較小的bias和MSE.因此貝葉斯P-樣條方法可用在部分線性模型中對未知曲線進(jìn)行擬合.

卵巢癌大多發(fā)生在40歲以上的女性身上，且未婚未育女性發(fā)生卵巢癌的幾率更大.卵巢癌若早發(fā)現(xiàn)，早治療，痊愈的概率很大，但老年人由于年齡的增長，身體各項機(jī)能減弱，患病后很難痊愈，所以在實例中得出的結(jié)論是正確的.即對卵巢癌患者生存時間影響顯著的兩個變量是治療方案與患者年齡.其中患者年齡與卵巢癌患者的生存時間呈現(xiàn)非線性關(guān)系：當(dāng)年齡小于60歲時，年齡與生存時間成正比；當(dāng)年齡大于60歲時，年齡與生存時間成反比.