崔俊峰

（隴南師范高等?？茖W(xué)校，數(shù)信學(xué)院，甘肅 成縣 742500）

威布爾型共享脆弱模型的貝葉斯實(shí)證研究

崔俊峰

（隴南師范高等?？茖W(xué)校，數(shù)信學(xué)院，甘肅 成縣 742500）

本文主要借助 MCMC方法,利用 Gibbs抽樣對威布爾型共享脆弱模型進(jìn)行貝葉斯分析,利用數(shù)據(jù)模擬說明了模型的有效性.

威布爾型共享Frailty模型,Gibbs抽樣,Bayesian分析

1 引言

近20年來,脆弱模型在描述子組中個體壽命之間的相關(guān)性方面尤為活躍,國外踴現(xiàn)出了大量關(guān)于脆弱模型改進(jìn)和優(yōu)化的研究文獻(xiàn).20年前，人們經(jīng)常聽到的一句話是：“貝葉斯分析在理論上確實(shí)很完美,但遺憾的是在實(shí)際應(yīng)用過程中不能計(jì)算出結(jié)果”.隨著計(jì)算機(jī)技術(shù)的高速發(fā)展,再復(fù)雜的模型也可通過貝葉斯方法進(jìn)行處理.貝葉斯方法主要集中在后驗(yàn)期望的計(jì)算上，一種常見的貝葉斯計(jì)算類型是計(jì)算后驗(yàn)分布的眾數(shù).計(jì)算后驗(yàn)分布期望的傳統(tǒng)數(shù)值計(jì)算方法是數(shù)值積分、Laplace近似計(jì)算和 Monte-Carlo重要抽樣［1］.數(shù)值積分在中等維數(shù)(最大為10)問題上非常有效，而 MonteCarlo重要抽樣方法可以計(jì)算維數(shù)很大的問題，并且有著很高的計(jì)算精度.目前，MCMC方法已經(jīng)變成了非常流行的貝葉斯計(jì)算方法.一方面是由于它處理非常復(fù)雜問題的效率，另一方面是因?yàn)樗木幊谭椒ㄏ鄬θ菀祝?］.需要強(qiáng)調(diào)的是，并不是說 MCMC方法已經(jīng)完全替代傳統(tǒng)的方法，在一些特別的場合(如要求精度),傳統(tǒng)方法還具有它的優(yōu)勢.

共享脆弱模型是 COX比例風(fēng)險(xiǎn)模型［3］的擴(kuò)展,是脆弱模型中較為簡單的一種類型,因個體在同一組內(nèi)享有相同的異質(zhì)性而得名,下面在文獻(xiàn)［4］給出的部分定義基礎(chǔ)上討論這一模型.

2 威布型模型及其Bayesian分析

其中 i=1，2，…，n， j=1，2，…，mi，wi為第 i組中個體的脆弱性，h0（yij）為基準(zhǔn)風(fēng)險(xiǎn)率函數(shù)，xij為 p×1的伴隨變量，獨(dú)立于時(shí)間變量，代表了i組中影響第 j個個體壽命分布的主要因素,β為p×1的回歸系數(shù)，刻畫了伴隨變量影響程度的大小.

式(2-1)中，wi為來自均值為 0，方差通常未知的獨(dú)立同分布的樣本，和COX比例風(fēng)險(xiǎn)模型相比，我們可以清楚地看到，當(dāng)wi的取值大于 1時(shí)，該組中個體失效的速度會傾向于比獨(dú)立模型中以wi概率取1時(shí)的速度更快，反之，若 wi取值小于 1，該組中個體的預(yù)期存活時(shí)間將大于采用獨(dú)立模型的預(yù)期存活時(shí)間［5］.現(xiàn)有的文獻(xiàn)對wi服從的分布討論較多,尤其以單參數(shù)伽瑪分布更為常用，假設(shè) wi服從單參數(shù)的伽瑪分布,記為Γ（k-1，k-1），即有-1

設(shè) w=（w1，w2，…，wn）'，則 w的先驗(yàn)分布形式為：

設(shè) Ti為獨(dú)立同分布的隨機(jī)變量，具有分布函數(shù) F與密度函數(shù) f，Li為獨(dú)立同分布的隨機(jī)變量，具有分布函數(shù)G與密度函數(shù)g，且與 Ti獨(dú)立，由于截?cái)嗟陌l(fā)生，觀察到的僅是min（Ti，Li），令

因此似然函數(shù)［6］為：

假設(shè){Li}為第 i個產(chǎn)品的隨機(jī)截尾時(shí)間，序列{Ti}為第 i個產(chǎn)品的壽命，當(dāng)壽命與其隨機(jī)截尾時(shí)間相互獨(dú)立時(shí)，分布函數(shù) G中不含 F中的信息，因此似然函數(shù)可寫成

設(shè) X=（X1，X2，…，Xn），其中 Xi為第 i子組中m×p的協(xié)變量矩陣.δ=（δ11，δ12，…，δnmn）'為刪失示性變量集,其中 δij為第 i個子組中第 j個個體壽命的示性函數(shù).D=（X，δ，y，w）為完備數(shù)據(jù)集, Dobs=（X，δ，y）為觀測數(shù)據(jù)集,

由隨機(jī)截?cái)鄶?shù)據(jù)的似然函數(shù)式可知該模型的完備數(shù)據(jù)集似然函數(shù)為：

理論上講，可以通過對上式的積分獲得觀測數(shù)據(jù)集 Dobs=（X，δ，y）的似然函數(shù)，但此式復(fù)雜,利用積分方法計(jì)算（β，α，γ）的后驗(yàn)分布維數(shù)高、難度大［7］.因而對后驗(yàn)分布的計(jì)算借助MCMC方法,利用 Gibbs抽樣方法得出參數(shù)聯(lián)合后驗(yàn)分布的抽樣.令 π（·）表示有關(guān)參數(shù)的驗(yàn)前分布,則 wi（i=1，2，…，n）的條件密度為：

