國家癌癥中心/中國醫(yī)學科學院北京協(xié)和醫(yī)學院腫瘤醫(yī)院(100021) 楊 召 王少明 梁 赫 喬友林 范金虎
·綜述·
競爭風險型數(shù)據(jù)統(tǒng)計分析理論研究進展
國家癌癥中心/中國醫(yī)學科學院北京協(xié)和醫(yī)學院腫瘤醫(yī)院(100021) 楊 召 王少明 梁 赫 喬友林 范金虎△
在醫(yī)學隨訪研究中,受試者在觀察期間常常由于出現(xiàn)其他結(jié)局事件而阻礙目標結(jié)局事件的發(fā)生或改變目標結(jié)局事件發(fā)生的概率。例如,在人群原因別死亡率的研究中目標結(jié)局事件為腫瘤相關(guān)死亡,但隨訪過程中常出現(xiàn)心腦血管相關(guān)死亡而阻止了腫瘤相關(guān)死亡事件的發(fā)生,即競爭風險(competing risk)[1]。在此理論框架下可能出現(xiàn):(1)目標結(jié)局事件(event of interest);(2)競爭事件(competing event);(3)刪失事件(cersored)。競爭風險在醫(yī)學研究中十分普遍[2-5],其最早可以追溯到1760年伯努利接種“疫苗”根除天花對人群死亡率的影響[6]。
本研究系統(tǒng)性地回顧競爭風險型數(shù)據(jù)的分析方法,主要從下面幾個方面對其研究進展進行系統(tǒng)性的描述:(1)競爭風險(competing risk)的定義;(2)常見競爭風險型數(shù)據(jù)分析方法的簡介:原因別風險模型(cause-specific hazard model)、累積風險模型(subdistribution hazard model)、加性模型(additivemodel)、基于虛擬觀測的回歸模型、混合模型(mixture)和縱向模型(verticalmodel);(3)競爭風險型數(shù)據(jù)統(tǒng)計分析的研究進展。
在隨訪研究中,目標結(jié)局事件短時間內(nèi)可能無法確定,因此為了評價研究因素對特定人群中目標結(jié)局事件的影響,所有受試者從基線開始隨訪,到發(fā)生目標結(jié)局事件或試驗結(jié)束。若隨訪期間未觀察到受試者出現(xiàn)目標結(jié)局事件則被定義為刪失(censored),如失訪(lost to follow-up),且在經(jīng)典的生存分析中假設刪失為“非信息性刪失(non-informative censored)”,即相同條件下受試者在某時刻發(fā)生刪失事件的概率和發(fā)生結(jié)局事件的概率相等[7]。若隨訪期間,受試者只出現(xiàn)目標結(jié)局事件,則可以通過經(jīng)典的分析方法,如KM法[8]和Cox等比例風險模型[7]進行分析,詳見圖1a。若隨訪期間出現(xiàn)多個結(jié)局事件,則為競爭風險[1],其根據(jù)對目標結(jié)局事件的影響程度分為兩類:
1.經(jīng)典競爭風險(classical competing risk)
經(jīng)典的競爭風險,指隨訪期間受試者出現(xiàn)多種互斥結(jié)局事件,即某一結(jié)局事件的發(fā)生常阻止其他結(jié)局事件的發(fā)生,詳見圖1b。譬如,在人群死亡研究中,若受試者出現(xiàn)腫瘤相關(guān)的死亡,則不會出現(xiàn)心腦血管疾病相關(guān)的死亡;反之亦然。
2.半式競爭風險(semi-competing risk)
半式競爭風險,也叫狀態(tài)轉(zhuǎn)移風險,指隨訪期間受試者由于出現(xiàn)其他結(jié)局事件而導致目標結(jié)局事件出現(xiàn)的概率發(fā)生改變,詳見圖1c。例如,在乳腺癌預后研究中,若目標結(jié)局事件為死亡,則受試者可能會因為出現(xiàn)復發(fā)轉(zhuǎn)移而使其死亡的風險發(fā)生改變。

