王寧寧 徐淑一 方積乾

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

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

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

模型設(shè)計(jì)和參數(shù)估計(jì)過(guò)程

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

λ(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)

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

(6)

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

(7)

(8)

(9)

(10)

下面詳細(xì)討論模型參數(shù)的估計(jì)步驟。我們?cè)O(shè)第i組的隨機(jī)效應(yīng)的方差陣為：

(11)

(12)

重復(fù)這一步直到收斂。

(13)

重復(fù)這一步直到收斂。其中，

(14)

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

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

證明：

(15)

證明：

(16)

(17)

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

應(yīng)用實(shí)例

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

表1 CGD數(shù)據(jù)的AR(1)隨機(jī)效應(yīng)模型的估計(jì)結(jié)果

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

結(jié)論和進(jìn)一步的研究方向

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

參 考 文 獻(xiàn)

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.王寧寧,徐淑一,方積乾.脆弱性模型的隨機(jī)效應(yīng)檢驗(yàn)和推廣的擬合檢驗(yàn).中國(guó)衛(wèi)生統(tǒng)計(jì)，2009，26(1):2-6,10.

7.徐淑一,王寧寧.競(jìng)爭(zhēng)風(fēng)險(xiǎn)下縱列數(shù)據(jù)的隨機(jī)效應(yīng)建模和估計(jì).中山大學(xué)學(xué)報(bào)(自然科學(xué)版)，2007，46(1):7-10.

8.徐淑一,王寧寧.競(jìng)爭(zhēng)風(fēng)險(xiǎn)下我國(guó)住房抵押貸款風(fēng)險(xiǎn)的實(shí)證研究.統(tǒng)計(jì)研究，2011，28(2):45-52.

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

10.徐淑一,王寧寧,王美今.基于持續(xù)期數(shù)據(jù)的我國(guó)住房抵押貸款違約和提前還款風(fēng)險(xiǎn)分析.南方經(jīng)濟(jì)，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.