廣義指數(shù)分布下區(qū)間刪失數(shù)據(jù)貝葉斯回歸分析

李 群, 董小剛, 王純杰, 趙 波

(長(zhǎng)春工業(yè)大學(xué) 基礎(chǔ)科學(xué)學(xué)院， 吉林 長(zhǎng)春 130012)

廣義指數(shù)分布下區(qū)間刪失數(shù)據(jù)貝葉斯回歸分析

李 群, 董小剛, 王純杰*, 趙 波

(長(zhǎng)春工業(yè)大學(xué) 基礎(chǔ)科學(xué)學(xué)院， 吉林 長(zhǎng)春 130012)

研究了在兩參數(shù)廣義指數(shù)分布下的區(qū)間刪失壽命時(shí)間的貝葉斯回歸分析模型。生存時(shí)間在服從廣義指數(shù)分布的情況下，假定形狀參數(shù)的先驗(yàn)分布來(lái)自伽馬分布，建立了尺度參數(shù)與生存時(shí)間貝葉斯回歸模型，從而得到生存時(shí)間的變化。選取MCMC算法對(duì)參數(shù)進(jìn)行估計(jì)，并運(yùn)用R軟件進(jìn)行了模擬。

廣義指數(shù)分布; 區(qū)間刪失; 貝葉斯回歸; MCMC算法

0 引 言

在可靠性壽命試驗(yàn)中，兩參數(shù)廣義指數(shù)分布可簡(jiǎn)稱廣義指數(shù)分布或GE分布。作為指數(shù)分布的推廣，由于廣義指數(shù)分布對(duì)于刪失時(shí)間數(shù)據(jù)有很好的分析效果，而且還可以作Gamma分布和Weibull分布的替代分布，因而在壽命試驗(yàn)和可靠性工程中有著重要的應(yīng)用[1-2]。壽命數(shù)據(jù)分析已經(jīng)成為航空、工程、醫(yī)學(xué)和生物科學(xué)等多個(gè)領(lǐng)域中統(tǒng)計(jì)學(xué)家和實(shí)際工作者十分關(guān)心的一個(gè)問(wèn)題，因此，對(duì)廣義指數(shù)分布的研究有著十分重要的實(shí)際意義。同時(shí)，在生存分析中也經(jīng)常研究感興趣的時(shí)間與哪些因素有密切的關(guān)系，也會(huì)研究不同的藥物類型中，哪種藥物對(duì)于患者更有效果等。文中將通過(guò)建立貝葉斯回歸模型來(lái)進(jìn)行研究感興趣的時(shí)間與相關(guān)因素的關(guān)系及影響[3]。

李榮[4]于2006年給出了一篇?jiǎng)h失實(shí)驗(yàn)壽命的貝葉斯威布爾生存回歸模型，建立了威布爾分布下關(guān)于參數(shù)λ的回歸模型，并給回歸系數(shù)賦予先驗(yàn)分布。在刪失壽命實(shí)驗(yàn)的條件下，給出了貝葉斯威布爾回歸模型的似然函數(shù)，基于Gibbs抽樣得出參數(shù)的后驗(yàn)分布，利用WinBUGS軟件包求解威布爾回歸模型的貝葉斯估計(jì)的過(guò)程。朱惠明[5]等于2007年給出了刪失試驗(yàn)壽命的貝葉斯生存極值回歸模型，同樣引入?yún)?shù)λ的協(xié)變量,并建立了貝葉斯回歸模型，用MCMC方法和Gibbs抽樣獲得參數(shù)后驗(yàn)分布，同樣利用WinBUGS軟件包求解極值回歸模型的貝葉斯估計(jì)的過(guò)程。Upadhyay[6]發(fā)表了基于Gibbs抽樣下對(duì)數(shù)正態(tài)回歸的后驗(yàn)分析，分別建立對(duì)數(shù)正態(tài)分布的均值、方差兩個(gè)參數(shù)關(guān)于協(xié)變量影響的貝葉斯回歸模型。Puja Makkar[7]給出了頭頸癌在對(duì)數(shù)正態(tài)模型下的貝葉斯生存分析，在不知道先驗(yàn)信息的情況下，采用Gibbs抽樣的方法得到參數(shù)的后驗(yàn)分布，并分析不同的治療方案對(duì)患頭頸癌患者壽命的影響。

