肖 翔，吳懿祺，古 晞

（1.上海工程技術(shù)大學(xué) 數(shù)理與統(tǒng)計(jì)學(xué)院,上海 201620；2.上海對外經(jīng)貿(mào)大學(xué) 統(tǒng)計(jì)與信息學(xué)院,上海 201620；3.同濟(jì)大學(xué) 數(shù)學(xué)科學(xué)學(xué)院,上海 200092）

混合分布模型常用于對混合計(jì)數(shù)數(shù)據(jù)的研究[1?4].其中混合治愈模型在醫(yī)學(xué)上有著廣泛的應(yīng)用，它最早是在癌癥臨床研究中提出的，用于解決治愈率估計(jì)問題.目前，國內(nèi)外學(xué)者在對混合治愈模型的研究中取得豐富的成果.Boag[5]基于對數(shù)正態(tài)分布構(gòu)建易感群體的生存函數(shù)和混合治愈模型，對癌癥復(fù)發(fā)的可能性進(jìn)行預(yù)測.Farewell[6]利用威布爾分布建立關(guān)于治愈率的Logistic 回歸模型，用極大似然法估計(jì)參數(shù).Ghitany 等[7]引入?yún)f(xié)變量建立基于指數(shù)分布混合治愈模型，采用極大似然法估計(jì)參數(shù)的相合性和漸進(jìn)正態(tài)性.Taylor[8]使用Kaplan-Meier 估計(jì)量對易感群體的失效時(shí)間建模，引入?yún)f(xié)變量對治愈率建立回歸模型，并利用EM 算法得到參數(shù)的估計(jì)值.Peng 等[9]采用非參數(shù)估計(jì)方法建立基于比例風(fēng)險(xiǎn)的混合治愈模型，并對乳腺癌數(shù)據(jù)進(jìn)行實(shí)證分析.Zhou 等[10]使用多重插值補(bǔ)獲得治愈概率和未治愈的患者生存函數(shù)的參數(shù)和方差估計(jì).Diao 等[11]針對具有存活率的當(dāng)前狀態(tài)數(shù)據(jù)，提出一類半?yún)?shù)變化治愈模型，并證明模型回歸系數(shù)的極大似然估計(jì)具有一致性、漸進(jìn)正態(tài)性和漸進(jìn)有效性.Lam 等[12]提出聚類加權(quán)廣義估計(jì)方程方法，使用基于Bernstein 多項(xiàng)式的偽似然法估計(jì)模型的參數(shù)和非參數(shù)分量.

上述文獻(xiàn)一般采用傳統(tǒng)的或者基于EM 算法的極大似然估計(jì)法、多重插值補(bǔ)和偽似然法對治愈模型或者混合治愈模型中的參數(shù)進(jìn)行估計(jì)，很少采用Jeffreys 先驗(yàn)和reference 先驗(yàn)等客觀先驗(yàn)，這一領(lǐng)域研究幾乎空白.因此，本研究使用客觀先驗(yàn)對基于指數(shù)分布且含有刪失數(shù)據(jù)的混合治愈模型進(jìn)行客觀貝葉斯分析.

1 基于指數(shù)分布的混合治愈模型

假設(shè)被研究總體分成兩個(gè)群體：一類是治愈人群，其治愈率為p(0t)，T為研究對象的生存時(shí)間.一般意義下的混合治愈模型為

假設(shè)易感人群的生存時(shí)間服從指數(shù)分布，其概率密度為

式中：λ>0.對應(yīng)易感人群的生存函數(shù)為S0(t)=exp(?λt)，代入式(1)得到基于指數(shù)分布混合治愈模型為

假設(shè){(ti,di)}(i=1,2,···,n)是樣本容量為n的觀測值，參數(shù)(p,λ)的似然函數(shù)為

式中：f(ti;p,λ)為生存時(shí)間的概率密度函數(shù)；S(ti;p,λ)為生存函數(shù).

若生存時(shí)間為服從指數(shù)分布的混合治愈模型式(2)，則 (p,λ)的似然函數(shù)可以進(jìn)一步寫成

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

