威布爾脆弱性比例危險模型的等級似然估計和應用*

王寧寧，徐淑一

(1.廣州醫科大學公共衛生學院，廣東 廣州 510182；2.中山大學嶺南學院，廣東 廣州 510275)

Cox[1]比例危險模型自從1972年提出以來一直是占據生存數據分析的主要模型，在面對同一個體的重復觀測時(或者數據按照不同來源分組時)，對生存數據的建模需要考慮個體的異質性影響。Keyfitz和Littman發現忽略個體的異質性會導致對協變量系數的估計有偏。Hougarrd與Aalen也發現，忽略個體的異質性時，估計出的相對危險率都有偏高的傾向，Lancaster對失業持續時間建模時，發現模型中忽略脆弱性因子會導致低估協變量效應。Pickles等[2]對脆弱性模型的回顧性總結中表明，一般來說，模型中忽略異質性的影響會使得協變量系數的估計向0這個方向接近。將脆弱性(Frailty)因子引入模型就是考慮到這種不可忽略的異質性。脆弱性模型是Cox提出的比例危險模型的一個推廣，脆弱性模型允許同一個體的重復觀測數據之間或者同一個組內的個體之間具有相關性，廣泛應用于多元生存數據分析中。

經典的脆弱性比例危險模型h(tij|vi)=h0(tij)exp{xijβ+vi}=h0(tij)exp{xijβ}ui，下標ij表示第i個體的第j次重復觀測。除了參數向量β之外，v的分布函數F(v)以及基本危險率h0都是未知的。在很多脆弱性比例危險模型的應用中，基本危險率通常認為是不規則的形式，通過非參數方法估計。有很多學者對基本危險率是非參數設定的脆弱性模型的估計進行了研究，McGilchrist 和Aisbett以及McGilchrist研究了對數正態脆弱性模型，利用Cox的部分似然(Partial Likelihood)方法估計參數，并且使用這種模型來分析腎感染的數據，但是沒有處理數據打結的情況。近年來的研究主要集中于非正態脆弱性模型的估計，比如，出于數學上處理的方便，Klein及Nielsen等建立了伽馬(Gamma)脆弱性模型的EM算法。Hougaard建議了三種脆弱性因子的分布：伽馬分布、逆高斯(Inverse Gaussian)分布和正穩定(Positive stable)分布，Hougaard[3]詳細討論了脆弱性因子的分布的選擇。EM算法在估計帶有不可觀測的協變量模型是很有用的技術，見Mclachlan等[4]。EM算法的收斂對初始值的選擇和停止計算的規則非常敏感，而且EM算法需要脆弱性因子在給定觀測數據下的條件期望，除了伽馬分布，逆高斯分布等幾種特定的分布，一般情況下，無法寫出解析表達式，而是需要數值算法，而且參數估計的方差也不能直接獲得(Louis，Jamshidian等[5])。

由Lee等[6]提出的帶隨機效應的廣義線性模型的等級似然(Hierarchicacl Likelihood)估計方法，允許隨機效應因子是任何分布，等級似然方法近年來開始應用于生存數據的混合模型估計，展現了良好的應用前景。Ha等[7]將等級似然估計方法用于脆弱模型的估計，估計了隨機效應(脆弱性因子)呈正態分布和Gamma分布的情形，并證明了在給定隨機效應分布參數的情況下，最大化邊際似然得到的協變量系數估計和最大化等級似然得到的協變量系數估計是一致的。自此，等級似然估計方法在縱列數據的生存分析與混合模型方面的應用研究進一步展開，如Ha等[8-9]、徐淑一等[10]的研究。近年來學術界關于等級似然方法的應用研究仍在繼續，如Lee等[11]使用等級似然方法估計預測疾病測繪時相對風險的的區間，模擬顯示其不弱于貝葉斯方法。Noh等[12]使用等級似然建立非線性混合效應模型處理縱向數據。Noh等[13]使用等級似然方法提供了一種減弱對數據缺失機制不正確假定的影響并提供了實例和模擬研究。Lee等[14]把等級似然函數方法用于缺失數據的建模并應用于厚尾分布的縱向數據分析。Wu等[15]使用等級似然解決因子分析模型中二元相應變量的情形，使用簡單而且高效。王寧寧等將該方法推廣到AR(1)相依結構的脆弱性模型研究。