廣義指數(shù)分布是由Gupta R D和Kundu D于1999年提出的。此外Gupta R D[8-9]等給出了廣義指數(shù)分布的一些統(tǒng)計(jì)推斷的性質(zhì)。Kundu D[10]等于2008年給出了廣義指數(shù)的貝葉斯估計(jì)的相關(guān)理論。此外，郭環(huán)[11]研究了兩參數(shù)廣義指數(shù)分布的一些參數(shù)估計(jì)方法和優(yōu)良性質(zhì)，給出了在一定條件下兩個(gè)參數(shù)的貝葉斯估計(jì)。但是上述文獻(xiàn)均沒(méi)有涉及廣義指數(shù)分布的貝葉斯生存回歸模型。

MCMC方法是一種簡(jiǎn)單易行、廣泛應(yīng)用的計(jì)算隨機(jī)模擬方法。該方法的核心思想是構(gòu)建一個(gè)概率轉(zhuǎn)移矩陣，建立一個(gè)以分布π(x)為平穩(wěn)分布的Markov鏈，得到π(x)的樣本，通過(guò)隨機(jī)抽樣得到的樣本就可以進(jìn)行各種統(tǒng)計(jì)推斷和估計(jì)[12]。MCMC方法中最常用的一種方法是Metropolis-Hastings，該方法最早由Metropolis于1953年給出的，后來(lái)Metropolis的算法由Hastings改進(jìn)，合稱為M-H算法[13-14]。M-H算法是MCMC的基礎(chǔ)方法，由M-H算法演化出了許多新的抽樣方法，包括目前在MCMC中最常用的Gibbs抽樣也可以看做M-H算法的一個(gè)特例[15]。

1 符號(hào)描述和模型介紹

文中主要研究的是區(qū)間刪失下的廣義指數(shù)分布模型的建立及貝葉斯回歸分析的應(yīng)用。下面假設(shè)第i個(gè)個(gè)體滿足以下關(guān)系：

假定每個(gè)個(gè)體都可以觀測(cè)兩次，其中U、V代表兩個(gè)隨機(jī)變量，并且以概率1滿足U

文中采用的是廣義指數(shù)分布對(duì)區(qū)間刪失數(shù)據(jù)進(jìn)行建模[8-9]。廣義指數(shù)分布的密度函數(shù)為：

(1)

其分布函數(shù)為：

(2)

生存函數(shù)為:

(3)

風(fēng)險(xiǎn)函數(shù)為：

(4)

式中： α----形狀參數(shù);

λ----尺度參數(shù)。

形狀參數(shù)為α，尺度參數(shù)為λ的廣義指數(shù)分布記為GE(α,λ)。

其對(duì)應(yīng)的全數(shù)據(jù)的似然函數(shù)為：

(5)

文中研究的是區(qū)間刪失情況下的貝葉斯回歸模型，則區(qū)間刪失情況下的似然函數(shù)為：

(6)

故區(qū)間刪失數(shù)據(jù)的對(duì)數(shù)似然函數(shù)可以表示為：

(7)

接下來(lái)建立貝葉斯層次模型如下：

(8)

α～gamma(1,0.001)

其中,λi指每個(gè)個(gè)體生存時(shí)間所服從的廣義指數(shù)分布的尺度參數(shù)，βj,j=0,1,…,m的先驗(yàn)分布為正態(tài)分布，α的先驗(yàn)分布為gamma分布。