圖1 競爭風險型數(shù)據(jù)模式圖(λi為風險函數(shù))
與經(jīng)典的生存分析方法相同,競爭風險模型數(shù)據(jù)的分析主要從以下三個方面:(1)時間效應,即目標結(jié)局事件的發(fā)病率如累積發(fā)病率(cumulative incidence)或發(fā)病密度(incidence densities);(2)干預效應,即不同干預措施對目標結(jié)局事件發(fā)病率的影響;(3)存在混雜因素時,干預因素對目標結(jié)局事件的累積發(fā)病率的影響。
為了能夠更加清晰明了地介紹下述常見競爭風險型數(shù)據(jù)的統(tǒng)計分析方法,做出以下假設:隨訪過程中會出現(xiàn)2個競爭事件(J,J=1提示發(fā)生結(jié)局事件1;J=2提示發(fā)生結(jié)局事件2),則該研究中所收集得到的競爭風險型數(shù)據(jù)為{ti,εi,εi×J,Xi},其中:ti為第i名受試者的最短觀察時間ti=min(TiJ=1,TiJ=2,Ci);εi為第i名受試者是否出現(xiàn)結(jié)局事件(εi=1,發(fā)生結(jié)局事件;εi=0,未發(fā)生結(jié)局事件);εi×J為第i名受試者發(fā)生結(jié)局事件的類型(εi×J=1,發(fā)生結(jié)局事件1;εi×J=2,發(fā)生結(jié)局事件2);Xi為第i名受試者協(xié)變量的狀態(tài)。那么,受試者發(fā)生結(jié)局事件的風險函數(shù)(hazard function)為λ(t);生存函數(shù)(survival function)為S(t);累積發(fā)病率函數(shù)(cumulative incidence function,CIF)為F(t)。
在經(jīng)典的生存分析中,統(tǒng)計學家將干預因素對目標結(jié)局事件發(fā)病率/死亡率的影響,通過Kaplan-Meier經(jīng)驗估計法轉(zhuǎn)化為干預因素對目標結(jié)局事件發(fā)生風險的影響,即:將原始無效假設H0:F(g=1)(t)=F(g=2)(t)轉(zhuǎn)化為H0′:λ(g=1)(t)=λ(g=2)(t)或H0′:S(g=1)(t)=S(g=2)(t)。這種轉(zhuǎn)化在競爭風險理論中卻不是等價的[9]。下文將著重從風險函數(shù)、生存函數(shù)和累積發(fā)病率函數(shù)的角度介紹幾種常見的競爭風險型數(shù)據(jù)分析方法:
1.風險函數(shù)為基礎的分析方法
(1)原因別風險模型(cause-specific hazard model)
原因別風險模型,作為Cox等比例風險模型的衍生模型,最早由Prentice[10]引進到競爭風險型數(shù)據(jù)的分析中,隨后在爭議中迅速發(fā)展并得到廣泛的應用[1]。該理論直接將競爭結(jié)局事件定義為刪失事件,然后依次對每個結(jié)局事件輪流擬合經(jīng)典的Cox等比例風險模型。此時,t時刻發(fā)生目標結(jié)局事件j風險函數(shù)

則t時刻的總體生存函數(shù)SCS(t)和t時間內(nèi)累積發(fā)病率函數(shù)分別為:

當存在干預因素時,該模型與經(jīng)典分析方法一樣需要做出以下假設:(1)干預因素對目標結(jié)局事件的影響不隨時間的變化而改變,即滿足等比例風險假設(proportional hazard assumption,PHA);(2)各結(jié)局事件之間的發(fā)生互不影響,即獨立性。根據(jù)經(jīng)典Cox等比例風險模型[7]的定義,則原因別風險模型即為:

