參數(shù)正則穩(wěn)態(tài)脆弱模型在臨床復(fù)發(fā)事件數(shù)據(jù)分析中的應(yīng)用*

郭 強(qiáng) 王晶晶 劉桂芬 羅天娥△

參數(shù)正則穩(wěn)態(tài)脆弱模型在臨床復(fù)發(fā)事件數(shù)據(jù)分析中的應(yīng)用*

郭 強(qiáng)1王晶晶2劉桂芬1羅天娥1△

目的 探討參數(shù)正則穩(wěn)態(tài)脆弱模型在臨床復(fù)發(fā)事件數(shù)據(jù)中的應(yīng)用。方法 以慢性肉芽腫(CGD)患者為例構(gòu)建參數(shù)正則穩(wěn)態(tài)脆弱模型，假定復(fù)發(fā)事件的基線分布服從威布爾分布、指數(shù)分布或Gompertz分布，脆弱項(xiàng)服從正則穩(wěn)態(tài)脆弱分布，采用邊際極大似然估計(jì)(MMLE)來(lái)實(shí)現(xiàn)參數(shù)估計(jì)并進(jìn)行不同條件下模型對(duì)比分析。結(jié)果 參數(shù)正則穩(wěn)態(tài)脆弱模型既考慮了個(gè)體內(nèi)復(fù)發(fā)時(shí)間的相關(guān)性，又考慮了不同個(gè)體間的異質(zhì)性，可以用來(lái)分析慢性肉芽腫患者反復(fù)住院的影響因素，結(jié)果解釋合理，軟件實(shí)現(xiàn)便捷。結(jié)論 參數(shù)正則穩(wěn)態(tài)脆弱模型適合分析臨床復(fù)發(fā)數(shù)據(jù)并進(jìn)行臨床療效評(píng)價(jià)或者影響因素研究。

復(fù)發(fā)事件數(shù)據(jù) 脆弱模型 正則穩(wěn)態(tài)脆弱分布 威布爾分布

在臨床醫(yī)學(xué)研究中，Cox[1-2]比例風(fēng)險(xiǎn)模型通常是進(jìn)行時(shí)間-事件數(shù)據(jù)研究的經(jīng)典模型，其理論假設(shè)是研究對(duì)象間是同質(zhì)的。但對(duì)于同一個(gè)體多次重復(fù)觀測(cè)所獲得的復(fù)發(fā)事件數(shù)據(jù)(recurrent event data)，其往往具有非獨(dú)立性或異質(zhì)性的特點(diǎn)，若直接采用Cox比例風(fēng)險(xiǎn)回歸模型分析，會(huì)造成參數(shù)估計(jì)的偏差。Hougarrd和Aalen[3-4]在相關(guān)研究中發(fā)現(xiàn)，忽略個(gè)體的異質(zhì)性將導(dǎo)致估計(jì)的相對(duì)危險(xiǎn)度有偏高的傾向。Pickles等[5]在對(duì)脆弱模型的回顧性研究中發(fā)現(xiàn)，模型中忽略個(gè)體間異質(zhì)性的影響會(huì)導(dǎo)致協(xié)變量系數(shù)的估計(jì)值趨向于0。目前脆弱模型(frailty model)[6-7]廣泛應(yīng)用于復(fù)發(fā)事件數(shù)據(jù)的分析中，此類模型充分考慮了復(fù)發(fā)事件數(shù)據(jù)間的非獨(dú)立性特點(diǎn)，可較好地解釋潛在協(xié)變量所引起的異質(zhì)性問(wèn)題，提高了檢驗(yàn)精度[8]。本研究擬構(gòu)建不同基線分布的正則穩(wěn)態(tài)脆弱模型來(lái)對(duì)慢性肉芽腫患者反復(fù)住院的影響因素進(jìn)行分析并用R3.3.1軟件實(shí)現(xiàn)。

原理和方法

1.參數(shù)正則穩(wěn)態(tài)脆弱模型

脆弱模型[9-10]是在一般比例風(fēng)險(xiǎn)模型中引入脆弱因子，用于解釋異質(zhì)性或相關(guān)性問(wèn)題的統(tǒng)計(jì)分析模型。正則穩(wěn)態(tài)脆弱分布是Hougaard 于1986年提出的一種脆弱分布類型，若令ui為獨(dú)立同分布的正則穩(wěn)態(tài)脆弱變量，則該變量的概率密度函數(shù)如下[11]：

(1)

(2)

結(jié)合前面正則穩(wěn)態(tài)脆弱分布Laplace轉(zhuǎn)換式，可以得出第i組威布爾正則穩(wěn)態(tài)脆弱模型的聯(lián)合生存函數(shù)為：

