二元相依生存數(shù)據(jù)的貝葉斯統(tǒng)計診斷

佘思琪，唐安民

(云南大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，云南 昆明 650500)

在生存分析中，同一個個體（比如病人或設(shè)備）可能會有2 種類型的事件發(fā)生，即我們可以同時觀測到同一個體的2 個事件時間.在很多研究中假設(shè)個體的這2 個生存過程彼此獨立，則可以使用Kaplan-Meier乘積限估計方法[1]估計這2 個生存函數(shù).而實際情況中，這2 個生存過程往往不是獨立的，比如一個人2個不同器官的衰竭時間、癌癥復(fù)發(fā)時間和死亡的時間、第一次和第二次感染的時間等.忽略這種相依關(guān)系，而作為2 個獨立的生存過程來研究，得到的結(jié)果往往會是有偏的，甚至是錯誤的[2].

為了解決這個問題，在相依二元生存分析的研究中，較為流行的方法是Oakes[3]在1989 年提出的脆弱模型（frailty model），在給定脆弱性變量下，各個邊際生存時間條件獨立，從而實現(xiàn)對生存數(shù)據(jù)相關(guān)性建模.自從Sklar 提出Sklar 定理[4]后，Copula 函數(shù)受到了廣泛的關(guān)注，Copula 函數(shù)具有相對簡單的形式，可以基于各種參數(shù)或非參數(shù)模型的邊際分布建立依賴結(jié)構(gòu)，Hougaard[5]、Gustafson 等[6]、Romeo 等[7]、Muhammed[8]、Almetwally 等[9]都分別將Copula 函數(shù)應(yīng)用到二元相依生存分析研究中.然而，在二元相依生存數(shù)據(jù)模型統(tǒng)計診斷方面的國內(nèi)外研究仍然相對較少，盡管Cho 等[10]、Suzuki 等[11]和Lynch 等[12]均基于Kullback-Leibler 散度(K-L 散度)提出了相依二元生存數(shù)據(jù)的統(tǒng)計診斷方法，但邊緣分布都是基于全參數(shù)模型，全參數(shù)模型的假定可能偏離實際，導(dǎo)致有偏且不穩(wěn)健的參數(shù)估計.

本研究主要考慮使用FGM Copula 函數(shù)[13]度量二元生存過程之間的相依性，B 樣條擬合邊際對數(shù)基本危險函數(shù)，從而實現(xiàn)二元相依生存分析的半?yún)?shù)建模.在貝葉斯理論框架下，基于Gibbs 抽樣和MH（Metropolis-Hastings）抽樣的混合算法，對所擬合的二元生存模型進(jìn)行參數(shù)估計，并提出統(tǒng)計診斷方法，包括數(shù)據(jù)刪除影響分析和局部影響分析.

1 基于FGM Copula 的相依二元生存半?yún)?shù)模型

本節(jié)主要在貝葉斯理論框架下，假定事件的生存時間和刪失時間相互獨立的條件下，對具有相依關(guān)系的二元生存數(shù)據(jù)建立基于FGM Copula 的相依二元生存半?yún)?shù)模型.

1.1 模型建立根據(jù)Sklar 定理，一個Copula 函數(shù)將能夠?qū)⒉煌呺H分布連接起來，構(gòu)造成一個多元聯(lián)合分布函數(shù)，Copula 的參數(shù)則描述了變量之間的相關(guān)性.設(shè)二元聯(lián)合分布函數(shù)H(x,y) 的邊際分布為F(x) 和G(y)，則存在一個Copula 函數(shù)C，對任意的x,y∈[?∞,+∞]，有

如果F(x) 和G(y) 連續(xù)，則C是唯一的，并且該式很容易擴展到多元情況：H(x1,x2,···,xk)=C(F1(x1),F2(x2),···,Fk(xk)).