本文考慮Cox比例危險的脆弱性模型，允許數據打結以及刪失的情況下，假設脆弱性因子服從威布爾分布。對固定參數的估計和隨機效應實現的預測采用聯合最大化擴展似然方法，實際上擴展似然即為等級似然，但是等級似然要求特定的隨機效應尺度，見Lee和Nelder等；而隨機效應分布參數則采用調整的輪廓等級似然進行估計。本文推導了模型的等級似然函數和估計過程，并將不同分布的脆弱性模型做對比研究。

1 隨機效應是威布爾分布的半參數比例危險模型的估計

威布爾分布廣泛應用于生存時間以及可靠性試驗的統計分析中，Peto和Lee認為在一定條件下失效時間應該是威布爾分布。威布爾分布的危險率函數是h(t)=γtγ-1(尺度參數設為1)，當形狀參數γ大于1時，危險率隨時間而增加；當γ等于1時，危險率函數是常數；當γ小于1時，危險率是減函數。威布爾分布靈活的危險率形式使得它非常適合作為生存時間的模型，在生命科學領域，已經有許多的統計學家采用威布爾分布來分析數據。由于威布爾分布的危險率函數的簡單性和靈活性，威布爾分布在多元生存分析中也非常有用。Kimber和Crowder在對異質性建模時使用威布爾分布作為生存分布，Lancaster在考慮異質性的持續數據的比例危險模型中采用威布爾分布。

令Tij(i=1,2,…,q，j=1,2,…,ni)是第i個個體(組)的第j次重復觀測，假定給定Ui=ui的的條件下，數據對(Tij,Cij)是條件獨立的，而且Tij與Cij條件獨立；并且給定Ui=ui，{Cij,j=1,2,…,ni}不含有ui的信息。Cij為刪失時間，設Tij的條件危險率為比例危險形式，Ui之間獨立同分布，設其服從威布爾分布，考慮危險率函數

hi(t|ui)=h0(t)exp{xiβ}uiexp(-E(logui))

(1)

ui服從威布爾分布，其密度函數為

f(ui)=λρ(λui)ρ-1exp{-(λui)ρ}

(2)

對脆弱性因子ui進行變換，使ui的對數為零均值，令ui=ewi，那么wi=log(ui)具有極值分布

(3)

此時，

(4)

hi(t|ui)=h0(t)exp{xiβ+vi}

(5)

于是，等級似然函數可以寫成

(6)

其中

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

(7)

(8)

(9)

對h*關于β,v一起最大化，可以得到β和v的估計，估計過程采用Newton-Raphson迭代。

對隨機效應分布參數b的估計，Lee和Nelder建議使用最大化調整的輪廓似然估計，調整的輪廓似然為

(10)

2 各種隨機效應模型的比較和模擬研究

為了進一步研究不同的脆弱性模型，下面對隨機效應的不同分布生成模擬數據，采用不同的脆弱性模型進行估計。為了進行比較，也把帶有隨機效應生成的生存數據用一般的Cox 比例危險模型進行估計并和隨機效應模型估計結果進行比較，用以檢查忽略隨機效應對模型估計的影響，模型設為

(11)

其中，i=1,2,…,q，j=1,2,…,ni，也就是說假定q個病人，每個病人有ni個重復觀察；或者觀察數據按照其自然狀態有q個組，每個組內有ni個個體。

表1 正態隨機效應生成的生存數據的估計結果1)

1)Cox是指Cox比例危險模型，LNM指對數正態脆弱性模型，GAM指伽馬脆弱性模型，WEL指威布爾脆弱性模型；MEAN指200次模擬估計的平均值，SD是200次模擬估計的標準離差，SE是指平均的標準差，95％CP為95％的置信區間包含參數真值的比例。下面表2、表3相同

