條件脆弱模型在評價癲癇復發(fā)事件數(shù)據(jù)中的應用及實現(xiàn)*

羅天娥 趙晉芳 薩 建 劉成芳 劉桂芬△

條件脆弱模型在評價癲癇復發(fā)事件數(shù)據(jù)中的應用及實現(xiàn)*

羅天娥1趙晉芳1薩 建1劉成芳2劉桂芬1△

目的 探討條件脆弱模型在癲癇復發(fā)事件數(shù)據(jù)中的應用及軟件實現(xiàn)。方法 利用癲癇復發(fā)數(shù)據(jù)構(gòu)建條件脆弱模型，采用懲罰偏似然函數(shù)進行參數(shù)估計、并與共享脆弱模型和Cox比例風險回歸模型進行對比分析。結(jié)果 條件脆弱模型分析復發(fā)數(shù)據(jù)既考慮復發(fā)時間的非獨立性也考慮了異質(zhì)性問題，可以用來評價癲癇復發(fā)臨床療效，結(jié)果解釋合理，軟件容易實現(xiàn)。結(jié)論 條件脆弱模型可以較好分析復發(fā)數(shù)據(jù)，值得推廣應用。

條件脆弱模型 懲罰偏似然 癲癇復發(fā) 復發(fā)事件數(shù)據(jù)

*:國家青年科學基金項目資助(編號81001294);山西醫(yī)科大學青年基金資助(02200913)

1．山西醫(yī)科大學衛(wèi)生統(tǒng)計教研室(030001)

2．山西醫(yī)科大學解剖教研室

△通信作者:劉桂芬，E-mail:liugf66@126．com

生物醫(yī)學研究中常常對反復發(fā)作的醫(yī)學事件及其影響因素感興趣，如反復感染、反復住院、腫瘤的復發(fā)等，這類同一事件在同一受試對象上反復發(fā)生的數(shù)據(jù)，稱為復發(fā)事件數(shù)據(jù)(recurrent event data)〔1－3〕，該資料不僅記錄事件是否反復發(fā)作，還需記錄反復發(fā)作的時間間隔。復發(fā)事件資料的主要特點是同一結(jié)局反復發(fā)生、復發(fā)時間有刪失情況、具有非獨立性和異質(zhì)性的特點。因此，常用的針對生存資料的Cox比例風險回歸模型就不再適用。目前的復發(fā)數(shù)據(jù)的統(tǒng)計分析方法主要有脆弱模型和方差校正模型。本研究擬采用條件脆弱模型來評價癲癇復發(fā)事件數(shù)據(jù)的臨床療效及R軟件實現(xiàn)，并與共享脆弱模型和Cox比例風險回歸模型進行對比分析。

原理和方法

1.條件脆弱模型

脆弱模型(frailty models)〔4〕是包含隨機效應的比例風險模型，脆弱性即模型中的隨機效應，模型利用脆弱因子來描述生存時間的相關(guān)性。共享脆弱模型(shared frailty models)〔5－6〕假定同組個體或同一個體在不同時點具有相同的脆弱，即“共享”脆弱，在此條件下滿足比例危險假定。

條件脆弱模型(conditional frailty model)〔7－8〕是基于事件進行分層，同時引入隨機效應即脆弱項，來分析復發(fā)事件數(shù)據(jù)的非獨立性和異質(zhì)性;同一個體復發(fā)事件的相依性可用事件分層，即通過改變第k次事件的基線風險來反映，數(shù)據(jù)模式是把總時間按事件的復發(fā)次數(shù)分為多個時間段，風險集的定義是第k－1次事件發(fā)生后，k次事件才有發(fā)生的風險，參數(shù)估計解釋為即k－1次事件開始發(fā)生到第k次事件發(fā)生的風險估計;個體間異質(zhì)性用隨機效應來反映，針對右截尾失效時間個體i發(fā)生第k次事件的風險，用hik表示，建立條件脆弱模型如下:

式中，k表示第k次復發(fā)事件數(shù)，h0k是基線風險函數(shù)，可以隨k不同而改變，解釋事件依賴性，如果h0k=h0，表示事件間無依賴性，如果h0k=kh0，表示事件間有依賴性存在;t－tk－1表示從前一次即k－1次事件開始發(fā)生到第k次事件發(fā)生的間隔時間，X表示自變量，可以是分類變量也可以是隨時間變化的時依協(xié)變量，β是各協(xié)變量對應的回歸系數(shù)，ui=exp(ωi)表示個體i的隨機效應或脆弱項，即個體i在所有的復發(fā)時間里共享一個脆弱項，解釋異質(zhì)性的存在，如果ω=0，表示沒有未觀測的異質(zhì)性存在，如果ω～N(0，1)則有異質(zhì)性存在，θ是u的方差，它越大表示異質(zhì)性程度越大。

2.參數(shù)估計

懲罰偏似然(penalized partial likehood，PPL)方法〔9－12〕與EM算法相似，脆弱項被認為是受懲罰函數(shù)約束，避免不同組脆弱間差異太大，計算上與其他懲罰回歸的收縮估計(shrinkage estimator)方法相似，如嶺回歸(ridge regression)和光滑樣條(smoothing splines)，該方法收斂速度較快。