在二元分析中，對于任意的u,v∈[0,1]，參數(shù)為 φ 的Copula 分布函數(shù)表示為C(u,v|φ)，對應(yīng)的密度函數(shù)為c(u,v|φ)=條件分布函數(shù)分別表示為：常用的Copula 函數(shù)有Gauss Copula，Clayton Copula，F(xiàn)rank Copula，Gumbel Copula，F(xiàn)GM Copula 等，不同的Copula 刻畫了不同的相關(guān)關(guān)系，根據(jù)Nelson[14]對Copula 的定義，將這些常用Copula 函數(shù)的形式和參數(shù)取值范圍列舉如表1 所示.

表1 常用Copula 函數(shù)形式和參數(shù)取值范圍Tab.1 Common Copula function and parameter range

由于FGM Copula 具有簡單的結(jié)構(gòu)形式，有利于減輕計算負(fù)擔(dān)，其參數(shù)取值范圍為[?1,1]，可以同時度量生存過程之間的正相依性和負(fù)相依性，因而本文考慮使用FGM Copula建立二元生存時間之間的相關(guān)關(guān)系.

考慮n個獨立個體，令和Ci j分別是第i個個體在第j個生存過程的真實時間和刪失時間，其中i=1,2,···,n，j=1,2，本文考慮右刪失情形，即只能觀測到事件時間并引入示性指標(biāo)當(dāng) δij=1 時，觀測到真實生存時間，當(dāng) δij=0 時，觀測到刪失時間，而不是真實生存時間，但知道真實生存時間大于刪失時間.令Sij(tij) 和fij(tij) 分別表示第i個個體在第j個生存過程中在tij的邊際生存函數(shù)和邊際密度函數(shù)，則基于FGM Copula，可得 (Ti1,Ti2) 的聯(lián)合生存函數(shù)和聯(lián)合密度函數(shù)如下：

其中，參數(shù) φ 衡量了生存時間之間相依性的強度，并且當(dāng) φ=0 時，Ti1和Ti2獨立，即有Si(ti1,ti2|φ)=Si1(ti1)Si2(ti2).

接下來，基于邊際危險函數(shù) λi j(ti j) 對包含在（2）式中的邊際生存函數(shù)Sij(tij) 和邊際密度函數(shù)fij(ti j) 建模.通常地，邊際危險函數(shù) λi j(ti j) 可以建模如下：

其中 βj是協(xié)變量Xi=(xi1,···,xip)T系數(shù)向量，非負(fù)函數(shù) λ0j(ti j) 表示第j個生存過程中不考慮個體差異的基本危險函數(shù)，則邊際生存函數(shù)為：

由(3)，(4)式易得第i個個體的第j個邊際密度函數(shù)為fij(tij)=Sij(tij)λij(tij).在(4)式中未知的基本危險函數(shù) λ0j(ti j)，用B 樣條[15]建模其對數(shù)基本危險函數(shù)如下：

這里的k=1,2,···,qj+Kj+1.

由于(4)式涉及到復(fù)雜的積分運算，本文基于矩形近似積分方法提出一種高效計算方法[17]，將第j個生存過程中的事件時間軸用Lj+1 個時刻分為首尾相連但不相交的Lj個小區(qū)間：(c0j,c1j],(c1j,c2j],···,并且有c0j=0，則結(jié)合(4)，(5)式，邊際生存函數(shù)Sij(tij)可以表示如下：

顯然，當(dāng)tij＜cl?1,j，Iijl=0，該方法可以高效計算第i個個體的第j個生存過程的邊際生存函數(shù)和邊際密度函數(shù)

令Di=表示第i個個體觀測到的事件時間和刪失指標(biāo)，D=(D1,···,Dn)，在假定事件的生存時間和刪失時間相互獨立及個體間相互獨立的條件下，得到參數(shù)的聯(lián)合似然函數(shù)為：

其中第i個個體的似然函數(shù)Li(Di|φ,β1,β2,φ1,φ2) 可以表示如下：

其中Ri表示與參數(shù) φ 有關(guān)的項，具體表示如下：

其中 θ[?φ]表示 θ 中除了參數(shù) φ 以外的參數(shù)集合，表示類似，顯然各個參數(shù)的條件后驗不是常見的概率分布形式，不能直接從常見的分布形式中得到，本文借助Gibbs 抽樣和MH 抽樣的混合算法，對各個參數(shù)（向量）進(jìn)行估計，步驟如下：

