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

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

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

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

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

*:國(guó)家青年科學(xué)基金項(xiàng)目資助(編號(hào)81001294);山西醫(yī)科大學(xué)青年基金資助(02200913)

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

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

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

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

原理和方法

1.條件脆弱模型

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

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

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

2.參數(shù)估計(jì)

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

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

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

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

第二項(xiàng)是Gamma脆弱的懲罰項(xiàng):

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

3.基于Schoenfeld殘差比例風(fēng)險(xiǎn)假定的作圖法和檢驗(yàn)法

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

實(shí)例分析

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

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

表2 變量賦值表

用R 2.15.1版免費(fèi)軟件來(lái)擬合Cox比例危險(xiǎn)模型、共享Gamma脆弱模型和條件共享Gamma脆弱模型，程序編輯如下:

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

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

圖1 Schoenfeld殘差圖

討 論

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

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

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

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

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ù)發(fā)的應(yīng)用及實(shí)現(xiàn)．中國(guó)衛(wèi)生統(tǒng)計(jì)，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．余松林，向惠云．重復(fù)測(cè)量資料分析方法與SAS程序．北京:科學(xué)出版社，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

(責(zé)任編輯:郭海強(qiáng))