為了能夠直觀形象地描述協(xié)變量與目標結(jié)局事件發(fā)生風險,exp(βj)=HRCS表示協(xié)變量每改變1個單位結(jié)局事件j發(fā)生風險的變化。此時,該模型不能直接評價協(xié)變量對目標結(jié)局事件累積發(fā)病率影響,究其原因主要是因為:結(jié)局事件1的累積發(fā)病率函數(shù)不僅僅與結(jié)局事件1的原因別風險函數(shù)有關(guān)系,而且與結(jié)局事件2的原因別風險函數(shù)有關(guān)系。當結(jié)局事件1的原因別風險函數(shù)增加時,其累積發(fā)病率函數(shù)不一定增加;反之亦然[9]。因此,原因別風險模型不能直接用來評價競爭風險型數(shù)據(jù)中目標結(jié)局事件的累積發(fā)病率;當利用HRCS描述協(xié)變量與目標結(jié)局事件發(fā)生風險之間的關(guān)聯(lián)時,需謹慎下結(jié)果。此外,該模型中回歸系數(shù)的估計采用部分似然函數(shù)(partial likelihood function)用最大似然估計方法得到,對模型的假設檢驗常采用似然比檢驗(maximum likelihood ratio test),其擬合過程可以借助標準統(tǒng)計分析軟件包如SAS、R或SPSS實現(xiàn)。當存在時依性協(xié)變量時,該模型和Cox模型一樣可拓展為時依性協(xié)變量的原因別風險模型(time-dependent causespecific hazard model),具體理論介紹詳見[11]。
(2)累積風險模型(subdistribution hazard model)
考慮到原因別風險模型存在前提假設嚴格且不能直接估計競爭風險型數(shù)據(jù)中目標結(jié)局事件累積發(fā)病率等缺點,F(xiàn)ine和Cray[12]提出利用累積風險模型直接估計協(xié)變量對目標結(jié)局事件累積發(fā)病率的影響。該理論充分考慮了競爭風險事件對目標結(jié)局事件的影響,其定義t時刻目標結(jié)局事件j的風險函數(shù)

則,t時間內(nèi)目標結(jié)局時間j的累積風險函數(shù)與經(jīng)典生存分析相同:

當存在干預因素或協(xié)變量時,該模型也需要滿足等比例風險假設和非信息性刪失。此時,累積風險模型為:

(3)加性模型(additivemodel)
在經(jīng)典的競爭風險理論假設中,結(jié)局事件之間的發(fā)生是互斥事件。Klein[13]參照Aalen加性模型提出各結(jié)局事件之間的累積發(fā)病率具有可加性,則在累積風險模型中,定義t時間內(nèi)結(jié)局事件j={1,2}的累積發(fā)病率(詳見公式(8))分別為:

根據(jù)加性模型理論,結(jié)局事件j={1,2}在t時間內(nèi)的累積發(fā)病率函數(shù)為:

當存在協(xié)變量時,根據(jù)累積風險模型理論:

式中,g(·)為轉(zhuǎn)換函數(shù)。若進一步將公式(14)代入(12):

此時,在公式(14)中出現(xiàn)了悖論,即在t=0時刻,公式兩側(cè)是不相等的。但是,該加性模型同時將各個競爭事件納入同一個模型中進行研究,可以直觀地比較協(xié)變量同時對多個結(jié)局事件的影響。該模型多采用偏似然估計法(partial likelihood approach)進行擬合,且可以在R中“survival”軟件包實現(xiàn)。此外,實際應用中研究者更注重協(xié)變量對目標結(jié)局事件遠期的影響[1]。
2.虛擬觀測(pseudo-observations)為基礎的分析方法
Andersen[14-15]認為競爭風險本質(zhì)上為多狀態(tài)模型的一種特殊形式,其將原因別風險看作為轉(zhuǎn)移風險,并引入潛在失效時間的分析方法。其主要思想是構(gòu)造虛擬觀測替換刪失個案,然后通過廣義線性模型擬合協(xié)變量對目標結(jié)局事件的影響。該理論定義t時間內(nèi)目標結(jié)局事件j累積發(fā)病率函數(shù)為

若假設隨訪期間有i={1,2,…,n}個刪失觀測,隨訪時間點為τt={τ1,τ2,…,τT},則虛擬觀測定義為:

其中,θit為τt時刻第i個刪失的虛擬觀測值,則分別為τt時間內(nèi)結(jié)局事件j的累積發(fā)病率函數(shù)、剔除刪失事件i之后τt時間內(nèi)結(jié)局事件j的累積發(fā)病率函數(shù)。當不存在刪失事件時,可以按照經(jīng)典廣義線性模型的方式分析各個協(xié)變量對目標結(jié)局事件的影響。定義連接函數(shù)g(·)為logit函數(shù),則在協(xié)變量X存在情況下,虛擬觀測是否發(fā)生結(jié)局事件的擬合過程即為多元logistic回歸分析的過程:

此時,針對各協(xié)變量對結(jié)局事件影響的解釋與logistic回歸分析相同。此外,Klein[16]發(fā)表了基于該理論的SAS和R統(tǒng)計分析代碼,通過案例研究證實該模型的擬合結(jié)果與累積風險模型的擬合效果相近。
3.聯(lián)合分布函數(shù)為基礎的分析方法
考慮到競爭風險型數(shù)據(jù)包含多結(jié)局事件及其對應的觀測時間,前者構(gòu)成結(jié)局事件的分布函數(shù),后者構(gòu)成觀測時間的分布函數(shù)。Larson和Dinse[17]提出利用聯(lián)合分布P(J,T)的來擬合目標結(jié)局事件P(J)的邊際分布及其觀測時間P(T|J)的條件分布,即混合模型(mixturemodel):

該聯(lián)合分布函數(shù)要求隨訪已經(jīng)結(jié)束,且需要利用EM(expectation-maximization)法修復刪失觀測,最后通過廣義線性模型進行估計。鑒于混合模型理論過于復雜且要求研究已經(jīng)完成,Nicolaie[18]提出使用觀測時間分布P(T)和觀察期內(nèi)結(jié)局事件J的條件分布P(J|T)擬合聯(lián)合分布P(T,J),即縱向模型(vertical model):

該理論定義結(jié)局事件J發(fā)生的條件概率P(J=j|T=t)為相對原因別風險(relative cause-specific hazards,csRH),則在上述原因別風險模型的基礎上csRH(t)可描述為:

值得注意的是,csRH與風險比(hzarad ratio,HR)是相同的,其也要求刪失為非信息性刪失,即:
式中,P(J=j|T=t,C≥t)和P(J=j|T=t,T≤t)分別表示t時刻受試者發(fā)生右刪失和左刪失的概率。若假設任意結(jié)局事件發(fā)生的時間為τt={τ1,τ2,…,τT,根據(jù)公式(1)和公式(20)csRHj(τt)可以表達為:

若假定csRHj(τt)為連續(xù)變量,則需要引入平滑函數(shù)B(τt)擬合每個時間點的csRHj(τt),此時通過多元logistic模型建立csRHj(τt)和B(τt)之間的關(guān)聯(lián),即為:

此外,在原因別風險模型的理論框架下,建立csRHj(τt)與結(jié)局事件j在t時間內(nèi)累積發(fā)病率的關(guān)系如下:

式中,f′為全部結(jié)局事件發(fā)生時間的密度函數(shù)。此時,csRHj(t)可以描述為全部結(jié)局時間的分布密度函數(shù)在累積發(fā)病率函數(shù)中的比例。當存在協(xié)變量時,該模型可以進一步拓展,考慮到其擬合過程較為復雜,因此此處不再詳述[18]。
縱觀競爭風險型數(shù)據(jù)統(tǒng)計分析發(fā)展史,我們不難發(fā)現(xiàn)競爭風險的本質(zhì)是將刪失事件進行更細致的歸類和處理,其主要的理論分析框架仍為生存分析。然而,傳統(tǒng)的生存分析,沒有充分地認識到競爭事件對目標結(jié)局事件的影響或直接將競爭事件視為刪失事件,導致結(jié)局事件的風險函數(shù)λ(t)出現(xiàn)有偏估計,繼而引起生存函數(shù)S(t)和累積發(fā)病率函數(shù)F(t)出現(xiàn)有偏估計。現(xiàn)有競爭風險型數(shù)據(jù)分析方法,在考慮競爭事件的同時,直接研究干預因素對目標結(jié)局事件發(fā)生率的影響,改變了經(jīng)典生存分析中將目標結(jié)局事件發(fā)生率的研究轉(zhuǎn)換為發(fā)生風險研究的思路,更加直觀、真實地評價干預對目標結(jié)局事件發(fā)生率的應用。但是,復雜的理論研究基礎導致其缺少標準的統(tǒng)計分析軟件包,進而阻礙了其在實際研究中的應用。值得慶幸的是,近年來國內(nèi)越來越多的研究也開始關(guān)注競爭風險型數(shù)據(jù)的分析[19-20]。
[1]Lau B,Cole SR,Gange SJ.Competing risk regressionmodels for epidem iologic data.Am JEpidem iol,2009,170:244-256.
[2]Austin PC,Lee DS,F(xiàn)ine JP.Introduction to the Analysis of Survival Data in the Presence of Competing Risks.Circulation,2016,133:601-609.
[3]de Glas NA,Kiderlen M,Vandenbroucke JP,et al.Perform ing Survival Analyses in the Presence of Competing Risks:A Clinical Example in Older Breast Cancer Patients.JNatl Cancer Inst,2016:108.
[4]Suri RM,ClavelMA,Schaff HV,etal.Effectof RecurrentM itral Regurgitation Follow ing Degenerative M itral Valve Repair:Long-Term Analysis of Competing Outcomes.J Am Coll Cardiol,2016,67:488-498.
[5]Ryser MD,Worni M,Turner EL,et al.Outcomes of Active Surveillance for Ductal Carcinoma in Situ:A Computational Risk Analysis.J Natl Cancer Inst,2016:108.
[6]Chiang CL.Competing risks in mortality analysis.Annu Rev Public Health,1991,12:281-307.
[7]David CR.Regression models and life tables(w ith discussion).Journal of the Royal Statistical Society,1972,34:187-220.
[8]Kaplan EL,Meier P.Nonparametric estimation from incomplete observations.Journal of the American statistical association,1958,53:457-481.
[9]Gray RJ.A class of K-sample tests for comparing the cumulative incidence of a competing risk.The Annals of statistics,1988:1141-1154.
[10]Prentice RL,Kalbfleisch JD,Peterson AV,Jr.,et al.The analysis of failure times in the presence of competing risks.Biometrics,1978,34:541-554.
[11]Sun Y,Hyun S,Gilbert P.Testing and estimation of time-varying cause-specific hazard ratios w ith covariate adjustment.Biometrics,2008,64:1070-1079.
[12]Fine JP,Gray RJ.A proportional hazardsmodel for the subdistribution of a competing risk.Journal of the American statistical association,1999,94:496-509.
[13]Klein JP.Modelling competing risks in cancer studies.Stat Med,2006,25:1015-1034.
[14]Andersen PK,Abildstrom SZ,Rosthφj S.Competing risks as a multistatemodel.Statisticalmethods inmedical research,2002,11:203-215.
[15]Andersen PK,Klein JP,Rosthφj S.Generalised linearmodels for correlated pseudo-observations,with applications to multi-state models.Biometrika,2003,90:15-27.
[16]Klein JP,Gerster M,Andersen PK,et al.SAS and R functions to compute pseudo-values for censored data regression.Comput Methods Programs Biomed,2008,89:289-300.
[17]Larson MG,Dinse GE.A mixturemodel for the regression analysisof competing risks data.Applied statistics,1985:201-211.
[18]Nicolaie MA,van Houwelingen HC,Putter H.Vertical modeling:a pattern mixture approach for competing risks modeling.Stat Med,2010,29:1190-1205.
[19]江一濤,胡海蘭,魏巧玲,等.競爭風險模型的發(fā)展與應用.中國衛(wèi)生統(tǒng)計,2009,26:445-447.
[20]肖媛媛,許傳志,趙耐青.常用生存分析模型及其對時依性協(xié)變量效應的估計方法.中國衛(wèi)生統(tǒng)計,2016:543-547+552.
(責任編輯:劉 壯)
△通信作者:范金虎,E-mail:fanjh@cicams.ac.cn