王寧寧 徐淑一 方積乾

在生物科學中，經(jīng)常遇到同一個體的多次觀測，比如同一個病人手術后疾病經(jīng)常反復發(fā)作，一個文獻中經(jīng)常被引用并研究的數(shù)據(jù)是慢性肉芽腫病(chronic granulomatous disease，CGD)數(shù)據(jù)。Yau和McGILCHRIST (1998)[1]對這種數(shù)據(jù)建立AR(1)結構的Frailty模型，進行了ML估計和REML估計，他們的估計相當于利用調(diào)整的集中化擴展似然函數(shù)分別對方差和相關系數(shù)進行估計，沒有考慮二者的聯(lián)合信息。

由Lee和Nelder(1996)[2]提出的帶隨機效應的廣義線性模型的等級似然(hierarchicacl likelihood)估計方法，允許隨機效應因子是任何分布，等級似然方法近年來開始應用于生存數(shù)據(jù)的混合模型估計，展現(xiàn)了良好的應用前景。LI.DO.HA等(2001)[3]將等級似然估計方法用于脆弱模型的估計，估計了隨機效應(脆弱性因子)呈正態(tài)分布和γ分布的情形，并證明了在給定隨機效應分布參數(shù)的情況下，最大化邊際似然得到的協(xié)變量系數(shù)估計和最大化等級似然得到的協(xié)變量系數(shù)估計是一致的。LI.DO.HA等(2002)[4]將等級似然估計應用于存在刪失數(shù)據(jù)的混合線性模型，假設隨機效應呈正態(tài)分布。LI.DO.HA等(2005)[5]將其應用于多重混合線性模型的估計。王寧寧等(2009)[6]將其應用于脆弱性模型的檢驗。徐淑一和王寧寧(2007，2009，2011)[7-10]將其應用于競爭風險模型的估計與應用研究中。王寧寧等(2011)[11]使用這種方法對威布爾形式的脆弱性模型進行估計和應用。近年來，文獻上關于等級似然估計方法在生存分析數(shù)據(jù)中的研究逐漸增多，基于篇幅，不再一一列舉。

本文針對個體重復觀測的生存數(shù)據(jù)，采用AR(1)形式的frailty結構，建立比例危險模型。本文采用等級似然估計方法估計協(xié)變量參數(shù)并預測隨機效應的實現(xiàn)值，采用MAPHL方法估計隨機效應分布的參數(shù)，并且將ML、REML方法和MAPHL方法比較，并證實了REML估計等價于調(diào)整的輪廓等級似然關于θ和ψ的一階導數(shù)等于零的迭代公式。

模型設計和參數(shù)估計過程

假設在給定的可觀測的協(xié)變量X以及不可觀測的異質(zhì)性因子U的條件下，建立比例危險模型如下：λ(t|u)=λ0(t)exp(Xβ)u,記v=log(u)，考慮有q個病人(q個組)，第i個病人重復觀測次數(shù)為ni，則第i個病人的第j次生存時間的危險率模型為：

λ(tij|u)=λ0(tij)exp(ηij)

(1)

(2)

其中，

(3)

在給定(X,V)條件下，Tij的條件生存函數(shù)是：

S(t|X,V)=exp(-Λ(s)exp(ηij))

(4)

那么Tij的條件密度是：

f(t|X,V)=h(s|X,V)exp(-Λ(s)exp(ηij))

(5)

于是，等級似然函數(shù)可以寫成：

(6)

l1ij=δij{logλ0(yij)+ηij}-{Λ0(yij)exp(ηij)}

(7)

(8)

(9)

(10)

下面詳細討論模型參數(shù)的估計步驟。我們設第i組的隨機效應的方差陣為：

(11)

(12)

重復這一步直到收斂。

(13)

重復這一步直到收斂。其中，

(14)

隨機效應分布參數(shù)的MAPHL估計與REML以及ML估計的比較

Yau和McGILCHRIST (1998)推導了隨機效應分布參數(shù)的REML估計和ML估計。我們采用Lee 和Nelder(1996)的MAPHL估計方法，和REML及ML方法比較，經(jīng)過推導，得到以下結論：

證明：

(15)

證明：

(16)

(17)

上述REML或者ML估計都相當于調(diào)整的等級似然對θ與φ單獨求導令其為零得到，這么做顯然沒有考慮θ與φ的聯(lián)合信息，因此我們認為聯(lián)合估計θ與φ更加合理。本文采用SAS的IML語言分別編寫了MAPHAL估計和REML以及ML估計的程序，完整的實現(xiàn)了前述估計過程。

應用實例