第二步，對于第s(s=1,···,S) 步迭代，各個參數(shù)值迭代更新如下：

并從均勻分布U(0,1)中產(chǎn)生隨機數(shù)uφ，如果uφ≤αφ，則取φ(s)=φ*，否則φ(s)=φ(s?1).

第三步，不斷迭代第二步，直到產(chǎn)生的馬爾科夫鏈依分布收斂到目標(biāo)分布.

2 統(tǒng)計診斷

統(tǒng)計診斷是數(shù)據(jù)分析中的重要組成部分，本節(jié)在貝葉斯框架下，提出一種適合二元相依生存模型的數(shù)據(jù)刪除影響分析和局部影響分析方法，檢測出數(shù)據(jù)中影響模型的潛在異常點或強影響點.

2.1 貝葉斯數(shù)據(jù)刪除影響分析在貝葉斯分析中，基于數(shù)據(jù)刪除的K-L 距離是對有特定數(shù)據(jù)點和沒有特定數(shù)據(jù)點的后驗分布之間差異的度量.π(θ|D) 表示給定所有個體數(shù)據(jù)下 θ 的條件后驗概率密度，π(θ|D[?i])表示刪除第i個個體數(shù)據(jù)后 θ 的條件后驗概率密度，基于K-L 距離考慮如下貝葉斯診斷測度：

該值度量了條件后驗概率 π(θ|D) 和 π(θ|D[?i]) 的距離，較大的KL(i) 值意味著刪除第i個個體的數(shù)據(jù)可能對參數(shù) θ 后驗概率的計算產(chǎn)生較大影響.(11)式需要計算 π(θ|D[?i]) 的值，如果每刪除一個個體就計算一次刪除該個體數(shù)據(jù)后 θ 的條件后驗概率密度，計算量將不堪重負(fù)，還好，在個體間相互獨立假定下，可以證明 π(θ|D) 和 π(θ|D[?i]) 有如下關(guān)系：

其中Li(θ) 表示第i個個體的似然函數(shù)值，因此(11)式可以變換為：(11)式在求期望時涉及到難以精確計算的積分，因此采用蒙特卡洛積分可近似如下：

其中 θ(h) 表示從 π(θ|D) 抽取的第h次觀測值.

2.2 貝葉斯局部影響分析局部影響有助于識別數(shù)據(jù)對模型中的微小影響[18]，在生存數(shù)據(jù)的貝葉斯分析中，局部影響分析被發(fā)展為一類敏感性分析方法，基于Zhu 等[19]的局部影響方法，本文發(fā)展了二元相依生存模型的貝葉斯局部影響診斷方法.令 ω=(ω1,···,ωW) 表示擾動向量，表示擾動后的模型，刻畫了對原假定模型 M={p(D,θ)} 的不同擾動，滿足合適的擾動模型中，總存在擾動向量使得即這一擾動向量對模型沒有擾動.在貝葉斯框架下，擾動模型主要靠W×W維Fisher 信息矩陣來刻畫，記作和 ωk分別為 ω 的第w個分量和第k個分量，令則有：

其中Eω表示對p(D,θ|ω) 求期望，gww(ω) 可以解釋為擾動 ωw的擾動程度，gwk(ω) 表示 ωw和 ωk的關(guān)聯(lián)，因此刻畫了擾動分量 ωw和 ωk的相關(guān)程度，值越接近于1，認(rèn)為 ωw和 ωk在擾動模型中作用越相似.合適的擾動模型總要求 G(ω0) 是對角矩陣，當(dāng) G(ω0) 不是對角矩陣時，我們總能選擇一個新的擾動向量=ω0+G(ω0)12(ω?ω0)，并且

考慮對 ω 和 ω0的貝葉斯因子取對數(shù)，有：B(ω)=?(ω)??(ω0)，令 ω(t) 表示擾動模型 M={p(D,θ|ω)},ω ∈RW的擾動向量，取 ω(t)=ω0+th，并且 ω(0) 表示無擾動，得到一階Taylor 展開式：B(ω(t))=B(ω(0))+其中在沿著h的方向上定義一階局部影響測度為：

