李海彬 李 霞 王安心 陶麗新 劉 龍 陳斯鵬 郭 晉 郭秀花△

·計算機應(yīng)用·

競爭風(fēng)險模型及其在Stata軟件實現(xiàn)*

李海彬1，2＃李 霞1，2＃王安心1，2陶麗新1，2劉 龍1，2陳斯鵬1，2郭 晉1，2郭秀花1，2△

目的介紹競爭風(fēng)險模型及在Stata軟件中的實現(xiàn)。方法使用競爭風(fēng)險模型對白血病患者接受骨髓干細胞移植進行建模，展示該模型在Stata軟件中的實現(xiàn)過程并對最終結(jié)果進行解釋。結(jié)果Stata軟件“Competing-risks regression”菜單和“stcrreg”命令可以實現(xiàn)競爭風(fēng)險模型。實例應(yīng)用中，當(dāng)考慮移植相關(guān)死亡事件作為疾病復(fù)發(fā)的競爭事件后，基于競爭風(fēng)險模型得到病人的疾病進展階段與復(fù)發(fā)有關(guān)。結(jié)論競爭風(fēng)險事件廣泛存在于臨床研究中，Stata軟件可以簡單靈活地實現(xiàn)這一過程。

競爭事件 競爭風(fēng)險模型 Stata軟件

經(jīng)典生存分析中，研究者往往只關(guān)注研究對象隨訪過程中特定事件的發(fā)生。然而在現(xiàn)實世界中，研究對象不僅經(jīng)歷一種類型事件，不同類型結(jié)局事件相互影響，即形成競爭事件。處理含有競爭事件數(shù)據(jù)的統(tǒng)計模型稱作“競爭分析模型”。1999年Fine和Gray提出部分分布的半?yún)?shù)比例風(fēng)險模型，該模型在考慮到競爭事件的條件下，探討其他影響因素對結(jié)局事件的作用［1］，彌補了傳統(tǒng)生存分析的不足。目前競爭風(fēng)險模型被廣泛應(yīng)用于臨床試驗、流行病學(xué)調(diào)查等研究中。然而國內(nèi)文獻中介紹競爭風(fēng)險模型及Stata程序?qū)崿F(xiàn)相對較少，本文將結(jié)合臨床實例介紹競爭風(fēng)險模型及如何使用Stata軟件進行建模，旨在更好地為研究中提高對競爭事件的認識和通過Stata軟件實現(xiàn)這一過程操作提供幫助。

模型原理和方法

1．競爭事件

競爭風(fēng)險（competing risks）是指研究對象出現(xiàn)感興趣事件（interesting event）的同時還會出現(xiàn)其他結(jié)局事件，這些結(jié)局事件將阻止感興趣事件的出現(xiàn)或使其發(fā)生的概率降低，各結(jié)局事件間形成所謂的＂競爭＂關(guān)系，這一系列事件稱作競爭事件（competing risk events）［2］。例如在白血病患者接受骨髓干細胞移植后復(fù)發(fā)風(fēng)險的研究中，應(yīng)當(dāng)將移植相關(guān)死亡事件作為競爭事件處理［3］，倘若忽視了該競爭事件的存在，而簡單的將競爭事件等同于右刪失（right censoring），使用傳統(tǒng)K-M法估計終點事件發(fā)生概率、Cox模型進行多因素分析，將產(chǎn)生偏差［4］。

2．競爭風(fēng)險模型

1999年Fine和Gray［1］提出的部分分布的半?yún)?shù)比例風(fēng)險模型（Fine-Gray模型）就是一種競爭風(fēng)險模型（competing risksmodel），使用累積風(fēng)險函數(shù)（cumulative incidence function，CIF）來估計結(jié)局事件的累積發(fā)生概率。在該模型中，t時刻發(fā)生事件j風(fēng)險定義如下［1］：

