基于指數分布混合治愈模型的客觀貝葉斯分析

肖 翔，吳懿祺，古 晞

（1.上海工程技術大學 數理與統計學院,上海 201620；2.上海對外經貿大學 統計與信息學院,上海 201620；3.同濟大學 數學科學學院,上海 200092）

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

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

1 基于指數分布的混合治愈模型

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

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

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

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

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

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

其對數似然函數為

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

其對數完全似然函數為

2 客觀貝葉斯分析

2.1 Fisher 信息矩陣

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

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

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

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

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

2.2 客觀先驗

相比于Laplace 先驗，Jeffreys 先驗能夠在參數變換下保持不變性，比Laplace 先驗具有更廣泛的應用場合[13].

定理2參數 (p,λ)的Jeffreys 先驗為

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

reference 先驗是基于信息量準則推導出的先驗分布，能夠使參數的先驗分布和后驗分布之間的Kullback-Liebler 距離最大.reference 先驗是Jeffreys 先驗的推廣，當參數是一維時，reference先驗就變成了Jeffreys 先驗.

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

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

證 明：對于參數組合 {p,λ}，p為感興趣參數，Fisher 信息矩陣I可寫為

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

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

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

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

2.3 后驗分析

定理5基于先驗分布 πJ和 πR，數據擴充策略下(p,λ)的后驗分布都是恰當的.

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

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

2.4 抽樣算法設計

基于完全似然函數式(5)，本節設計后驗樣本的Gibbs 抽樣機制.

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

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

對logπR(p|λ,ti,di,zi)求關于p的二階導數

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

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

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

最后，實施Gibbs 抽樣，具體步驟如下.

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

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

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

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

3 數值模擬

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

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

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

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

從表中可以看出，隨著樣本容量的增加，客觀貝葉斯估計下的估計偏差與均方誤差都在降低，說明參數的點估計效果在逐步提高.同時，隨著樣本容量的增加，估計效果的差異在降低，說明客觀先驗在樣本量小時優勢比較明顯.從參數估計95%覆蓋率來看，兩種客觀貝葉斯先驗下的估計結果都相對比較穩定，基于reference 先驗的估計效果比基于Jeffreys 先驗的效果更穩定，這是因為reference 先驗比Jeffreys 先驗更擅長處理多維參數的估計.

4 結語

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