其中 FIB,h的最大值等價于的最大特征值，記作hmax，該值量化了擾動模型中 ω 局部影響的最大程度，并與最大特征值的特征向量相對應(yīng)，揭示了如何擾動模型從而得到貝葉斯因子的最大局部改變，進(jìn)而識別出潛在的影響觀測.在某些平滑條件下，可以得到：

其中為 G(ω0) 表示在 ω0處的Fisher 矩陣.計算關(guān)于B(ω) 的局部影響度量，我們僅需要計算 ?B和 G(ω0)，可以用蒙特卡洛方法近似.在我們所提出的相依二元生存模型中，根據(jù)Ibrahim 等[20]的建議，考慮對第w個個體觀測值的生存過程進(jìn)行如下的擾動：

3 數(shù)值模擬與實例分析

本節(jié)基于數(shù)值模擬和實例數(shù)據(jù)分析，驗證我們建議模型的可行性和統(tǒng)計診斷方法能否以較大概率找出潛在的異常點或強影響點.

3.1 數(shù)值模擬產(chǎn)生300 個生存數(shù)據(jù)的樣本，其中FGM Copula 中相關(guān)系數(shù) φ=0.8，2 個邊際分布的基本危險函數(shù)假定如下：1+0.5cos(1.5πt) 和 1+0.7cos(2.5πt).為了模型可識別，不考慮截距項，分別從參數(shù)為0.5的二項分布和2 個標(biāo)準(zhǔn)正態(tài)分布獨立產(chǎn)生3 個協(xié)變量xi1,xi2,xi3，i=1,···,300，未知參數(shù)向量 β1=(1.5,1,1)T和 β2=(2,?1,1)T，如此設(shè)置使得刪失率在20%～25%，樣條次數(shù)設(shè)置為2 并取5 個內(nèi)部節(jié)點.進(jìn)行積分近似運算時，劃分100 個小積分區(qū)間.為了檢驗馬爾科夫鏈的收斂性，從不同的初值出發(fā)產(chǎn)生3 條平行的馬爾科夫鏈，每條鏈迭代10 000 次，通過計算感興趣參數(shù)（φ,βj,j=1,2）的EPSR（Estimated Potential Scale Reduction）值的方法[21]來判斷生成的鏈?zhǔn)欠袷諗?，如果迭代大? 000 次時所有參數(shù)的EPSR 值都小于1.2，則認(rèn)為5 000 次以后的馬爾科夫鏈已經(jīng)收斂了.為了判斷基于B 樣條建模的基本生存函數(shù)是否收斂，在其定義域中均勻插入20 個點，以這20 個點擬合的基本生存函數(shù)擬合值作為“虛擬參數(shù)”值，計算這些“虛擬參數(shù)”的EPSR 值，與感興趣參數(shù)的EPSR 值一起作圖，得到圖1.

圖1 數(shù)值模擬中參數(shù)EPSR 值Fig.1 EPSR values of parameters in simulation

當(dāng)?shù)螖?shù)為5 000 時，所有參數(shù)的EPSR 值都小于1.2，因此本文進(jìn)行100 次重復(fù)試驗，每次試驗中丟掉前5 000 次迭代值，為避免自相關(guān)，使用后5 000 次間距為2 的迭代值的平均值作為未知參數(shù) θ 在第v次重復(fù)試驗的貝葉斯估計值最終得到各參數(shù)的貝葉斯估計值：均方根誤差：RMS=以及參數(shù)偏差：Bias=匯總得到表2.計算基本生存函數(shù)的均方根誤差：RASE=并在圖2 中畫出RASE中位數(shù)對應(yīng)那次試驗的基本生存函數(shù)擬合圖.

表2 FGM Copula 半?yún)?shù)二元生存模型的參數(shù)估計Tab.2 Parameters estimation of FGM Copula semiparametric bivariate survival model