其中λ0（t）是事件j的部分分布基準風(fēng)險（baseline hazard of the sub-distribution）；Z是協(xié)變量；βj′Z協(xié)變量的回歸系數(shù)；exp｛βj′Z｝是部分分布相對風(fēng)險（sub-distribution hazard ratio）。部分分布風(fēng)險模型偏似然估計定義如下：

風(fēng)險集Rj定義如下：

風(fēng)險集是由在時間t沒有發(fā)生感興趣事件個體和在時間前發(fā)生競爭事件的個體構(gòu)成。因此經(jīng)歷其他類型事件的個體仍在風(fēng)險集中。權(quán)重定義如下：

G（t）是Kaplan-Meier方法估計的生存函數(shù)。時間t前沒有發(fā)生競爭事件的個體，在時間t發(fā)生感興趣事件有相等的權(quán)重（wi＝1）；ti＜t時，發(fā)生競爭事件個體權(quán)重wi＜1。

競爭風(fēng)險模型通過以下公式建立λj（t，Z）與CIF的關(guān)系：

因此，通過競爭風(fēng)險模型能夠估計出協(xié)變量Z相應(yīng)的系數(shù)β′j，聯(lián)合累積部分分布基準風(fēng)險預(yù)測特定時間個體發(fā)生感興趣事件發(fā)生概率，即絕對風(fēng)險。

3．競爭風(fēng)險模型軟件實現(xiàn)情況

競爭風(fēng)險模型可以通過R軟件“cmprsk”程序包實現(xiàn)［5］，可參考文獻［3、6］［3，6］，或SAS軟件宏“%Cum Inc”［7］、“%CIF”［8］實現(xiàn)以及SAS9．4版本“PRCO PHREG”［9］。考慮到使用上述過程需要有一定的編程基礎(chǔ)，而Stata軟件可以使用“stcrreg”命令擬合競爭風(fēng)險模型［9］，也可以簡單通過Stata軟件菜單操作完成，極大的方便了臨床醫(yī)生及其他研究者對競爭風(fēng)險模型的實現(xiàn)。本文結(jié)合Stata 14．0軟件菜單操作和命令語句結(jié)合臨床實例擬合競爭分析模型，為更好地通過Stata軟件實現(xiàn)競爭風(fēng)險模型分析提供幫助。

實例分析

本文的實例數(shù)據(jù)是關(guān)于白血病患者接受骨髓干細胞移植后復(fù)發(fā)風(fēng)險的研究［3］，下載地址http：//www．stat．unipg．it/luca/R/，在Competing Risks Analysis目錄，下載bmtcrr．xls文件。研究的主要目的是考慮移植相關(guān)死亡即競爭事件后，探索病人復(fù)發(fā)的累積發(fā)生率以及復(fù)發(fā)的影響因素。該研究包含有177例接受骨髓干細胞移植的急性白血病患者。表1為變量命名、賦值及研究對象的基本特征。

表1 變量命名、賦值及統(tǒng)計描述

操作步驟

1．?dāng)?shù)據(jù)集導(dǎo)入與變量命名、賦值

Stata軟件可以直接導(dǎo)入Excel文件，點擊菜單File—Import—Excel spreadsheet（*．xls*．xlsx）—Browse—選擇bm tcrr．xls文件—勾選Import first row as variable names—OK。數(shù)據(jù)導(dǎo)入結(jié)果如圖1所示。原始數(shù)據(jù)是字符型，需要將變量賦值。命令如下：

圖1 數(shù)據(jù)導(dǎo)入及原始數(shù)據(jù)結(jié)構(gòu)

2．?dāng)M合競爭風(fēng)險模型

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

圖2 生存數(shù)據(jù)定義及競爭風(fēng)險模型菜單操作

結(jié)果與分析

1．結(jié)果

圖3為生存數(shù)據(jù)定義結(jié)果，顯示177名病人中56人發(fā)生復(fù)發(fā)事件，共觀察到299．08人年（3588．94/12），最長隨訪時間為131．77個月。