由于式(3)、式(4)比較復(fù)雜，不利于客觀先驗(yàn)的計(jì)算，因此引入隱變量Z=(z1,z2,···,zn). 當(dāng)zi=0時(shí)，個(gè)體i來源于易感群體；當(dāng)zi=1時(shí)，個(gè)體i來源于治愈群體，該群體占總體的比例為p.zi是服從治愈率為p的伯努利分布，即zi～Bernoulli(p)，且E(zi)=p，i=1,2,···,n.因此，基于擴(kuò)充數(shù)據(jù){(ti,di,zi)}得到參數(shù)(p,λ)的完全似然函數(shù)為

其對數(shù)完全似然函數(shù)為

2 客觀貝葉斯分析

2.1 Fisher 信息矩陣

由于原始的對數(shù)似然函數(shù)式(4)非常復(fù)雜，本節(jié)采用對數(shù)完全似然函數(shù)式(6)進(jìn)行推導(dǎo)，得到Fisher 信息矩陣是對角矩陣，極大地簡化了客觀先驗(yàn)的計(jì)算.

定理1假設(shè)對研究對象的觀測時(shí)間足夠長，則基于指數(shù)分布混合治愈模型式(2)，參數(shù) (p,λ)的Fisher 信息矩陣為

證明：對數(shù)完全似然函數(shù)式(6)的一階和二階偏導(dǎo)數(shù)為

記 為失效人數(shù)，假設(shè)對研究對象的觀測時(shí)間足夠長，這樣就可以識(shí)別出治愈者，通過樣本數(shù)據(jù)觀測到的未治愈率無限趨向于總體未治愈率 1?p，為后續(xù)計(jì)算方便，本研究取r≈n(1?p)，得到Fisher 信息陣的組成元素為

則在數(shù)據(jù)擴(kuò)充策略下，參數(shù) (p,λ)的Fisher 信息矩陣為

2.2 客觀先驗(yàn)

相比于Laplace 先驗(yàn)，Jeffreys 先驗(yàn)?zāi)軌蛟趨?shù)變換下保持不變性，比Laplace 先驗(yàn)具有更廣泛的應(yīng)用場合[13].

定理2參數(shù) (p,λ)的Jeffreys 先驗(yàn)為

證明：參數(shù)(p,λ)的Jeffreys先驗(yàn)與Fisher信息矩陣行列式的平方根成正比，且有

reference 先驗(yàn)是基于信息量準(zhǔn)則推導(dǎo)出的先驗(yàn)分布，能夠使參數(shù)的先驗(yàn)分布和后驗(yàn)分布之間的Kullback-Liebler 距離最大.reference 先驗(yàn)是Jeffreys 先驗(yàn)的推廣，當(dāng)參數(shù)是一維時(shí)，reference先驗(yàn)就變成了Jeffreys 先驗(yàn).

reference 先驗(yàn)根據(jù)參數(shù)的重要性，可得到不同重要次序下的參數(shù)組合.例如，記號(hào){p,λ}表示p為感興趣參數(shù)，λ為討厭參數(shù)，p比λ更重要.反之，記號(hào){λ,p}表示λ為感興趣參數(shù)，p為討厭參數(shù)，λ比p更重要.

定理3對于參數(shù)組合{p,λ}，(p,λ) reference 先驗(yàn)為

證 明：對于參數(shù)組合 {p,λ}，p為感興趣參數(shù)，F(xiàn)isher 信息矩陣I可寫為

式中：p?和 λ?為參數(shù)空間中事先給定的值.

定理4對于參數(shù)組合{λ,p}，(p,λ)的reference先驗(yàn)為

式中：p?和 λ?為參數(shù)空間中事先給定的值.

從定理3 和4 的結(jié)論可以看出，無論p為感興趣參數(shù)，還是 λ為感興趣參數(shù)，它們的reference 先驗(yàn)相同，為后續(xù)討論方便，記 πR=πR1= πR2.

2.3 后驗(yàn)分析

定理5基于先驗(yàn)分布 πJ和 πR，數(shù)據(jù)擴(kuò)充策略下(p,λ)的后驗(yàn)分布都是恰當(dāng)?shù)?

結(jié)合式(8)和(9)，可得

因此，基于先驗(yàn)分布πR，(p,λ)的后驗(yàn)分布πR(p,λ|ti,di,zi)是恰當(dāng)?shù)?，同理，基于先?yàn)分布πJ，(p,λ)的后驗(yàn)分布πJ(p,λ|ti,di,zi)也是恰當(dāng)?shù)?從而，保證了后驗(yàn)樣本抽樣的可行性，如果后驗(yàn)分布不是恰當(dāng)?shù)?，則無法抽樣得到后驗(yàn)樣本.