這樣就可以建立起區(qū)間刪失數(shù)據(jù)的廣義指數(shù)分布貝葉斯回歸模型。接下來(lái)可根據(jù)貝葉斯層次模型寫出后驗(yàn)的聯(lián)合密度函數(shù),即后驗(yàn)似然函數(shù)[3,16]為：

(9)

故得到后驗(yàn)對(duì)數(shù)似然函數(shù)為：

(10)

接著，運(yùn)用MCMC算法對(duì)參數(shù)進(jìn)行估計(jì)。

2 數(shù)值模擬

用數(shù)值模擬過(guò)程來(lái)評(píng)價(jià)文中建立的模型性能，給出模擬步驟如下：

1)產(chǎn)生來(lái)自均勻分布U[-2,2]的N個(gè)獨(dú)立同分布的x1,x2,…,xN。

2)設(shè)定β0=1，β1=1，α=1.5，并令λi=exp(β0+β1xi)。

3)產(chǎn)生N個(gè)服從廣義指數(shù)分布的失效時(shí)間T,形狀參數(shù)α=1.5，尺度參數(shù)λi=exp(β0+β1xi)。

4)產(chǎn)生N個(gè)服從參數(shù)為θ1=6的指數(shù)分布的第一次觀測(cè)時(shí)間U,產(chǎn)生N個(gè)服從參數(shù)為θ2=0.2指數(shù)分布的第二次觀測(cè)時(shí)間V,并滿足U

5)比較U、V和失效時(shí)間T的大小關(guān)系，若TV，則令δ3=1，否則δ3=0。令δ2=1-δ1-δ3。

6)給出β和α的先驗(yàn)分布。并寫出先驗(yàn)似然函數(shù)(LL)和后驗(yàn)似然函數(shù)(LP)。

7)應(yīng)用MCMC算法估計(jì)參數(shù)β和α。

按照上述算法步驟，循環(huán)500次計(jì)算出待估參數(shù)β和α的均值和方差。樣本量設(shè)定為N分別為200、300、500，模擬結(jié)果見(jiàn)表1。

表1 樣本量為200,300,500的模擬結(jié)果

由表1 可以看出，在樣本量不同,且左刪失比例約為0.2，右刪失比例約為0.4的情況下，模擬參數(shù)的估計(jì)值較真值偏差較小，能夠給出對(duì)應(yīng)參數(shù)較好的估計(jì)結(jié)果，并且精度會(huì)隨著樣本量的增加而增加，樣本標(biāo)準(zhǔn)差也會(huì)隨著樣本量的增加而減小。由此可見(jiàn)，該模型用來(lái)進(jìn)行后驗(yàn)估計(jì)是可行的。在算法的選擇上也可采用其他的算法進(jìn)行估計(jì)。

3 實(shí)證分析

對(duì)一個(gè)實(shí)際數(shù)據(jù)例子進(jìn)行研究分析，選取的數(shù)據(jù)是1976年到1980年之間在波士頓進(jìn)行乳腺癌早期治療的回顧性研究數(shù)據(jù)。該數(shù)據(jù)由Finkelstein和Wolfe在1985年展現(xiàn)出來(lái)，數(shù)據(jù)是由94位病人組成，其中分為給予放射性治療組(RT)和放射性療法加輔助性化學(xué)治療組(RCT)。放射治療組共計(jì)46位病人，放療加輔助化療組共有48位病人[17]。