圖3 生存數(shù)據(jù)定義結(jié)果

圖4 為競爭風(fēng)險模型的結(jié)果。右側(cè)“No．of obs、No．of subjects、No．failed、No．competing、No．censored＂分別是對觀察個體、研究對象、出現(xiàn)結(jié)局事件、競爭事件、刪失的個數(shù)的統(tǒng)計?！癓og pseudolikelihood＝－266．52444，Wald chi2（7）＝22．93，Prob＞chi2＝0．0018＂是對競爭風(fēng)險模型整體評價，P＜0．05，具有統(tǒng)計學(xué)意義。圖4中的下半部分是多因素回歸結(jié)果。從左列到右列分別是變量、部分分布風(fēng)險比（SHR）、標準誤、統(tǒng)計量（Z）、P值和SHR的95%置信區(qū)間（95%CI）。

圖4 競爭風(fēng)險模型的結(jié)果

2．結(jié)果解釋

以性別變量為例，在考慮移植相關(guān)死亡為競爭風(fēng)險事件的前提下，控制年齡、疾病類型、進展階段、移植類型后，女性白血病患者骨髓移植后發(fā)生復(fù)發(fā)的風(fēng)險是男性的0．97倍（95%CI：0．55－1．71）。95%CI包含1，同時P＝0．904＞0．05，差異不具有統(tǒng)計學(xué)意義。其它變量的解釋與性別變量解釋類似，但疾病進展階段變量需要注意，為無序變量，采用啞變量處理。以狀態(tài)“Relapse”為參考組，考慮競爭事件后“CR1”組發(fā)生復(fù)發(fā)的風(fēng)險降低67%（SHR，0．33，95%CI：0．16－0．70，P＝0．004）；“CR2”組發(fā)生復(fù)發(fā)的風(fēng)險降低64%（SHR，0．36，95%CI：0．18－0．73，P＝0．004）“CR3”組發(fā)生復(fù)發(fā)的風(fēng)險降低52%（SHR，0．48，95%CI：0．15－1．49，P＝0．206）。圖5為疾病進展階段的累積發(fā)生曲線。

圖5 疾病進展階段（phase）的累積發(fā)生曲線

討論與結(jié)論

競爭風(fēng)險模型在Stata軟件中主要是采用“stcrreg”命令實現(xiàn)。分析前需要“stset”命令對生存數(shù)據(jù)中生存結(jié)局和生存時間進行定義。同時可以將使用的語句保存為Stata軟件DO文件，方便以后使用。值得注意的是“phase”變量是無序多分類變量，應(yīng)以啞變量處理，“ib（first）．phase”即以第一個組作為參考組。輸入“stcrreg，noshr”輸出回歸系數(shù)。關(guān)于競爭風(fēng)險模型的“stcrreg”命令的其他選項，可以在Stata中輸入“help stcrreg”進行詳細學(xué)習(xí)。本文側(cè)重該方法的軟件實現(xiàn)，競爭風(fēng)險模型的詳細理論請參考相應(yīng)的文獻。

競爭風(fēng)險模型在臨床研究和流行病調(diào)查中越來越受到研究者的認可和重視。相比需要編程基礎(chǔ)的SAS、R軟件，Stata軟件的菜單和命令可以更簡單、靈活地實現(xiàn)競爭風(fēng)險模型的擬合。

［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軟件分析競爭風(fēng)險模型簡明攻略．中國衛(wèi)生統(tǒng)計，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．

（責(zé)任編輯：鄧 妍）

國家自然科學(xué)基金項目（81302516，81502885）；教育部人文社會學(xué)研究項目（13YJCZH090）

1首都醫(yī)科大學(xué)公共衛(wèi)生學(xué)院（100069）

2臨床流行病學(xué)北京市重點實驗室

＃第一作者：李海彬、李霞（具有相同貢獻）

△通訊作者：郭秀花，E-mail：guoxiuh＠ccmu．edu．cn