由表1可以發現，刪失率為0.1時，三種隨機效應模型的估計結果差別不大，尤其是伽馬脆弱性模型的估計最為精確,絕對誤差為0.001 4，似乎最為合適的脆弱性模型應該是伽馬脆弱性模型。刪失率的增加到0.3時，不再是伽馬脆弱性模型估計最精確，而是威布爾脆弱性模型估計最精確，絕對誤差為0.012 4，而此時伽馬脆弱性模型最差。刪失率增加到0.5時，這種變化趨勢已經很明顯了，此時威布爾脆弱性估計最精確，絕對誤差僅為0.011 6，而且和相對低的刪失率情況下估計結果差別很小，而其它兩種脆弱性模型誤差都比較大，尤其伽馬脆弱性精確度很差。從這三種脆弱性模型的估計結果來看，假如真實的隨機效應是對數正態的，估計時分布的誤設帶來的誤差不是很明顯，然而，當刪失率提高時，威布爾脆弱性模型卻能大大改善回歸系數的估計。進一步地，我們發現，模擬估計的標準差和平均的SE基本一致，說明用等級似然函數關于參數的二階導數矩陣地相反數地逆作為估計的方差也很精確。

表2 伽馬隨機效應生成的生存數據的估計結果

由表2的結果可以看出，刪失率為0.1時，三種隨機效應模型估計結果差不多，絕對誤差都比較小，伽馬模型估計的最精確，但是刪失率增加到0.3，再增加到0.5的時候，伽馬脆弱性模型估計結果不如另外兩種模型精確，這種變化趨勢很明顯，尤其刪失率比較高時，比如說0.5時，威布爾脆弱性模型估計最為精確，對數正態模型次之。

表3 威布爾隨機效應生成的生存數據的估計結果

上述三種分布的隨機效應生存數據的模擬估計還顯示：如果不采用隨機效應模型，協變量系數的估計將會嚴重偏離，整體上均呈現高估的趨勢，絕對值上是低估，這說明，如果模型中忽略異質性(Frailty)的影響，會使得估計結果存在很大謬誤。模擬還發現，即使采用了錯誤的脆弱性模型，也比不用脆弱性模型估計結果好得多，這說明，實際中，對縱列生存數據建模和估計，要非常小心數據是否存在個體或者組別的異質性，以及注意模型是否忽略了重要因素而導致估計不準確。在實際中，我們推薦使用威布爾脆弱性模型。

模擬研究還發現，用不正確的脆弱性模型以及隨著刪失率的增加，估計也呈現高估趨勢(協變量參數在絕對之上是低估趨勢)。究其原因，我們認為是由于隨機效應的存在，一方面增加了參數估計的不確定性，另一方面，由于隨機效應問題，使得生存數據的波動存在于兩個方面，一是協變量的影響，二是隨機效應的影響，如果隨機效應高估，那么協變量參數會呈現低估現象。觀察表1到表3的模擬結果，我們發現，當隨機效應方差低估時，協變量系數會高估；當隨機效應方差高估時，協變量系數則低估。

威布爾隨機效應生成的生存數據，用威布爾脆弱性模型估計的隨機效應分布參數也相當精確，δ的真實值為1，刪失率為0.1、0.3、0.5時估計的參數分別為0.978 69(0.184 45)、0.962 53(0.170 6)、0.915 83(0.154 32)，括號內為模擬估計的標準差，可以看出，他們也是顯著的，而且，對威布爾隨機效應模擬的數據，用調整的輪廓等級似然函數的二階導數估計方差，也是比較精確的，因為估計的標準差和用二階導數計算的平均的SE差別不大，僅有一點高估。其它兩種脆弱性模型的隨機效應分布參數的估計也很精確，但是估計量的方差用調整的輪廓似然的二階導數計算則偏誤較大。由于不同隨機效應分布的參數涵義不一致，不同脆弱性模型隨機效應參數的直接比較沒有意義，但是按照和數據生成過程一致的脆弱性模型估計的隨機效應分布參數與真值非常接近。

模擬還顯示，對帶有隨機效應的生存數據，不能忽略隨機效應的存在，否則，Cox模型會有很大偏誤。對不同隨機效應的生存數據，當然理論上是能夠使用真實分布的脆弱性模型估計會很好，但是模擬研究發現，對正態分布、伽馬分布、威布爾分

