參數(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分布，脆弱項服從正則穩(wěn)態(tài)脆弱分布，采用邊際極大似然估計(MMLE)來實現(xiàn)參數(shù)估計并進(jìn)行不同條件下模型對比分析。結(jié)果 參數(shù)正則穩(wěn)態(tài)脆弱模型既考慮了個體內(nèi)復(fù)發(fā)時間的相關(guān)性，又考慮了不同個體間的異質(zhì)性，可以用來分析慢性肉芽腫患者反復(fù)住院的影響因素，結(jié)果解釋合理，軟件實現(xiàn)便捷。結(jié)論 參數(shù)正則穩(wěn)態(tài)脆弱模型適合分析臨床復(fù)發(fā)數(shù)據(jù)并進(jìn)行臨床療效評價或者影響因素研究。

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

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

原理和方法

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

脆弱模型[9-10]是在一般比例風(fēng)險模型中引入脆弱因子，用于解釋異質(zhì)性或相關(guān)性問題的統(tǒng)計分析模型。正則穩(wěn)態(tài)脆弱分布是Hougaard 于1986年提出的一種脆弱分布類型，若令ui為獨立同分布的正則穩(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)險函數(shù)H(t)不同外，形式上與公式(3)完全相同。

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

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

L(λ,υ,β,θ)=

(4)

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

L(λ,β,θ)=

(5)

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

L(λ,φ,β,θ)=

(6)

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

實例分析

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

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

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

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

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

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

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

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

相關(guān)數(shù)據(jù)采用R 3.3.1軟件進(jìn)行分析，需加載R包有survival，parfm和eha。分別用Weibull比例風(fēng)險(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ù)估計結(jié)果見表3。

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

*：*表示P<0.05。

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

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

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

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

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

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

討 論

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

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

3.待深入研究的問題：對于臨床復(fù)發(fā)事件數(shù)據(jù)構(gòu)建脆弱模型進(jìn)行分析時，脆弱分布的選擇往往基于數(shù)理特點以及分布函數(shù)Laplace轉(zhuǎn)換式是否簡單，而不是根據(jù)數(shù)據(jù)資料本身的相關(guān)結(jié)構(gòu)特點去選擇模型。因此，還需進(jìn)一步探討選擇脆弱分布是否合理的檢驗方法；對于正則穩(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]肖媛媛,許傳志,趙耐青.常用生存分析模型及其對時依性協(xié)變量效應(yīng)的估計方法.中國衛(wèi)生統(tǒng)計,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)險Frailty模型簡介與軟件實現(xiàn).中國衛(wèi)生統(tǒng)計,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)用及實現(xiàn).中國衛(wèi)生統(tǒng)計,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]王晶晶.非獨立生存數(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)典似然到等級似然的理論概述和應(yīng)用.中國衛(wèi)生統(tǒng)計,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;

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

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

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

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