這一部分對本文開始提到的CGD數(shù)據(jù)做應用研究，總共128個病人，來自14個不同的醫(yī)療機構，主要考察γ干擾素對慢性肉芽腫病的治療效果，一部分病人采取了γ干擾素措施，其余的病人采用安慰劑作對比。變量rIFN-g=0表示安慰劑，rIFN-g=1表示采用該療法。為了充分考慮患者個體的不同狀況以及盡量控制其他因素的影響，模型還包括了病人身體特征的協(xié)變量：病人進入研究時的年齡(Age)、性別(Sex)、身高(Height)以及體重(Weight)；病人的遺傳模式(Pattern of Inheritance)，用Inheritance表示；病人在進入研究時是否使用皮質(zhì)類固醇(Corticosteroids)；病人在進入研究時是否采用了預防性的抗生素(prophylactic antibiotics)。除此之外，模型還考慮了不同的醫(yī)療中心的因素，將14個醫(yī)療中心分為兩類，一類來自于美國，一類來自于歐洲的醫(yī)療中心，用Institution category表示。由于很多病人都有重復觀測，多的達到7次，需要考慮個體的異質(zhì)性影響，我們估計九個協(xié)變量的AR(1)結構的隨機效應模型，結果列在表1中(括號內(nèi)表示估計的標準差)。

表1 CGD數(shù)據(jù)的AR(1)隨機效應模型的估計結果

我們發(fā)現(xiàn)各種估計方法的結果是接近的，而且，協(xié)變量參數(shù)REML估計和MAPHL估計結果很接近，總的說來，三種方法差別不大，各個變量的顯著性很接近。隨機效應因子方差的REML和MAPHL估計結果很接近,然而φ的ML和MAPHL估計比較接近，而且非常小，接近于零。不同病人的異質(zhì)性影響我們認為是存在的，但是φ的估計結果卻值得懷疑。我們的估計結果和Yau和McGILCHRIST (1998)的回歸系數(shù)估計結果比較接近，但是θ和φ的估計差別比較大，原因還有待于進一步研究。

結論和進一步的研究方向

本文針對包含右刪失的重復觀測的縱列生存數(shù)據(jù)，考慮到同一個體的重復觀測的生存時間問題，在Cox比例危險模型的基礎上，建立了AR(1)結構的隨機效應模型，推導了等級似然函數(shù)，并且和REML、ML方法進行了比較，證明了θ和φ的REML的迭代公式正是調(diào)整的集中化等級輪廓似然對θ和φ的一階導數(shù)為零形成的迭代公式。本文進一步的研究工作是將本文提出的針對AR(1)隨機效應的比例危險模型進行模擬研究，并探討檢驗隨機效應是否存在AR(1)結構的合適的方法。

參 考 文 獻

1.Yau K,McGilchrist C.ML and REML estimation in survival analysis with time dependent correlated frailty.Statistics in Medicine，1998，17(11):1201-1213.

2.Lee Y,Nelder JA.Hierarchical generalized linear models.Journal of the Royal Statistical Society.Series B (Methodological)，1996:619-678.

3.Ha I D,Lee Y,Song JK.Hierarchical likelihood approach for frailty models.Biometrika，2001，88(1):233-233.

4.Ha I D,Lee Y,Song JK.Hierarchical-likelihood approach for mixed linear models with censored data.Life time data analysis，2002，8(2):163-176.

5.Ha I D,Lee Y.Multilevel mixed linear models for survival data.Lifetime Data Analysis，2005，11(1):131-142.

6.王寧寧,徐淑一,方積乾.脆弱性模型的隨機效應檢驗和推廣的擬合檢驗.中國衛(wèi)生統(tǒng)計，2009，26(1):2-6,10.

7.徐淑一,王寧寧.競爭風險下縱列數(shù)據(jù)的隨機效應建模和估計.中山大學學報(自然科學版)，2007，46(1):7-10.

8.徐淑一,王寧寧.競爭風險下我國住房抵押貸款風險的實證研究.統(tǒng)計研究，2011，28(2):45-52.

9.徐淑一,王寧寧,王美今.競爭風險下縱列持續(xù)數(shù)據(jù)隨機效應模型的估計與模擬研究.數(shù)理統(tǒng)計與管理，2009，28(6):1013-1023.

10.徐淑一,王寧寧,王美今.基于持續(xù)期數(shù)據(jù)的我國住房抵押貸款違約和提前還款風險分析.南方經(jīng)濟，2009(6):34-43.

11.Wang N,Xu S,Fang J.Hierarchical likelihood approach for the Weibull frailty model.Journal of Statistical Computation and Simulation，2011，81(3):343-356.

12.Breslow N.Covariance analysis of censored survival data.Biometrics，1974:89-99.