布這三種隨機效應的生存數據來說，威布爾脆弱性模型估計結果顯示出很大的優勢，尤其當刪失率提高時，這種優勢最為明顯。那么，我們可以得出這樣的結論，當數據存在較大比例的刪失時，不妨采用威布爾脆弱性模型。至于為何威布爾脆弱性模型顯示出這么明顯的優勢，還有待于進一步研究。本文本文模擬采用的樣本容量為200，而實際上樣本容量為50的情況下，本文的模型也給出了非常類似的結果，這表明等級似然估計方法應用于脆弱性模型估計的時候對于小樣本情形有很好的近似。

3 應用實例

下面對文獻上常見的兩個縱列生存數據的例子用前面討論的各種脆弱性模型進行估計，一個是腎病感染數據，另一個是燒傷病人皮膚移植數據。

腎病感染數據來自于McGilchrist和Aisbett，總共有38個腎病患者做便攜式腎透析，記錄下在導管插入處發生感染的時間，從導管插入開始時記錄時間，每個患者均被記錄了兩次發生感染的時間，如果某次感染發生時，導管已經由于感染或其它原因被移除，則這次感染時間被記錄為刪失。研究中主要考慮五個協變量的影響：年齡，用x1表示；性別，用x2表示，x2=1為女性，x2=0為男性；是否為血球性腎炎用x3表示，x3=1表示是血球性腎炎，否則x3=0；是否為急性腎炎，用x4表示，x4=1表示是急性腎炎，否則x4=0；是否為多囊腎病，用x5表示，x5=1表示是多囊腎病，否則x5=0否則?？紤]每個病人的異質性，建立脆弱性模型為

H(tij|x1,x2,x3,x4,x5;ui)=h0(tij)·

exp(β1x1ij+β2x2ij+β3x3ij+β4x5ij+β5x5ij+vi)，

i=1,2,…,38，j=1,2

每個病人有一個異質性因子vi，以此來體現每個病人兩次感染之間的聯系。對腎感染數據，采用一般Cox比例危險模型、對數正態脆弱性模型、伽馬脆弱性模型、以及威布爾脆弱性模型進行估計，結果列在表4中，可以看出，性別比較顯著，是個危險因素，這說明男性比女性更容易發生感染，這個結果與McGilchrist和Aisbett使用的邊際似然估計的結果一致。不同脆弱性模型估計結果有一定的差別，而且都顯示存在顯著的隨機效應，也就是說不同病人除了性別以及疾病類型之外，還有明顯的個體差異。

表4 腎感染數據的各種模型估計結果

但是Klein的伽馬邊際似然方法給出的結果和一般Cox模型結果類似，而且利用Klein的檢驗說明該例不存在個體異質性。通過上一章的模擬研究，我們發現用調整的等級輪廓似然的二階導數估計隨機效應分布參數估計值的方差會偏低，但偏低并不是很嚴重，本例的估計結果顯示存在個體異質性。Klein的異質性檢驗方法是個很一般的檢驗，功效往往很弱，經常檢驗不出異質性。Yus T.BOENG在他的博士論文研究中對參數基本危險率的威布爾脆弱性模型估計結果也證實存在個體異質性，而且提出了針對威布爾分布的基本危險率的得分檢驗證實存在異質性。根據上一節的模擬研究結論，當數據存在適量刪失時，建議采用威布爾Frailty模型。

皮膚移植數據來自Batchelor 和hackett，研究了16個嚴重燒傷的病人，由于每一個病人多塊皮膚需要移植，研究每塊皮膚的治療效果，有的病人只有一塊皮膚移植，而有的病人有多塊，這是一個非平衡的觀測數據，不同病人做皮膚移植的塊數不同，最多的有四塊，因此同一個體可以有多個生存時間，適合建立脆弱性模型研究。這個數據只考慮了一個協變量因素，就是人白細胞抗原(human leukocyte antigen，HLA)的匹配好壞，模型設為：H(tij|x;ui)=h0(tij)exp(β×HLAij+vi)，HLA為1表示匹配的好，0表示匹配的不好，i=1,2,…,16，j=1,…,ni。采用一般Cox比例危險模型和三種脆弱性模型進行估計，結果列在表5中。