在研究過(guò)程中，病人每4～6個(gè)月隨訪一次，然而，病人的實(shí)際訪問(wèn)時(shí)間不同，每個(gè)病人的兩次隨訪時(shí)間也是不同的。在就診過(guò)程中醫(yī)生會(huì)根據(jù)乳腺收縮程度來(lái)評(píng)估病人情況。這項(xiàng)研究的目的是為了比較這兩組治療方式對(duì)患者的治療效果，看放療輔助化療方法是否可以提高患者的無(wú)復(fù)發(fā)率和總的生存率。但是有一些實(shí)驗(yàn)和臨床證據(jù)表明，化療加劇了正常組織對(duì)放射治療的急性反應(yīng)。這個(gè)數(shù)據(jù)包含了關(guān)于乳腺收縮的信息，但是沒(méi)有精確的觀測(cè)時(shí)間。這里有38例患者在研究期內(nèi)沒(méi)有明顯的乳腺收縮，所以這部分觀測(cè)設(shè)定為右刪失數(shù)據(jù)，即這個(gè)區(qū)間觀測(cè)沒(méi)有右側(cè)端點(diǎn)。對(duì)于其他患者，觀測(cè)時(shí)間的時(shí)間間隔代表著在這段時(shí)間內(nèi)發(fā)生過(guò)乳腺收縮。觀測(cè)時(shí)間的左端點(diǎn)是從第一次診所就診時(shí)間開始，到最后一次就診時(shí)發(fā)現(xiàn)乳腺收縮截止。例如，觀測(cè)到的(6,10]表示在第6個(gè)月隨訪時(shí)患者未出現(xiàn)乳腺收縮，但是在下一次隨訪，即第10個(gè)月時(shí)，患者出現(xiàn)了乳腺收縮。乳腺收縮情況出現(xiàn)在第6個(gè)月至第10個(gè)月兩次隨訪之間，但精確的時(shí)間未知。所以我們用區(qū)間的刪失時(shí)間數(shù)據(jù)來(lái)描述乳腺收縮。將這組數(shù)據(jù)進(jìn)行詳細(xì)地分析估計(jì)，觀測(cè)數(shù)據(jù)見(jiàn)表2。

在進(jìn)行數(shù)據(jù)分析的過(guò)程中，若第i個(gè)病人屬于放射治療組，定義協(xié)變量xi=0；若第i個(gè)病人屬于放療輔助化療組，定義協(xié)變量xi=1，并且假定乳腺癌發(fā)作時(shí)間服從廣義指數(shù)分布。估計(jì)結(jié)果見(jiàn)表3。

通過(guò)表3的實(shí)驗(yàn)結(jié)果可以求出

λ=exp(-15.873x)

可以判斷出：當(dāng)病人屬于放射治療組時(shí)，λ=1；當(dāng)病人屬于放療輔助化療組時(shí)，0<λ<1。從而根據(jù)生存函數(shù)可以判斷出,放療輔助化療方法可以提高患者的無(wú)復(fù)發(fā)率和總的生存率。

表2 乳腺癌觀測(cè)數(shù)據(jù)

表3 乳腺癌數(shù)據(jù)估計(jì)結(jié)果

4 結(jié) 語(yǔ)

在貝葉斯框架下建立了服從廣義指數(shù)分布的生存時(shí)間的尺度參數(shù)同影響生存時(shí)間的相關(guān)因素之間的回歸模型，并給出后驗(yàn)似然函數(shù)，采用MCMC方法對(duì)后驗(yàn)似然函數(shù)進(jìn)行求解最大值，同時(shí)解出了待估參數(shù)。并對(duì)該模型進(jìn)行了模擬，模擬效果較好。并將該方法應(yīng)用到乳腺癌數(shù)據(jù)例子中，結(jié)果表明，放射療法輔助化療方法對(duì)于提高患者的總的生存率有著一定的效果。

[1]GuptaRD,KunduD.Exponentiatedexponentialfamily;analternativetogammaandweibulldistributions[J].BiometricalJournal,2001,43(1):117-130.

[2] Gupta R D, Kundu D. Generalized exponential distribution: different method of estimations[J]. Journal of Statistical Computation and Simulation,2001,69(4):315-337.

[3] Ibrahim J G, Chen M H, Sinha D. Bayesian survival analysis[M]. New York: John Wiley & Sons Ltd.,2005.

