李海彬 李 霞 王安心 陶麗新 劉 龍 陳斯鵬 郭 晉 郭秀花△
·計算機應用·
競爭風險模型及其在Stata軟件實現*
李海彬1,2#李 霞1,2#王安心1,2陶麗新1,2劉 龍1,2陳斯鵬1,2郭 晉1,2郭秀花1,2△
目的介紹競爭風險模型及在Stata軟件中的實現。方法使用競爭風險模型對白血病患者接受骨髓干細胞移植進行建模,展示該模型在Stata軟件中的實現過程并對最終結果進行解釋。結果Stata軟件“Competing-risks regression”菜單和“stcrreg”命令可以實現競爭風險模型。實例應用中,當考慮移植相關死亡事件作為疾病復發的競爭事件后,基于競爭風險模型得到病人的疾病進展階段與復發有關。結論競爭風險事件廣泛存在于臨床研究中,Stata軟件可以簡單靈活地實現這一過程。
競爭事件 競爭風險模型 Stata軟件
經典生存分析中,研究者往往只關注研究對象隨訪過程中特定事件的發生。然而在現實世界中,研究對象不僅經歷一種類型事件,不同類型結局事件相互影響,即形成競爭事件。處理含有競爭事件數據的統計模型稱作“競爭分析模型”。1999年Fine和Gray提出部分分布的半參數比例風險模型,該模型在考慮到競爭事件的條件下,探討其他影響因素對結局事件的作用[1],彌補了傳統生存分析的不足。目前競爭風險模型被廣泛應用于臨床試驗、流行病學調查等研究中。然而國內文獻中介紹競爭風險模型及Stata程序實現相對較少,本文將結合臨床實例介紹競爭風險模型及如何使用Stata軟件進行建模,旨在更好地為研究中提高對競爭事件的認識和通過Stata軟件實現這一過程操作提供幫助。
1.競爭事件
競爭風險(competing risks)是指研究對象出現感興趣事件(interesting event)的同時還會出現其他結局事件,這些結局事件將阻止感興趣事件的出現或使其發生的概率降低,各結局事件間形成所謂的"競爭"關系,這一系列事件稱作競爭事件(competing risk events)[2]。例如在白血病患者接受骨髓干細胞移植后復發風險的研究中,應當將移植相關死亡事件作為競爭事件處理[3],倘若忽視了該競爭事件的存在,而簡單的將競爭事件等同于右刪失(right censoring),使用傳統K-M法估計終點事件發生概率、Cox模型進行多因素分析,將產生偏差[4]。
2.競爭風險模型
1999年Fine和Gray[1]提出的部分分布的半參數比例風險模型(Fine-Gray模型)就是一種競爭風險模型(competing risksmodel),使用累積風險函數(cumulative incidence function,CIF)來估計結局事件的累積發生概率。在該模型中,t時刻發生事件j風險定義如下[1]:

其中λ0(t)是事件j的部分分布基準風險(baseline hazard of the sub-distribution);Z是協變量;βj′Z協變量的回歸系數;exp{βj′Z}是部分分布相對風險(sub-distribution hazard ratio)。部分分布風險模型偏似然估計定義如下:

風險集Rj定義如下:

風險集是由在時間t沒有發生感興趣事件個體和在時間前發生競爭事件的個體構成。因此經歷其他類型事件的個體仍在風險集中。權重定義如下:

G(t)是Kaplan-Meier方法估計的生存函數。時間t前沒有發生競爭事件的個體,在時間t發生感興趣事件有相等的權重(wi=1);ti<t時,發生競爭事件個體權重wi<1。
競爭風險模型通過以下公式建立λj(t,Z)與CIF的關系:

因此,通過競爭風險模型能夠估計出協變量Z相應的系數β′j,聯合累積部分分布基準風險預測特定時間個體發生感興趣事件發生概率,即絕對風險。
3.競爭風險模型軟件實現情況
競爭風險模型可以通過R軟件“cmprsk”程序包實現[5],可參考文獻[3、6][3,6],或SAS軟件宏“%Cum Inc”[7]、“%CIF”[8]實現以及SAS9.4版本“PRCO PHREG”[9]。考慮到使用上述過程需要有一定的編程基礎,而Stata軟件可以使用“stcrreg”命令擬合競爭風險模型[9],也可以簡單通過Stata軟件菜單操作完成,極大的方便了臨床醫生及其他研究者對競爭風險模型的實現。本文結合Stata 14.0軟件菜單操作和命令語句結合臨床實例擬合競爭分析模型,為更好地通過Stata軟件實現競爭風險模型分析提供幫助。
本文的實例數據是關于白血病患者接受骨髓干細胞移植后復發風險的研究[3],下載地址http://www.stat.unipg.it/luca/R/,在Competing Risks Analysis目錄,下載bmtcrr.xls文件。研究的主要目的是考慮移植相關死亡即競爭事件后,探索病人復發的累積發生率以及復發的影響因素。該研究包含有177例接受骨髓干細胞移植的急性白血病患者。表1為變量命名、賦值及研究對象的基本特征。

表1 變量命名、賦值及統計描述
1.數據集導入與變量命名、賦值
Stata軟件可以直接導入Excel文件,點擊菜單File—Import—Excel spreadsheet(*.xls*.xlsx)—Browse—選擇bm tcrr.xls文件—勾選Import first row as variable names—OK。數據導入結果如圖1所示。原始數據是字符型,需要將變量賦值。命令如下:

圖1 數據導入及原始數據結構