(3)

Hx,c(tij)=H0(tij)exp(βXij)為生存函數(shù)中的非條件部分。指數(shù)分布和Gompertz分布正則穩(wěn)態(tài)脆弱模型的聯(lián)合生存函數(shù)，其公式除內(nèi)部累積風(fēng)險(xiǎn)函數(shù)H(t)不同外，形式上與公式(3)完全相同。

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

Fisher(1921)提出似然概念以來(lái)，關(guān)于參數(shù)估計(jì)的似然理論和方法得到了極大的發(fā)展和廣泛應(yīng)用[14]。本文采用邊際極大似然估計(jì)(MMLE)來(lái)實(shí)現(xiàn)對(duì)以上模型的參數(shù)估計(jì)[13]。假設(shè)有i=1,2,…,n個(gè)組，每個(gè)組有j=1,2,…,q個(gè)樣本，則可用q個(gè)樣本的概率乘積來(lái)表示第i個(gè)組的概率函數(shù)，而聯(lián)合密度函數(shù)可用Laplace轉(zhuǎn)換后的第n次求導(dǎo)得出。假定基線風(fēng)險(xiǎn)函數(shù)服從威布爾分布，對(duì)所有的個(gè)體或組求和后的邊際似然函數(shù)：

L(λ,υ,β,θ)=

(4)

若假定基線風(fēng)險(xiǎn)函數(shù)服從指數(shù)分布，對(duì)所有的個(gè)體或組求和后的邊際似然函數(shù)記作：

L(λ,β,θ)=

(5)

若假定基線風(fēng)險(xiǎn)函數(shù)服從Gompertz分布，對(duì)所有的個(gè)體或組求和后的邊際似然函數(shù)記作：

L(λ,φ,β,θ)=

(6)

L(D)表示Laplace轉(zhuǎn)換后的第D次求導(dǎo)。

實(shí)例分析

該資料是126例慢性肉芽腫患者反復(fù)感染的數(shù)據(jù)[15]，欲研究伽瑪干擾素對(duì)慢性肉芽腫疾病的療效，患者被隨機(jī)分為兩組，分別給予伽瑪干擾素和安慰劑治療；選入的變量有患者的年齡、性別、身高、體重、遺傳方式等因素。部分?jǐn)?shù)據(jù)見表1，賦值情況見表2。

表1 慢性肉芽腫患者資料部分?jǐn)?shù)據(jù)

表2 慢性肉芽腫患者資料變量賦值表

1.慢性肉芽腫患者資料的分析

由于該資料是慢性肉芽腫患者反復(fù)感染的數(shù)據(jù)，構(gòu)建模型時(shí)需要考慮個(gè)體重復(fù)觀測(cè)數(shù)據(jù)間的相關(guān)性，建立參數(shù)正則穩(wěn)態(tài)脆弱模型：

β6cortico+β7prophy+β8hospital+β9rIFN+ui)

其中，患者編號(hào)i=1,2,…,126，重復(fù)觀測(cè)最多的次數(shù)j=1,2,…,8。每個(gè)患者的異質(zhì)性因子為ui，假定為正則穩(wěn)態(tài)分布，以此來(lái)反映每位患者重復(fù)觀測(cè)值間的相關(guān)性；在給定脆弱因子ui的情況下，可以假定患者重復(fù)觀測(cè)的時(shí)間是相互獨(dú)立的；假定基線分布為威布爾分布、指數(shù)分布或Gompertz分布。

2.慢性肉芽腫(CGD)資料參數(shù)脆弱模型分析

相關(guān)數(shù)據(jù)采用R 3.3.1軟件進(jìn)行分析，需加載R包有survival，parfm和eha。分別用Weibull比例風(fēng)險(xiǎn)(Weibull PH)模型、威布爾正則穩(wěn)態(tài)脆弱(W-PS)模型、指數(shù)正則穩(wěn)態(tài)脆弱模型(exponential-PS)及Gompertz正則穩(wěn)態(tài)脆弱模型(Gompertz-PS)等模型進(jìn)行擬合并比較，參數(shù)估計(jì)結(jié)果見表3。

表3 慢性肉芽腫患者數(shù)據(jù)資料的幾種模型比較

*：*表示P<0.05。