表5 皮膚移植數據的各種模型估計結果

我們發現，HLA的影響是顯著的，是影響皮膚移植生存的重要因素，不同脆弱性模型估計結果差別很小，對數正態脆弱性模型和威布爾脆弱性模型估計結果比較接近，伽馬脆弱性模型和Cox模型估計結果比較接近，而且都顯示了存在顯著的隨機效應，尤其是威布爾模型，顯示個體存在明顯的異質性。

4 結 論

本文主要考察了并非指數族分布的威布爾脆弱性模型的建模，使用等級似然方法對其進行估計。結果顯示，等級似然方法可以很好給出威布爾脆弱性模型的估計結果，調整的輪廓等級似然可以用來估計得到精確的隨機效應的參數。模擬結果顯示在縱列持續數據分析和應用中，對帶有隨機效應的生存數據，不能忽略隨機效應的存在，否則，Cox模型會有很大偏誤。對不同隨機效應的生存數據，當然理論上是能夠使用真實分布的脆弱性模型估計會很好，但是模擬研究發現，對正態分布、伽馬分布、威布爾分布這三種隨機效應的生存數據來說，威布爾脆弱性模型估計結果顯示出很大的優勢，尤其當刪失率提高時，這種優勢最為明顯。本文模擬和應用研究均表明，等級似然估計在實際應用中對于中等以下規模(樣本容量50以下)的樣本數據具有很好的穩健性。因此，除了伽瑪脆弱性模型和對數正態脆弱性模型之外，威布爾脆弱性模型是也是一個非常值得應用的模型。

參考文獻：

[1]COX D R.Regression models in life tables (with discussion) [J]．J Roy Statist Soc：Ser B，1972,34:187-220.

[2]PICKLES R,CROUCHLEY.A comparison of frailty models for multivariate survival data [J].Statistics in Medicine,1995,14(13):1447-1461.

[3]HOUGAARD P.Frailty models for survival data [J].Lifetime data analysis,1995,1(3):255-273.

[4]MCLACHLAN G,PEEL D.Mixtures of factor analyzers [J].Finite Mixture Models,2000:238-256.

[5]JAMSHIDIAN M,JENNRICH R I.Acceleration of the EM algorithm by using quasi-Newton methods [J].Journal of the Royal Statistical Society:Series B (Statistical Methodology),1997,59(3):569-587.

[6]LEE Y,NELDER J A.Hierarchical generalized linear models [J].Journal of the Royal Statistical Society：Series B (Methodological),1996:619-678.

[7]HA IL DO,LEE Y,SONG J K.Hierarchical likelihood approach for frailty models [J].Biometrika,2001,88(1):233-233.

[8]HA IL DO,LEE Y,SONG J K.Hierarchical-likelihood approach for mixed linear models with censored data [J].Lifetime data analysis,2002,8(2):163-176.

[9]HA IL DO,LEE Y.Multilevel mixed linear models for survival data [J].Lifetime Data Analysis,2005,11(1):131-142.

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

[11]LEE Y,JANG M,LEE W.Prediction interval for disease mapping using hierarchical likelihood [J].Computational Statistics,2011,26(1):159-179.

[12]NOH M,LEE Y,KENWARD M G.Robust estimation of dropout models using hierarchical likelihood [J].Journal of Statistical Computation and Simulation,2011,81(6); 693-706.

[13]NOH M,WU L,LEE Y.Hierarchical likelihood methods for nonlinear and generalized linear mixed models with missing data and measurement errors in covariates [J].Journal of Multivariate Analysis,2012,109:42-51.

[14]LEE D,LEE Y,PAIK M C,et al.Robust inference using hierarchical likelihood approach for heavy-tailed longitudinal outcomes with missing data:An alternative to inverse probability weighted generalized estimating equations [J].Computational Statistics & Data Analysis,2013,59:171-179.

[15]WU J,BENTLER P M.Application of H-likelihood to factor analysis models with binary response data [J].Journal of Multivariate Analysis,2012,106:72-79.