根據(jù)表2 和圖2，各個感興趣參數(shù)的偏差的絕對值均小于0.1，均方誤差根RMS 小于0.2，估計標(biāo)準(zhǔn)差也均小于0.2，可認(rèn)為參數(shù)估計較為精確和穩(wěn)健，邊際基本生存函數(shù)的擬合曲線與真實曲線很接近，除尾部有輕微偏差外，總體擬合的效果較好.

圖2 數(shù)值模擬中邊際基本生存函數(shù)擬合圖Fig.2 Fitting plot of marginal baseline survival functions in simulation

接下來，基于貝葉斯數(shù)據(jù)刪除影響分析和貝葉斯局部影響分析診斷二元相依生存數(shù)據(jù)是否存在異常點.首先，數(shù)據(jù)產(chǎn)生同前，人為產(chǎn)生異常點或強影響點如下基于(13)式計算K-L 貝葉斯診斷統(tǒng)計量，作圖得到圖3，說明我們建議的數(shù)據(jù)刪除影響分析方法能診斷出潛在的異常點.其次，為驗證本文所提出的局部影響診斷方法，產(chǎn)生新的生存時間 {ti j|i=80,150,200,j=1,2} 作為異常點或強影響點，方法如下：t80j和t150j分別從邊際危險函數(shù)為

圖3 數(shù)值模擬中貝葉斯數(shù)據(jù)刪除影響診斷圖Fig.3 Diagnostic plot of Bayesian data deletion influence analysis in simulation

的生存分布中產(chǎn)生，t200,1和t200,2則從邊際危險函數(shù)是

的生存分布中產(chǎn)生，其余數(shù)據(jù)產(chǎn)生方法同前，然后基于(17)式計算 |hmax|.類似于Zhu 等[22]，選取作一基準(zhǔn)線（m通常取2），表示向量的均值和標(biāo)準(zhǔn)差，結(jié)果如圖4 所示，不難看出所有的潛在影響點都被正確診斷出來，說明本文所提出的統(tǒng)計診斷方法有良好的性能.

圖4 數(shù)值模擬中貝葉斯局部影響診斷圖Fig.4 Diagnostic plot of Bayesian local influence analysis in simulation

3.2 實例分析糖尿病視網(wǎng)膜病變是與糖尿病相關(guān)的并發(fā)癥，糖尿病視網(wǎng)膜病變研究(DRS)始于1971 年，目的是研究激光光凝治療在延緩糖尿病視網(wǎng)膜病變患者失明的有效性[23].每個患者的一只眼睛隨機選擇治療，另一只眼睛不治療，本文使用DRS 標(biāo)準(zhǔn)定義的高?；颊邩颖緸榉治鲎蛹?n=197)對所提出的模型進(jìn)行研究.對于該樣本數(shù)據(jù)，觀測得到的事件時間是患者失明時間或隨訪結(jié)束等原因造成的刪失時間，當(dāng)獲得的是失明時間時，標(biāo)記為1，否則獲得刪失時間，標(biāo)記為0.以患者年齡作為協(xié)變量，生存過程1 用來建模接受治療眼睛的生存數(shù)據(jù)，生存過程2 用來建模未接受治療眼睛的生存數(shù)據(jù)，用本文所提出的半?yún)?shù)二元相依生存模型對數(shù)據(jù)進(jìn)行擬合.設(shè)置樣條次數(shù)為2 并取1 個內(nèi)部節(jié)點，先計算各個感興趣參數(shù)和虛擬參數(shù)的EPSR 值，繪制圖5.顯然當(dāng)?shù)螖?shù)為5 000 時，所有參數(shù)的EPSR 值都小于1.2，因此燃燒前5 000 次，并且為避免自相關(guān)，取后5 000 次間距為2 的迭代樣本的均值作為參數(shù)估計值，得到參數(shù)估計值見表3.

圖5 實例數(shù)據(jù)分析的EPSR 圖Fig.5 EPSR values of parameters in real data analysis

表3 實例數(shù)據(jù)分析的參數(shù)估計Tab.3 Parameters estimation of real data analysis