從模型構(gòu)建上來(lái)說(shuō)，Cox模型為半?yún)?shù)模型，與參數(shù)模型比較不太合理，因?yàn)槟Ｐ凸烙?jì)的參數(shù)個(gè)數(shù)不同，所以模型的AIC和BIC的大小沒有可比性。而Weibull比例風(fēng)險(xiǎn)(Weibull PH)模型是假定Cox模型中基線風(fēng)險(xiǎn)服從Weibull分布，為參數(shù)模型。因此選用Weibull PH模型與以上三種參數(shù)正則穩(wěn)態(tài)脆弱模型進(jìn)行比較。結(jié)果顯示，與脆弱模型相比，Weibull PH模型參數(shù)估計(jì)的絕對(duì)值和方差較脆弱模型普遍偏低；這是由于Weibull PH模型沒有考慮復(fù)發(fā)事件資料中不同時(shí)間上存在的相關(guān)性和異質(zhì)性問(wèn)題，導(dǎo)致了參數(shù)估計(jì)的偏差。這也從側(cè)面說(shuō)明，脆弱模型至少可以解釋資料中的部分未檢測(cè)到的相關(guān)性和異質(zhì)性。

將Weibull PH模型與其他三種脆弱模型擬合結(jié)果比較，-2LL分別為3.90，5.872和5.874，似然比異質(zhì)性檢驗(yàn)有統(tǒng)計(jì)學(xué)意義(P值分別為0.048，0.015和0.015)，資料存在異質(zhì)性，提示對(duì)該資料擬合脆弱模型更好；個(gè)體間相關(guān)性τ分別為0.068、0.072和0.117，子組間異質(zhì)性θ=1-τ，θ越小異質(zhì)性越大。

考慮三個(gè)協(xié)變量age、cortico和rIFN的脆弱模型擬合結(jié)果表明，幾種脆弱模型估計(jì)結(jié)果基本一致。根據(jù)AIC和BIC值越小模型擬合越好，該數(shù)據(jù)選用指數(shù)正則穩(wěn)態(tài)脆弱模型(exponential-PS)，結(jié)果解釋為：患者的年齡(age)和是否使用干擾素(rIFN)是影響該病患者反復(fù)感染的主要因素，使用干擾素治療的病人不容易出現(xiàn)反復(fù)感染，其反復(fù)感染的風(fēng)險(xiǎn)是安慰劑治療者的0.355倍，其相對(duì)危險(xiǎn)度(RR)的95%可信區(qū)間為(0.192，0.656)；另外患者年齡越小，發(fā)生反復(fù)感染的可能性越大，在其他變量不變的情況下，患者年齡每增一歲，其發(fā)生反復(fù)感染的風(fēng)險(xiǎn)變?yōu)樵瓉?lái)的0.968倍，RR的95%可信區(qū)間為(0.938，0.998)。具體結(jié)果見表4。

表4 慢性肉芽腫患者資料exponential-PS模型參數(shù)估計(jì)

以患者編號(hào)為橫軸，每個(gè)個(gè)體的脆弱值為縱軸繪制圖1，可見該資料個(gè)體間脆弱因子變異較大，即說(shuō)明資料具有異質(zhì)性，構(gòu)建exponential-PS模型是合適的。

圖1 指數(shù)正則穩(wěn)態(tài)脆弱模型的脆弱預(yù)測(cè)圖

討 論

1.脆弱模型是Cox回歸模型的擴(kuò)展，可有效處理臨床研究復(fù)發(fā)時(shí)間數(shù)據(jù)中的相關(guān)性和異質(zhì)性問(wèn)題。本文通過(guò)構(gòu)建參數(shù)正則穩(wěn)態(tài)脆弱模型，對(duì)參數(shù)采用MMLE方法進(jìn)行估計(jì)。并將該方法應(yīng)用于慢性肉芽腫(CGD)患者反復(fù)感染的影響因素分析，結(jié)果顯示參數(shù)正則穩(wěn)態(tài)脆弱模型在分析中既考慮了個(gè)體內(nèi)復(fù)發(fā)時(shí)間的相關(guān)性，又考慮了不同個(gè)體的異質(zhì)性，其結(jié)果解釋合理，軟件實(shí)現(xiàn)便捷。

2.參數(shù)正則穩(wěn)態(tài)脆弱模型可以指定不同的基線風(fēng)險(xiǎn)函數(shù)，而且正則穩(wěn)態(tài)脆弱分布也有其優(yōu)勢(shì)，若個(gè)體條件風(fēng)險(xiǎn)比和總體風(fēng)險(xiǎn)比都成比例，并且隨著時(shí)間變量的不斷延伸，總體風(fēng)險(xiǎn)比逐漸趨近于1，選擇此分布較好。參數(shù)脆弱模型由于事先假定了回歸方程的分布類型，其參數(shù)估計(jì)簡(jiǎn)單方便，結(jié)果可靠性較高，而且計(jì)算結(jié)果可以外延。