懲罰偏似然函數(shù)是Cox回歸偏似然函數(shù)的修飾函數(shù)，回歸系數(shù)和脆弱項被包括到風險函數(shù)里進行優(yōu)化。該模型的懲罰偏似然函數(shù)為:

第一項是條件脆弱模型的偏似然函數(shù):

式中，k表示復發(fā)事件數(shù)，δ是截尾變量(δ=1表示完全數(shù)據(jù)，δ=0表示截尾數(shù)據(jù))，Y是風險指示變量(Y=1表示個體在第k次事件處于風險狀態(tài)，Y=0表示個體未處于風險狀態(tài))。

第二項是Gamma脆弱的懲罰項:

f是懲罰函數(shù)用來限制ω，θ是Gamma脆弱分布的方差，組內(nèi)相關(guān)性用τ=θ/(2+θ)表示。參數(shù)估計時，首先設定脆弱值為1進行模型擬合，然后將脆弱項作為固定或已知的參數(shù)，第一步通過迭代對偏似然函數(shù)進行優(yōu)化;第二步，利用現(xiàn)有觀測值將脆弱項作為條件均數(shù)進行評估，如此重復直到收斂。該法的優(yōu)點是計算簡單，可以快速達到優(yōu)化。

3.基于Schoenfeld殘差比例風險假定的作圖法和檢驗法

利用Schoenfeld殘差圖檢驗是否滿足比例風險假定，Schoenfeld(1982)〔13〕定義偏殘差ri=(ri1，ri2，…，rik)'，rik=Xik－E(Xk|Ri)，即第k個變量ti時刻的偏殘差Xik的觀察值與給定危險集Ri下條件期望之差;Schoenfeld殘差不是每個個體有一個獨立的殘差，而是每個個體在每個協(xié)變量下有一個單獨的殘差。

實例分析

在一項臨床試驗〔12〕研究中，將254例癲癇病人按隨機雙盲法分為陽性藥物組和安慰劑對照組，從癲癇病人進入試驗開始記錄其第一次癲癇發(fā)作和第二次發(fā)作的間隔時間，隨訪一年后中止試驗。復發(fā)事件數(shù)據(jù)以計數(shù)過程的輸入方式記錄，例如一個個體有三次復發(fā)和一次截尾共有四次觀測，每次觀測有一個起始時間和一個終止時間，需要一個事件指示變量說明終止時間是復發(fā)時間或是截尾時間。數(shù)據(jù)格式見表1及變量賦值見表2。

表1 254例癲癇病人反復發(fā)作部分數(shù)據(jù)

表2 變量賦值表

用R 2.15.1版免費軟件來擬合Cox比例危險模型、共享Gamma脆弱模型和條件共享Gamma脆弱模型，程序編輯如下:

表3 癲癇復發(fā)數(shù)據(jù)不同模型結(jié)果比較

結(jié)果整理見表3，可見共享伽瑪脆弱模型和條件脆弱模型的標準誤均大于Cox比例危險回歸模型，條件脆弱模型的標準誤最大，這表明由于Cox比例危險回歸模型沒有考慮復發(fā)事件相關(guān)性的存在，導致方差低估，可信區(qū)間變窄，得出不滿足比例風險假定的結(jié)論;共享伽瑪脆弱模型和條件脆弱模型均顯示脆弱方差有統(tǒng)計學意義，說明個體間異質(zhì)性存在，并且滿足比例風險假定;因 χ2= －163.2，P＜0.001，說明條件脆弱模型優(yōu)于共享脆弱模型，條件脆弱模型不但解釋了個體間的異質(zhì)性，同時也考慮個體內(nèi)的相關(guān)性。條件脆弱模型結(jié)論是不同組別癲癇病人復發(fā)風險差別有統(tǒng)計學意義( β=－0.990，P=0.0002)，陽性試驗組癲癇患者復發(fā)風險是對照組的0.372倍。脆弱項方差θ=2.27，P=0.0002，隨機效應有統(tǒng)計學意義，說明個體間異質(zhì)性存在，組內(nèi)相關(guān)性估計τ=θ/(2+θ)=0.53。Schoenfeld殘差圖(圖1)，說明該資料符合條件比例風險假定。

圖1 Schoenfeld殘差圖

討 論

臨床試驗研究中為探討不同藥物的治療效果，進行安全性評定癲癇病人的生存時間與生存質(zhì)量，常用事件發(fā)生等待時間來評價。含有刪失數(shù)據(jù)的生存分布比較中，對于單一生存事件分析已經(jīng)有各種各樣的參數(shù)和非參數(shù)的統(tǒng)計分析方法，但對于復發(fā)事件數(shù)據(jù)分析的模型仍處于研究中。