[4] 李榮,朱慧明.刪失試驗(yàn)壽命的貝葉斯威布爾生存回歸模型[J].統(tǒng)計(jì)與決策,2006(24):20-22.

[5] 朱慧明,李榮,方博文.刪失試驗(yàn)壽命的貝葉斯生存極值回歸模型[J].系統(tǒng)工程與電子技術(shù),2007,29(11):1988-1990.

[6] Upadhyay S K, Peshwani M. Posterior analysis of lognormal regression models using the Gibbs sampler[J]. Statist. Papers,2008,49:59-85.

[7] Puja Makkar, Puneet K, Srivastava R S Singh, et al. Bayesian survival analysis of head and neck cancer data using lognormal model[J]. Communication in Statistics-Theory and Methods,2014,43(2):392-407.

[8] Gupta R D, Kundu D. Generalized exponential distributions[J]. Australian and New Zealand Journal of Statistics,1999,41(2):173-188.

[9] Gupta R D, Kundu D. Generalized exponential distributions: statistical inferences[J]. Technical Report, The University of New Brunswick, Saint John.,1999,41(3):111-115.

[10] Kundu D, Gupta R D. Generalized exponential distribution: bayesian estimations[J]. Computational Statistics & Data Analysis,2008,52(4):1873-1883.

[11] 郭環(huán).兩參數(shù)廣義指數(shù)分布的參數(shù)估計(jì)與數(shù)值模擬[D].武漢：華中科技大學(xué),2013.

[12] 黃小艷.MCMC方法分析[J].中國(guó)市場(chǎng),2015(14):185-186.

[13] Metropolis N, Rosenbluth A W, Rosenbluth M N, et al. Equation of state calculations by fast computing machines[J]. The Journal of Chemical Physics,1953,21(6):1087-1092.

[14] Hastings W K. Monte carlo sampling methods using markov chains and their applications[J]. Biometrika,1970,57(1):97-109.

[15] Geman S, Geman D. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence,1984(6):721-741.

[16] Bernardo J, Smith A. Bayesian theory[M]. West Sussex: John Wiley & Sons.,2000.

[17] Finkelstein D M, Ra W. A semiparametric model for regression analysis of interval-censored failure time data[J]. Biometrics,1985,41(4):933-945.

Bayesian survival regression analysis of interval censored data with generalized exponential Model

LI Qun, DONG Xiaogang, WANG Chunjie*, ZHAO Bo

(School of Basic Sciences, Changchun University of Technology, Changchun 130012， China)

Bayesian regression analysis model of interval censored lifetime under two-parameter Generalized Exponential is studied. Provided that the lifetime comes from generalized exponential distribution, and the prior distribution of shape parameter derives from the gamma distribution, the Bayesian regression model influenced by scale parameter and survival time is established to obtain the variation of lifetime. MCMC algorithm is used to estimate the parameters, and R software is used for simulation.

generalized exponential distribution; interval censored; bayesian regression; MCMC algorithm.

2016-07-19

國(guó)家自然科學(xué)基金青年基金項(xiàng)目(11301037)； 國(guó)家自然科學(xué)基金資助項(xiàng)目(11571051)； 吉林省教育廳“十三五”規(guī)劃項(xiàng)目(2016317)

李 群(1991-)，女，漢族，山東菏澤人，長(zhǎng)春工業(yè)大學(xué)碩士研究生，主要從事生存分析方向研究,E-mail:liqun91@live.com. *通訊作者：王純杰(1978-)，女，漢族，遼寧遼陽(yáng)人，長(zhǎng)春工業(yè)大學(xué)副教授，博士，主要從事數(shù)理統(tǒng)計(jì)、生存分析方向研究,E-mail:wangchunjie@ccut.edu.cn.

10.15923/j.cnki.cn22-1382/t.2016.6.16

O 212.4

A

1674-1374(2016)06-0597-06