競(jìng)爭(zhēng)風(fēng)險(xiǎn)模型及其在Stata軟件實(shí)現(xiàn)*

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

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

競(jìng)爭(zhēng)風(fēng)險(xiǎn)模型及其在Stata軟件實(shí)現(xiàn)*

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

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

競(jìng)爭(zhēng)事件 競(jìng)爭(zhēng)風(fēng)險(xiǎn)模型 Stata軟件

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

模型原理和方法

1．競(jìng)爭(zhēng)事件

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

2．競(jìng)爭(zhēng)風(fēng)險(xiǎn)模型

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

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

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

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

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

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

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

3．競(jìng)爭(zhēng)風(fēng)險(xiǎn)模型軟件實(shí)現(xiàn)情況

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

實(shí)例分析

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

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

操作步驟

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

Stata軟件可以直接導(dǎo)入Excel文件，點(diǎn)擊菜單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合競(jìng)爭(zhēng)風(fēng)險(xiǎn)模型

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

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

結(jié)果與分析

1．結(jié)果

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

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

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

圖4 競(jìng)爭(zhēng)風(fēng)險(xiǎn)模型的結(jié)果

2．結(jié)果解釋

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

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

討論與結(jié)論

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

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

［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軟件分析競(jìng)爭(zhēng)風(fēng)險(xiǎn)模型簡(jiǎn)明攻略．中國(guó)衛(wèi)生統(tǒng)計(jì)，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é)任編輯：鄧 妍）

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

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

2臨床流行病學(xué)北京市重點(diǎn)實(shí)驗(yàn)室

＃第一作者：李海彬、李霞（具有相同貢獻(xiàn)）

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