微信直播

临床医生必修课:实例演示如何使用STATA 13.0做生存分析

Published at: 2015年第1卷第S1期

刘畅 , 贺永明
关键词:

 作者 | 刘畅,苏州大学14级医学研究生;贺永明,苏州大学附属第一医院心内科副教授。

1. 临床研究中存在大量的生存资料,相对于简单的“有或无”事件,生存资料中的事件与时间相联系,包含更多的信息,处理这种数据资料是每一个临床医生需要掌握的技能,其中最为我们熟知的即为Kaplan-Meier存活函数和COX回归比例风险模型。本文以简单实例逐步演示生存资料的分析过程,旨在帮助初学者以最快的速度掌握这一分析方法。

2. 生存分析概述:

生存分析(survival analysis)是既考虑结局,又考虑生存时间的一种分析方法,并可以充分利用不完全数据,对生存时间的分布特征进行统计描述和统计推断,也可以通过多因素模型对影响生存时间的主要因素进行分析。通过生存分析方法才能对临床随访资料进行全面和准确的评价。

1) 基本概念: 

生存时间(survival time):从起始事件到终点事件所经历的时间。

起始事件(initial event):是反映生存时间起始特征的事件,如疾病的确诊、治疗开始等。

终点事件(failure event):它是根据研究目的来确定的,一般是指反映治疗效果特征的事件,又称死亡事件或失效事件。

2) 生存时间资料的类型:

生存时间资料的类型一般分为完全数据和截尾数据,完全数据是指在整个随访研究期间能够观察到终点事件,即能够观察到从起点到终点的生存时间。截尾数据(censored data)指在随访过程中,由于某种原因未能观察到患者的明确结局(终点事件),又被称为删失数据。截尾的主要原因主要有三种:①失访;②退出;③研究终止。

3) 随访研究:

下图即为临床随访研究模式图:

随访前,首先需要明确以下事项:a)起始事件;b)终点事件;c)确定研究起点和研究终点,即明确研究的时间跨度,如2007年1月1日至2017年1月1日;d)记录影响个体终点事件发生的其他变量因素。

3.生存资料的数据分析:

获取生存资料后,下一步即开始数据分析,本文以Stata线上数据集“Patient Survival in Drug Trial”为例进行演示(Stata版本为13.0)。

1)首先打开该数据集:

命令:webuse drugtr

结果:

如图,Stata已载入线上数据“Patient Survival in Drug Trial”,并显示该数据名称。

2)通过list命令可将该数据集整体列出:

命令:list

结果:

如图,该数据集包括48名参与肿瘤药物试验的患者,年龄从47到67岁不等,其中28名患者给予药物治疗(drug=1),其余20名患者给予安慰剂治疗(drug=0)。变量studytime以月为单位记录患者治疗后的随访时间。变量died记录了随访过程中患者的随访结局状况,即终点事件死亡(died=1)或删失(died=0)。变量_st,_d,_t,_t0为Stata定义为生存资料后出现的变量。

3)生存资料的定义(该数据已被定义为生存资料):

命令stset用于帮助Stata识别数据资料中哪个变量用以记录随访时间,以及哪个变量表示观测个体是以完全数据或截尾数据结束随访的,即结局变量。    

命令:stset [时间变量] [结局变量]

对应本数据为:stset studytime died

结果:

其中studytime指随访时间变量,即发生终点事件或者删失时的时间减去纳入随访时的初始时间得到的月数。

died为结局变量,Stata视变量died不等于0的非缺失值为发生终点事件。

Stata会同时产生4个新的变量:

_st 代表:数据中该条记录是否被定义为生存资料。

_d 代表:数据中该条记录是否出现预期结果。

_t 代表:数据中观察对象被随访的时间。

_t0代表:数据中观察对象第一次被观察到的时间(开始的时间为0)。

4)生存资料的描述:

简要描述:

命令:stdescribe

结果:

相比list命令直接列出整个数据资料,stdescribe命令则为我们输出一个该数据资料的简要描述。如图,可以看出,该资料共有48个对象,且每个对象只有一条记录。此外还提供有进入和退出随访的时间、发生终点事件的例数等信息。

概要统计:

命令:stsum [if 表达式] ,[by(分组变量)  选择项] 

对应本数据:stsum,by(drug)

结果:

stdescribe命令的目的不是进行数据分析,只是对数据排列的描述。而stsum命令是对整个数据进行概要统计。如图,在744个人月中,共有31名个体发生终点事件,得到发生率为:0.0416667。得出不同drug组中25%、50%和75%生存时间,由于drug=1组中发生终点事件的个体数远小于其总体的75%,Stata对其75%生存时间无法估计,用缺失值表示

具体计算中位生存时间、平均生存时间、生存时间的百分数及其可信区间。命令:stci [if 表达式],[by(分组变量)  选择项]

其中选择项有:

median(中位生存时间);rmean(平均生存时间);P(#)(生存时间的百分数);level(#)(可信区间的可信度)

对应本数据:stci,by(drug) median

结果:

如图,可以得出不同drug组的50%生存时间,与stsum命令结果一致,但由于drug=1组中发生终点事件的例数接近50%,但并没有达到50%,Stata只能估计出其50%生存时间,但95%置信区间上限无法估计,stata以缺失值代替。

stci,by(drug) rmean

结果:

如图,drug=1组的平均生存时间大于drug=0组。对于观察队列中最后一例为删失者,未发生终点事件,平均生存时间的估计值偏低,Stata在相应数值后加“*”表示。

 stci,by(drug) p(25)

结果:

如图,计算得出不同drug组25%生存时间,与stsum命令结果一致,并得出25%生存时间的95%置信区间。

5)生存率的估计:

生存分析中,对于生存率的分析有非参数法、半参数法和参数法三种形式,临床中最常用的为非参数法的Kaplan-Meier存活函数和Nelson-Aalen累积风险函数,以及半参数法的COX回归比例风险模型,下面将对以上三种方法继续进行演示。  

Kaplan-Meier生存函数曲线:

命令:sts graph,[by(分组变量) 绘图命令选择项]

 其中主要选项有:

failure(指定绘制“死亡曲线”,与生存曲线相反)

gwood(绘制生存曲线的可信区间)

lost(在曲线上标出该时间点的删失值例数)

risktable(在曲线下增加X轴时间节点的风险例数)

对应本数据:

sts graph , by(drug) gwood lost 

结果:

如图,两组患者其生存率随随访时间逐渐下降,期中drug=0组中最后的个体在随访时间内以终点事件结束,即生存率从1随随访时间的延长下降到0,而drug=1组最后的个体以删失的形式结束,生存率最后未归于0,并且两条曲线没有交叉,可见drug=1组的生存率明显比drug=0组高。图中同时给出生存曲线的95%可信区间和各时间区间的的删失例数。

进一步以sts list命令可以具体得出曲线中每个随访时间对应的生存率以及标准误和95%置信区间。

命令:sts list,by(drug)

结果:

如图,经sts list命令可以算出,drug=1组的生存率,从1随随访时间下降到0.3476,而drug=0组的生存率由1下降到0。

Kaplan-Meier生存曲线曲线下常规增加一个风险例数表,即表示不同时间节点面临风险的个体数的表格:

命令:sts graph , by(drug) risktable

结果:(如图红色部分)

Nelson-Aalen累积风险函数曲线:

命令:sts graph,[by(分组变量) na 绘图命令选择项]

对应本数据:sts graph,by(drug) na

结果:

同样,结合sts list命令可以算出,drug=0组的累积风险,从0随随访时间增长到3.3989,而drug=1组的累积风险从0随随访时间的延长增长到0.9866,并且两条曲线没有交叉,drug=1组的死亡风险明显比drug=0组低。

对应sts list命令及结果如下:

命令:sts list,by(drug) na

结果:

以Log-rank检验进行生存率的比较:

以上是以图形的形式对生存率进行比较,进一步以统计计量检验两组或多组生存率是否相同,一般采用Log-rank检验。

命令:sts test [分组变量],[选择项]

主要选择项有:

logrank(进行Log-rank检验)

trend(检验生存率是否随分组变量取值水平的增高的变化趋势是否有意义)

对应本数据:

sts test drug,logrank

结果:

Logrank检验结果P值为:P<0.0001,按照a=005的检验水准认为两组病人的生存率不同。