根據(jù)表3，Copula 相關(guān)參數(shù)φ估計值為0.821 1，且在顯著性水平為0.05 時，95%的可信區(qū)間不包括0，因而我們認(rèn)為接受治療眼睛的生存數(shù)據(jù)和未接受治療眼睛的生存數(shù)據(jù)2 者存在顯著正相關(guān)關(guān)系.另外，β11=?0.018 8 表明在控制其他條件不變的情況下，接受激光光凝治療的患者隨著年齡的上升，失明的時間變短，即治療效果變差，同時得到2 個生存過程的基本生存函數(shù)的擬合圖，見圖6 所示.

圖6 實例數(shù)據(jù)的邊際基本生存函數(shù)模型擬合圖Fig.6 Fitting plot of marginal baseline survival function in real data analysis

為了進(jìn)一步說明所提出模型對實例數(shù)據(jù)的擬合效果，本文用DIC 信息準(zhǔn)則（Deviance Information Criterion）對所提出的二元相依生存半?yún)?shù)模型與二元獨立半?yún)?shù)生存模型（φ=0）、邊際生存分布為Weibull 分布的二元相依全參數(shù)生存模型和二元獨立全參數(shù)生存模型進(jìn)行比較，整理得到表4.

根據(jù)表4，顯然無論是否考慮生存過程的相依性，半?yún)?shù)生存模型的DIC 值均比全參數(shù)模型的DIC 值小，說明半?yún)?shù)模型對數(shù)據(jù)進(jìn)行生存分析要比全參數(shù)模型更好；同時，相依模型的擬合效果也優(yōu)于獨立模型，因而也側(cè)面說明了在對二元生存過程建模時不應(yīng)該忽略生存過程之間的相依性問題.另外，本文所提出的相依性半?yún)?shù)模型的DIC 值為3 370.402 9，比其他3 個模型都要低，說明該模型在這個實例數(shù)據(jù)中具有良好的表現(xiàn)，即本文所提出的模型具有較好的實用意義.

表4 實例數(shù)據(jù)分析模型比較Tab.4 Model Comparison in real data analysis

最后本文對實例數(shù)據(jù)進(jìn)行貝葉斯數(shù)據(jù)刪除影響分析和貝葉斯局部影響分析，得到圖7 和圖8，在數(shù)據(jù)刪除影響分析中第28、102、120、183 個個體被診斷為異常點，在局部影響分析中第28 和49 個個體被診斷為影響點，第23、32、111、120、122 個個體也對模型有一定的擾動影響.

圖7 實例數(shù)據(jù)的貝葉斯數(shù)據(jù)刪除影響診斷圖Fig.7 Diagnostic plot of Bayesian data deletion influence analysis in real data

圖8 實例數(shù)據(jù)的貝葉斯局部影響診斷圖Fig.8 Diagnostic plot of Bayesian local influence analysis in real data

4 總結(jié)與展望

本文主要通過引入FGM Copula 函數(shù)將個體的不同邊際生存過程聯(lián)系起來，基于B 樣條函數(shù)建立二元相依生存分析半?yún)?shù)模型，在貝葉斯理論框架下，對所提出的相依二元生存分析模型進(jìn)行貝葉斯數(shù)據(jù)刪除影響分析和貝葉斯局部影響分析的統(tǒng)計診斷，最后，將所建議模型應(yīng)用到糖尿病視網(wǎng)膜病變研究實例數(shù)據(jù)中，驗證了模型的可行性.

在本文的研究基礎(chǔ)上，下一步研究方向可以有以下方面：首先將相依二元生存分析通過Pair-Copula理論或者多元FGM Copula 進(jìn)一步拓展到多元的情形，但是隨著生存過程的增加，模型的復(fù)雜度也會進(jìn)一步增加，給參數(shù)估計造成困難；其次，本文僅考慮了右刪失的情況，為了更具普遍性，可以拓展到區(qū)間刪失、截尾數(shù)據(jù)，甚至是考慮缺失數(shù)據(jù)的情況.