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

崔俊峰

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

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

崔俊峰

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

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

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

1 引言

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

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

2 威布型模型及其Bayesian分析

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

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

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

設(shè) Ti為獨立同分布的隨機變量，具有分布函數(shù) F與密度函數(shù) f，Li為獨立同分布的隨機變量，具有分布函數(shù)G與密度函數(shù)g，且與 Ti獨立，由于截斷的發(fā)生，觀察到的僅是min（Ti，Li），令

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

假設(shè){Li}為第 i個產(chǎn)品的隨機截尾時間，序列{Ti}為第 i個產(chǎn)品的壽命，當壽命與其隨機截尾時間相互獨立時，分布函數(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ù)集,

由隨機截斷數(shù)據(jù)的似然函數(shù)式可知該模型的完備數(shù)據(jù)集似然函數(shù)為：

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

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

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

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

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

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

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

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

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

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

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

［1］劉樂平,張美英,李姣嬌.基于 WinBUGS軟件的貝葉斯計量經(jīng)濟學 ［J］.華東理工大學學報,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］,杭州：浙江科學技術(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］.北京：中國人民大學出版社,2004

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

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

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

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

［責任編輯：劉向紅］

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

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

崔俊峰（1978-），男，甘肅隴西人，隴南師范高等專科學校講師，研究方向：概率論及可靠性理論.