Trend趋势检验要求分组变量有3组及3组以上才可运算,假设有该分组变量num3_variable。

命令:sts test num3_variable,trend

因本数据集中无3组及3组以上的分组变量,此处不再做演示。

Cox比例风险模型:(本资料主要研究变量drug对生存时间的影响,变量age作为混杂因素需要矫正)

在临床研究中,对因变量进行研究时通常需要考虑多个协变量的影响,但由于生存资料的特殊性,一般的方法难以解决,COX比例风险模型的提出,很好地解决了多因素生存分析的问题。

设有影响生存时间t的k个危险因素。设hi(t)为第i名受试者在时刻t的风险率即t时刻外后一瞬间的死亡速率。又设h0(t)表示不受危险因素x的影响下在时刻t的风险率,又称为基准风险率或基准函数。其模型的具体形式为:hi(t)/ h0(t)= exp(β1xi1+β2xi2+…+βmxim)

式中hi(t)/ h0(t)为相对于基准风险,受试者在t时刻后的相对瞬时死亡风险。X=(xi1,xi2,…,xim)'是可能与生存时间有关的m个危险因素所构成的变量。

计算各影响因素对生存时间的影响:

命令:stcox [协变量],[选择项]

对应本数据:

stcox i.drug age,hr

选项hr用于计算相对风险比

结果:

结果显示,Stata以drug=0组作为参照,对年龄进行校正后,drug=1组患者瞬时死亡的相对风险比为:0.105,95%置信区间为(0.043-0.256),p=0.000。即相对于使用安慰剂资料的患者,以药物治疗的患者有着更小的死亡风险,药物有效。对于连续变量age,其HR值表示年龄每改变一个单位所引起的瞬时死亡风险的改变量。即年龄增加1岁,其死亡风险是之前的1.12倍。P=0.002,具有统计学意义。

 stcox i.drug age,nohr

选项nohr 使Stata不再给出HR,而是给出回归系数 。

结果:

回归系数用于反映因素对生存时间的影响性质及强度,变量drug其系数为负值,也就是所变量drug是一种保护性因素,而变量age的回归系数为正值,即变量age是一种危险因素。一般而言,回归系数绝对值越大,则对应因素对生存时间的影响越大,但并不绝对,需要进一步比较其标准回归系数。

由于Stata中没有直接的求得COX回归标准回归系数的命令,需先对各变量进行标准化,再行回归计算求得。

变量标准化:

命令:

egen std_drug=std(drug)

egen std_age=std(age)

结果:

计算标准回归系数:

命令:stcox std_drug std_age,nohr

结果:

如图,变量drug的标准化系数绝对值为:1.12,变量age的标准化系数为:0.64。即变量drug对生存时间的影响比变量age大。

绘制不同分组的COX模型生存曲线:

默认得,Stata以所有变量的平均值来模拟生存曲线,我们也可以通过at选项指定一个变量为任意具体数值来进行模拟,以此获得分类变量不同分组的曲线。未指定的变量仍以其平均值参与模拟。并且通过at1() at2() at3()的形式可同时获得多条曲线。

命令:stcurve, survival at1(varname=#) at2(varname=&)

对应本数据:stcurve, survival at1(drug=0) at2(drug=1)

结果:

如图,我们可以分别得到drug=0组和drug=1组两组病人的模拟生存曲线,变量age以其平均值参与模拟过程。两条曲线不相交,即矫正混杂因素年龄后,不同的治疗处理对生存率的影响仍存在差异。

趋势检验

对于研究变量为三组及以上的等级变量,多需要对随变量等级变化的相对风险比进行趋势检验,假设存在该变量num3_variable,直接将该变量以连续变量而不是哑变量的形式纳入模型,矫正后的P值即为趋势检验的结果。

命令:stcox num3_variable co_variable1 co_variable2,hr

结果:

总体上COX比例风险模型与Logistics模型的分析相类似,其间的统计思路和方法可以相互借鉴。

至此,基本的生存分析步骤演示结束。临床研究中,生存分析至关重要,可以更好地帮助我们理解和解决问题,希望本文对你的临床研究有所帮助!

comments powered by Disqus

附件