易證后驗(yàn)密度均是對數(shù)凹的.

3 模型可靠性應(yīng)用實(shí)證

本文仿真采用的數(shù)據(jù)為某統(tǒng)計(jì)部門測得的腎臟移植患者死亡時(shí)間數(shù)據(jù),來源于參考文獻(xiàn)［5］P8例 1.5,該文獻(xiàn)中討論了應(yīng)用核平滑方法估計(jì)危險(xiǎn)率的問題,并且特別討論了改變帶寬和核選擇所帶來的影響.本文對這一數(shù)據(jù)進(jìn)行了刪減,并利用脆弱模型在WinBUGS軟件包中對其進(jìn)行模擬.

假定患者的壽命服從二參數(shù)威布分布,為了考察年齡和性別對患者壽命的影響(即伴隨變量為年齡和性別,且假設(shè)患者僅受這兩種環(huán)境因素的影響),取樣本容量 N=38,在觀察周期內(nèi)記錄個體感染的時(shí)間 yij，將觀測結(jié)束時(shí)發(fā)生感染的時(shí)間(每個個體觀測了兩次感染的時(shí)間)視為有效數(shù)據(jù),觀測時(shí)尚未感染的個體被隨機(jī)截尾,只知道其感染時(shí)間不低于對其觀測的時(shí)間［9］.因而，模型中的比例風(fēng)險(xiǎn)因素

exp（x'ijβ）=exp（βsexsexi+βageageij）

sexi=1表示第 i個患者為女性,0為男性,ageij表示第 i（i=1，2，L…，38）個患者在第 j（j=1，2）次感染時(shí)的年齡.

基于參考文獻(xiàn)［5］P8的數(shù)據(jù),建立如(2-2)式所示的脆弱模型,在WinBUGS軟件包中進(jìn)行模擬,截取前 1000次迭代數(shù)據(jù)做為退火(burnin)以消除初值的影響,從第1001次開始進(jìn)行10000次迭代分析.

表3 .1 10000次抽樣迭代的參數(shù)后驗(yàn)估計(jì)統(tǒng)計(jì)量

由上表不難看出，α的均值為1.233，置信區(qū)間為 (0.9429,1.546)，βage的均值為 0.007157，置信區(qū)間為 (-0.01835,0.08268)，βsex的均值為-1.886，置信區(qū)間為(-3.019,-0.8577).

為了驗(yàn)證退火后模型的收斂性,在 WinBUGS中對相關(guān)參數(shù)進(jìn)行多條鏈迭代,即同時(shí)給出多組初值,經(jīng)過馬爾可夫鏈蒙特卡洛模擬形成多條迭代鏈,假如參數(shù)模型收斂(即符合馬爾可夫鏈蒙特卡洛方法穩(wěn)定條件時(shí))則迭代動態(tài)軌跡圖逐步趨于重合.本模型中輸入不同的初始值分別進(jìn)行迭代,收斂性診斷圖以及初始值的迭代鏈軌跡圖中均趨于重合,再次說明了模型的有效性.

［1］劉樂平,張美英,李姣嬌.基于 WinBUGS軟件的貝葉斯計(jì)量經(jīng)濟(jì)學(xué) ［J］.華東理工大學(xué)學(xué)報(bào),2007,6(2)：101-107.

［2］Berger,J.Bayesian Analysis：A Look at Today and Thoughts of Tomorrow［J］.Journal of the American Statistical Association.2000,95：1269-1276.

［3］黎子良,鄭祖康.生存分析［M］,杭州：浙江科學(xué)技術(shù)出版社,1993.

［4］Gupta,P.L.,Gupta,R.C.,Ageing classes of Weibull mixtures［J］.Probab.Engrg.Inform.Sci.1996,10：591-600.

［5］彭非,王偉.生存分析［M］.北京：中國人民大學(xué)出版社,2004

［6］王炳興.Weibull分布基于定數(shù)逐次截尾壽命數(shù)據(jù)的統(tǒng)計(jì)分析［J］.科技通報(bào)，2004,20(6)：488-490.

［7］崔俊峰，脆弱模型中變異系數(shù)的性質(zhì)及其應(yīng)用［J］.廣東技術(shù)師范學(xué)院學(xué)報(bào),2016.2：5-7.

［8］林靜.基于 MCMC的貝葉斯生存分析理論及其在可靠性評估中的應(yīng)用［博士學(xué)位論文］,江蘇：南京理工大學(xué),2008.

［9］李榮.基于刪失試驗(yàn)的貝葉斯生存回歸模型及其應(yīng)用［D］,湖南：湖南大學(xué),2006.

［責(zé)任編輯：劉向紅］

Bayesian Analysis of Weibull Shared-Frailty Model

CUI Junfeng
(The Department of Mathematics,Longnan Teachers’College,Chengxian Gansu 742500,China)

The Bayesian analysis of Weibull shared-Frailty model by using the MCMC method and Gibbs sampling was discussed.The validity of this model was showed by data simulation.

Weibull shared-Frailty model;Gibbs sampling;Bayesian analysis

O 212

A

1672-402X（2016）08-0011-03

2016-03-19

甘肅省教育科學(xué)“十二五”規(guī)劃課題：“Frailty模型及其可靠性應(yīng)用研究”（項(xiàng)目主持人：崔俊峰，項(xiàng)目編號：GS［2015］GHB0903）

崔俊峰（1978-），男，甘肅隴西人，隴南師范高等?？茖W(xué)校講師，研究方向：概率論及可靠性理論.