2.擬合競爭風險模型
選擇菜單Statistics—Survival analysis—Regression models—Competing-risks regression—Survival setting。對生存資料進行定義,如圖2所示,Time variable選擇生存時間Ftime,Failure event中Failure variable生存結局變量選擇Status。下一步是對競爭風險模型的擬合,Independent variables即是選擇自變量,在下拉菜單中選擇相應的變量(Age,sex,disease,source,phase)。Competing-risks events即是對競爭事件的定義,Variable選擇Status,Value輸入2,即移植相關死亡。采用“stcurve”命令可以繪制疾病進展階段(phase)的累積發生曲線。命令如下:


圖2 生存數據定義及競爭風險模型菜單操作
1.結果
圖3為生存數據定義結果,顯示177名病人中56人發生復發事件,共觀察到299.08人年(3588.94/12),最長隨訪時間為131.77個月。

圖3 生存數據定義結果
圖4 為競爭風險模型的結果。右側“No.of obs、No.of subjects、No.failed、No.competing、No.censored"分別是對觀察個體、研究對象、出現結局事件、競爭事件、刪失的個數的統計?!癓og pseudolikelihood=-266.52444,Wald chi2(7)=22.93,Prob>chi2=0.0018"是對競爭風險模型整體評價,P<0.05,具有統計學意義。圖4中的下半部分是多因素回歸結果。從左列到右列分別是變量、部分分布風險比(SHR)、標準誤、統計量(Z)、P值和SHR的95%置信區間(95%CI)。

圖4 競爭風險模型的結果
2.結果解釋
以性別變量為例,在考慮移植相關死亡為競爭風險事件的前提下,控制年齡、疾病類型、進展階段、移植類型后,女性白血病患者骨髓移植后發生復發的風險是男性的0.97倍(95%CI:0.55-1.71)。95%CI包含1,同時P=0.904>0.05,差異不具有統計學意義。其它變量的解釋與性別變量解釋類似,但疾病進展階段變量需要注意,為無序變量,采用啞變量處理。以狀態“Relapse”為參考組,考慮競爭事件后“CR1”組發生復發的風險降低67%(SHR,0.33,95%CI:0.16-0.70,P=0.004);“CR2”組發生復發的風險降低64%(SHR,0.36,95%CI:0.18-0.73,P=0.004)“CR3”組發生復發的風險降低52%(SHR,0.48,95%CI:0.15-1.49,P=0.206)。圖5為疾病進展階段的累積發生曲線。

圖5 疾病進展階段(phase)的累積發生曲線
競爭風險模型在Stata軟件中主要是采用“stcrreg”命令實現。分析前需要“stset”命令對生存數據中生存結局和生存時間進行定義。同時可以將使用的語句保存為Stata軟件DO文件,方便以后使用。值得注意的是“phase”變量是無序多分類變量,應以啞變量處理,“ib(first).phase”即以第一個組作為參考組。輸入“stcrreg,noshr”輸出回歸系數。關于競爭風險模型的“stcrreg”命令的其他選項,可以在Stata中輸入“help stcrreg”進行詳細學習。本文側重該方法的軟件實現,競爭風險模型的詳細理論請參考相應的文獻。
競爭風險模型在臨床研究和流行病調查中越來越受到研究者的認可和重視。相比需要編程基礎的SAS、R軟件,Stata軟件的菜單和命令可以更簡單、靈活地實現競爭風險模型的擬合。
[1]Fine JP,Gray RJ.A proportional hazards model for the subdistribution of a competing risk.Journal of the American statistical association,1999,94(446):496-509.
[2]Pintilie M.Competing Risks:A Practical Perspective.John Wiley&Sons:New York,2006.
[3]Scrucca L,Santucci A,Aversa F.Regression modeling of competing risk using R:an in depth guide for clinicians.Bone marrow transplantation,2010,45(9):1388-1395.
[4]Scheike TH,Maiers MJ,Rocha V,Zhang MJ.Competing risks with m issing covariates:effect of haplotypematch on hematopoietic cell transplant patients.Lifetime Data Anal,2013,19(1):19-32.
[5]Gray B.cmprsk:Subdistribution Analysis of Competing Risks.R package version 2.2-7.http://CRAN.R-project.org/package=cmprsk.
[6]陶莊.使用R軟件分析競爭風險模型簡明攻略.中國衛生統計,2008,25(6):638-639.
[7]Rosth?j S,Andersen PK,Abildstrom SZ.SASmacros for estimation of the cumulative incidence functions based on a Cox regression model for competing risks survival data.Computermethods and programs in biomedicine,2004,74(1):69-75.
[8]Lin G,So Y,Johnston G.Analyzing survival data with competing risks using SAS(r)software.SASGlobal Forum.
[9]“Proportional Subdistribution Hazards Model for Competing-Risks Data”,SAS Institute Inc.2013.SAS/STAT(r)13.1 User's Guide:pp5991-5995.Cary,NC:SAS Institute Inc.
[10]“stcrreg-Competing-risks regression”,StataCorp.2013.Stata 13 Base Reference Manual.College Station,TX:Stata Press.
(責任編輯:鄧 妍)
國家自然科學基金項目(81302516,81502885);教育部人文社會學研究項目(13YJCZH090)
1首都醫科大學公共衛生學院(100069)
2臨床流行病學北京市重點實驗室
#第一作者:李海彬、李霞(具有相同貢獻)
△通訊作者:郭秀花,E-mail:guoxiuh@ccmu.edu.cn