2.4 抽樣算法設(shè)計(jì)

基于完全似然函數(shù)式(5)，本節(jié)設(shè)計(jì)后驗(yàn)樣本的Gibbs 抽樣機(jī)制.

首先，設(shè)計(jì)條件分布，對隱變量Z=(z1,z2,···,zn)進(jìn)行抽樣，公式為

在給定Z=(z1,z2,···,zn)的條件下，從后驗(yàn)分布πR(p,λ|ti,di,zi)中分別抽取參數(shù)p和λ的后驗(yàn)樣本.從式(7)可以得到p的滿條件后驗(yàn)分布為

對logπR(p|λ,ti,di,zi)求關(guān)于p的二階導(dǎo)數(shù)

式(13)說明p的滿條件后驗(yàn)分布式(12)是對數(shù)上凸的，滿足自適應(yīng)拒絕抽樣法(ARS)的使用條件.

另外，從式(7)可以得到 λ的滿條件后驗(yàn)分布為

可見，λ的滿條件后驗(yàn)分布(14)服從形狀參數(shù)

最后，實(shí)施Gibbs 抽樣，具體步驟如下.

1) 設(shè)定參數(shù)初始值p(0),λ(0).

2) 對t=1,2,···，實(shí)施迭代更新:

(i)利用上一步的參數(shù)估計(jì)p(t?1),λ(t?1)及式(10)和式(11)，得到隱變量Z=(z1,z2,···,zn)的 樣本，i=1,2,···,n；

(ii) 對式(12)采用自適應(yīng)拒絕抽樣法(ARS)，抽樣得到樣本p(t)；

3 數(shù)值模擬

本節(jié)分別基于先驗(yàn)分布πJ和πR，對p和λ進(jìn)行貝葉斯估計(jì).設(shè)置3組不同的參數(shù)真值：1)p=0.3,λ=3；2)p=0.3，λ=4；3)p=0.4，λ=4，樣本容量分別設(shè)置為n=20和n=50.每次模擬均重復(fù)1000次，并從估計(jì)偏差(Bias)、均方誤差(RMSE)和95%覆蓋率(CP)3 個(gè)方面對估計(jì)效果進(jìn)行評價(jià)，見表1 至表3.

表1 客觀貝葉斯估計(jì)下的估計(jì)偏差Table 1 Bias under objective Bayesian estimation

表2 客觀貝葉斯估計(jì)下均方誤差Table 2 RMSE under objective Bayesian estimation

表3 客觀貝葉斯估計(jì)的95%覆蓋率Table 3 95% coverage probabilities under objective Bayesian estimation

從表中可以看出，隨著樣本容量的增加，客觀貝葉斯估計(jì)下的估計(jì)偏差與均方誤差都在降低，說明參數(shù)的點(diǎn)估計(jì)效果在逐步提高.同時(shí)，隨著樣本容量的增加，估計(jì)效果的差異在降低，說明客觀先驗(yàn)在樣本量小時(shí)優(yōu)勢比較明顯.從參數(shù)估計(jì)95%覆蓋率來看，兩種客觀貝葉斯先驗(yàn)下的估計(jì)結(jié)果都相對比較穩(wěn)定，基于reference 先驗(yàn)的估計(jì)效果比基于Jeffreys 先驗(yàn)的效果更穩(wěn)定，這是因?yàn)閞eference 先驗(yàn)比Jeffreys 先驗(yàn)更擅長處理多維參數(shù)的估計(jì).

4 結(jié)語

本研究討論了基于指數(shù)分布的混合治愈模型及其客觀貝葉斯分析方法，基于完全似然函數(shù)，寫出了參數(shù)的Fisher 信息矩陣，較為簡潔地計(jì)算出參數(shù)Jeffreys 先驗(yàn)和reference 先驗(yàn)，并驗(yàn)證了后驗(yàn)分布的恰當(dāng)性，這種方法可以推廣到更復(fù)雜的混合治愈模型.今后將引入?yún)f(xié)變量信息，建立關(guān)于治愈率的回歸模型，制定個(gè)性化診療策略.