3.待深入研究的問(wèn)題：對(duì)于臨床復(fù)發(fā)事件數(shù)據(jù)構(gòu)建脆弱模型進(jìn)行分析時(shí)，脆弱分布的選擇往往基于數(shù)理特點(diǎn)以及分布函數(shù)Laplace轉(zhuǎn)換式是否簡(jiǎn)單，而不是根據(jù)數(shù)據(jù)資料本身的相關(guān)結(jié)構(gòu)特點(diǎn)去選擇模型。因此，還需進(jìn)一步探討選擇脆弱分布是否合理的檢驗(yàn)方法；對(duì)于正則穩(wěn)態(tài)脆弱回歸模型，目前的一些軟件缺乏較好的模型診斷方法，是今后需要解決的難題。

[1]Cox D.Regression models and life-tables.The Journal of the Royal Statistical Society(B),1972,34(2):187-220.

[2]肖媛媛,許傳志,趙耐青.常用生存分析模型及其對(duì)時(shí)依性協(xié)變量效應(yīng)的估計(jì)方法.中國(guó)衛(wèi)生統(tǒng)計(jì),2016,33(3):543-547.

[3]Hougaard P.Frailty models for survival data.Lifetime Data Analysis,1995,1(3):255-273.

[4]Aalen O.Modeling heterogeneity in survival analysis by the compound Poisson distribution.Annals of Applied Probability,1992,2(4):951-972.

[5]Pickles A,Crouchley R.A comparison of frailty models for multivariate survival data.Statistics in Medicine,1995,14(13):1447-1461.

[6]Zhang Y,Chen MH,Ibrahim JG,et al.Bayesian gamma frailty models for survival data with semi-competing risks and treatment switching.Lifetime Data Analysis,2014,20(1):76-105.

[7]Gerster M,Madsen M,Andersen PK.Matched survival data in a co-twin control design.Lifetime Data Analysis,2014,20(1):38-50.

[8]朱玉,梅楊,李杰,等.Cox比例風(fēng)險(xiǎn)Frailty模型簡(jiǎn)介與軟件實(shí)現(xiàn).中國(guó)衛(wèi)生統(tǒng)計(jì),2014,31(3):527-529.

[9]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(1):141-151.

[10]羅天娥,劉成芳,趙晉芳,等.共享伽瑪脆弱模型在癲癇復(fù)發(fā)的應(yīng)用及實(shí)現(xiàn).中國(guó)衛(wèi)生統(tǒng)計(jì),2012,29(2):175-176.

[11]Hanagal DD,Dabade AD.Compound negative binomial shared frailty models for bivariate survival data.Statistics and Probability Letters,2013,83(83):2507-2515.

[12]王晶晶.非獨(dú)立生存數(shù)據(jù)正則穩(wěn)態(tài)脆弱模型分析及應(yīng)用.山西醫(yī)科大學(xué),2012:7-13.

[13]Duchateau L,Janssen P.The frailty model.Springer:New York,2008,33-286.

[14]王寧寧,徐淑一,方積乾.從經(jīng)典似然到等級(jí)似然的理論概述和應(yīng)用.中國(guó)衛(wèi)生統(tǒng)計(jì),2016,2:364-367+369.

[15]Therneau TM,Grambsch PM.Modeling Survival Data:Extending the Cox Model.Journal of the American Statistical Association,2002,44(457):85-86.

(責(zé)任編輯：劉 壯)

Application of Parametric Positive Stable Frailty Model for Analysis of Clinical Recurrent Event Data

Guo Qiang,Wang Jingjing,Liu Guifen,et al

(ShanxiMedicalUniversity(030001),Taiyuan)

Objective To explore the application of parametric positive stable frailty model in analysis of clinical recurrent event data.Methods Taking the data of patients with chronic granulomatous disease(CGD)as an example,establishing the parametric positive stable frailty model,assuming that the distribution of baseline is Weibull distribution,exponential distribution and Gompertz distribution,frailty factor is positive stable distribution.Parameters were estimated by using the Marginal Maximum Likelihood Estimation(MMLE).The models in different baseline distribution were compared and analyzed.Results The parametric positive stable frailty model not only regard the dependent within groups,but also consider the heterogeneity between groups in recurrent event data.This model can be used to analyze repeatedly hospitalized factors of patients with chronic granulomatous disease(CGD).The results are easy to implement with software and its explain is reasonable.Conclusions The parametric positive stable frailty model is suitable for analyzing clinical recurrent event data and studying clinical evaluation or influencing factors.

Recurrent event data;Frailty model;Positive stable distribution;Weibull distribution;

* 國(guó)家青年科學(xué)基金項(xiàng)目資助(81001294)；國(guó)家自然科學(xué)基金項(xiàng)目資助(81172774)

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

2.山西省眼科醫(yī)院

△ 通信作者：羅天娥，E-mail:luotiane1977@163.com