用于分析復發(fā)事件數(shù)據(jù)的模型有方差校正模型和脆弱模型等。方差校正模型〔14〕包括AG模型、條件模型和邊際模型，前提是假設數(shù)據(jù)由獨立事件組成，參數(shù)估計后用穩(wěn)健標準誤來校正復發(fā)事件帶來的標準誤的下偏。所以該模型不能解決由于個體間的異質(zhì)性引起的事件時間非獨立而導致的偏移效應。

脆弱模型是目前最受關(guān)注的一種解釋異質(zhì)性存在模型，假設異質(zhì)性是客觀存在的，用脆弱描述子組中個體生存時間之間存在的相關(guān)性，不同的個體有不同的脆弱，意味著某些個體或多或少地易于經(jīng)歷事件，相關(guān)個體被認為具有相同傾向，但是這種相同傾向的來源是未知的，不可測量或是未測量的;脆弱實際上是一個隨機效應，通過將隨機效應引入生存數(shù)據(jù)模型，就可以解釋個體間不能被觀測的協(xié)變量引起的異質(zhì)性，脆弱模型只能解釋異質(zhì)性而不能解釋個體內(nèi)的相關(guān)性。

本文探討用條件脆弱模型來分析復發(fā)事件數(shù)據(jù)，既可以解釋個體內(nèi)相關(guān)性，又能解釋個體間異質(zhì)性。條件脆弱模型可以比較客觀地分析癲癇復發(fā)數(shù)據(jù)，結(jié)果解釋客觀，軟件容易實現(xiàn)。研究表明如果事件依賴性存在，必須對不同的事件估計不同的基線風險，即根據(jù)事件數(shù)分層，若強行估計隨機效應，反而會帶來結(jié)果的偏移;如果異質(zhì)性存在，必須估計隨機效應，若估計了不同的基線風險函數(shù)同樣會帶來結(jié)果的偏移;對于一份實際數(shù)據(jù)，我們還沒有信息知道它是否存在異質(zhì)性和相依性，或者兩者同時存在時，建議采用條件脆弱模型進行分析，它可以捕獲相依性和異質(zhì)性效應。

1．Clayton DG．Some approaches to the analysis of recurrent event data．Statistical Methods in Medical Research，1994，3(3):244-262．

2．Kelly PJ，Lim LY．Survival analysis for recurrent event data:an application to childhood infectious diseases．Statistical in Medicine，2000，19(1):13-33．

3．Hyun JL，Liu JX，Marlene ML．Comparison of methods for analyzing recurrent events data:application to the emergency department visits of pediatric firearm victims．Accident Analysis ＆ Prevention，2007，39(2):290-299．

4．Clayton DG．A model for association in bivariate life tables and its application in epidemiological studies of familial tendency in chronic disease incidence．Biometrika，1978，65:141-151．

5．Roberto G．Gutierrez．Parametric frailty and shared frailty survival models．The Stata Journal．2002，2(1):22-44．

6．羅天娥，劉成芳，趙晉芳，等．共享伽瑪脆弱模型在癲癇復發(fā)的應用及實現(xiàn)．中國衛(wèi)生統(tǒng)計，2012，29(2):175-176．

7．Box-Steffensmeier JM，De Boef S．Repeated events survival models:the conditional frailty model．Statistical in Medicine，2006，25(20):3518-3533．

8．Janet MB，Suzanna DB，Kyle AJ．Event dependence and heterogeneity in duration:the conditional frailty model．Political Analysis，2007，15(3):237-256．

9．Duchateau L，Janssen P．The frailty model(Statistics for Biology and Health)．Springer:2007．

10．Rondeau V，Gonzalez JR．Frailtypack:a computer program for the analysis of correlated failure time data using penalized likelihood estimation．Computer Methods and Programs in Biomedicine，2005，80(2):154-164．

11．Janet M Box-Steffensmeier，Suzanna De Boef．Repeated events survival models:the conditional frailty model．Statistics in medicine，2006，25(20):3518-3533．

12．Philip Hougaard．Analysis of Multivariate Survival Data．New York:Springer-Verlag，2000．

13．Miscellanea．A note on scaled schoenfeld residuals for the proportional hazards model．Biometrika，2001，88(2):565-571．

14．余松林，向惠云．重復測量資料分析方法與SAS程序．北京:科學出版社，2003．

Application of Conditional Frailty Model for Epileptic Seizures and Implement

Luo Tiane，Zhao Jinfang，Sa Jian，et al．Department of Health Statistics，Shanxi Medical University(030001)，Taiyuan

ObjectiveTo explore the applications of conditional frailty models in recurrent event data about epileptic seizures and software implement.MethodsCollecting the epileptic seizures data，establishing the conditional frailty models，parameters of the estimate are calculated by penalized partial likehood function，and comparing with shared frailty models and Cox proportional hazard models.ResultsThe conditional frailty models are best account for the heterogeneity and event dependence，we can evaluate the clinical effect of epileptic seizures data using conditional frailty models，the explain is reasonable，the implement of software is easy.ConclusionConditional frailty models are useful for analyzing recurrent events data，it deserves to be widely applied．

Conditional frailty model;Penalized partial likehood;Epileptic seizures;Recurrent event data

(責任編輯:郭海強)