本站分享:AI、大数据、数据分析师培训认证考试,包括:Python培训Excel培训Matlab培训SPSS培训SAS培训R语言培训Hadoop培训Amos培训Stata培训Eviews培训

stata生存分析的详细步骤及解释_stata 生存分析

stata培训 cdadata 21733℃

stata生存分析的详细步骤及解释

关键词: stata 生存分析,stata 生存曲线,stata做生存分析

前言

生存分析(survive analysis)即是将终点事件的出现与否和出现终点事件所经历的时间结合起来的一种统计分析方法,生存分析通常研究的终点事件是死亡,生存分析由此得名。但生存分析可更广泛的运用于恶性肿瘤、慢性疾病或其他情况的随访研究中事件分析,比如疾病的发生、复发、转移、伤口的愈合、某种症状的消失等。生存资料的分析主要特点就是考虑每个研究对象出现某一结局所经历的时间。生存曲线即是以生存时间为横轴,生存率为纵轴,将各个时间点对应的生存率连接在一起的曲线图[1-2]。

1. 生存分析中几个重要的基本概念

1.1生存时间

生存时间(survival time)也是一个广义概念,泛指所关心的某现象的持续时间,即随访观察持续的时间,常用符号t表示。生存时间分为两种类型:1. 完全数据(complete data):指从观察起点到发生“死亡”事件所经历的时间。提供了观察对象确切的生存时间。2. 截尾数据(censored data):亦称截尾值(censored value)或终检值。指从观察起点到发生非“死亡”事件所经历的时间。

1.2截尾

生存结局分为“死亡”与“截尾”两类,“死亡”是感兴趣的终点时间,其他终点事件或结局都归为截尾。

1.3死亡概率

死亡概率(probability of death)表示单位时间段开始存活的个体,在该段时间内死亡的可能性。符号q表示。”q=某年内死亡人数÷某年年初人口数”

1.4生存概率

生存概率(probability of survival)表示单位时间段开始存活的个体,到该段时间结束时仍存活的可能性。符号p表示。p=某年活满一年人口数÷某年年初人口数。P=1-q。

1.5生存率

生存率(survival rate, survival function)表示观察对象经历tk个单位时间段后仍存活的可能性。若无截尾数据,则。其中1。若有截尾数据,须分时段计算生存概率。假定观察对象在各个时段的生存事件独立,应用概率乘法定理:,Pi为某时段的生存概率,故生存率又称累积生存概率(cumulative probability of survival)。

1.6生存曲线

生存曲线(survival curve):生存时间为横轴,将各时点所对应的生存率连接在一起的曲线图,样本量小时生存曲线呈阶梯形,样本量足够大时,形成光滑的曲线。

1.7中位生存时间

中位生存时间是指50%观察对象能存活的时间。

2. 生存分析的统计学方法

由于生存时间一般不呈正态分布,而且需要考虑截尾数据,生存分析有其独特的统计学方法。常用的统计学方法有以下几种。

2.1描述性分析

根据样本生存资料估计总体生存率及其他有关指标(如中位生存时间等)。常采用Kaplan-Meier法(乘积极限法)进行分析。对于频数表资料则采用寿命表法进行分析。计算生存率需要考虑时间顺序。

2.2比较分析的方法

对不同组生存率进行比较分析,常采用非参数的log-rank检验,检验无效假设使两组或多组总体生存时间分布相同。

2.3影响因素分析

通过生存分析模型来探讨影响生存时间的因素,常用的方法为COX比例风险模型。

3. 基于Stata软件的统计学实现生存分析(笔者注:以下所举实例数据全部来自于陈峰教授主编《现代医学统计方法与Stata应用(第2版)》,相关Stata命令及结果解释大部分来自于这本书,其中部分命令有少许改动。陈锋教授主编的这本书通俗易懂,感兴趣的读者可以找来一读)

3.1生存资料的定义

在对随防资料进行生存分析之前,需先将该数据库定义为生存资料数据库,其命令是:

stset 时间变量 [,failure(截尾变量==#)]

其中,选择项failure(截尾变量==#)规定截尾变量取值为“#”时研究对象出现预期结果,没有该选择项时,Stata 以所有不等于0 的非缺失值为出现预期结果。对数据库进行定义时必须注意变量顺序,命令stset 后的变量顺序依次为时间变量、截尾变量。定义数据库后,系统自动产生四个变量:

_st /* 数据库中该条记录是否被定义为生存资料

_d /* 该条记录是否出现预期结果

_t /* 观察对象被随访的时间

_t0 /* 观察对象第一次被观察到的时间(开始过程的时间为0) 

例1 某医院泌尿外科于1979-1982 年间作了19例肾移植手术,拟了解肾移植后病人

的生存时间(天)。规定随访开始时间为病人术后一天,预期结果为该病人因与肾移植有关的各种原因的死亡。后改进手术方式,于1983-1986 年又作了14 例,资料如下(有+的数据表示该病人截尾)。计算各组的生存率及可信区间(资料已存入文件1.dta)。

一般手术: 3 9 15 20 20 26 30 41 46 64+ 64 135 223 365 450 596+ 680+ 900+ 900+

改进手术: 10 390+ 518+ 70+ 70+ 120 475+ 225 801+ 366 647+ 1001+ 1045+ 1045+

stata生存分析的详细步骤及解释_stata 生存分析

键入命令如下:

. stset time, failure(outcome)

得结果如下:

stata生存分析的详细步骤及解释_stata 生存分析

数据库“例-1”被定义为生存分析数据库,变量“outcome”取值不等于0 且不等于缺失值时,该记录为完全数据,即出现预期结果。反之则为截尾值,表示未观测到病人出现预期结果。完成上述定义后,即可用下面介绍的命令作进一步分析。

3.2 生存资料的描述

用于计算中位生存时间的命令是:

stsum [if 表达式] [,by(分组变量)]

可用stci 命令计算中位生存时间、平均生存时间、生存时间的百分位数,及其可信区间:

stci [if 表达式], [by(分组变量) 选择项]

其中,选择项有:

median /* 计算中位生存时间

emean /* 计算平均生存时间时,如果生存时间最长一例为截尾值,emean 假设数据服从指数分

,并根据指数分布将该例后生存曲线部分延长至与横轴相交,曲线下面积即为所求

的平均生存时间。

rmean /* 计算平均生存时间时,如果生存时间最长一例为截尾值,rmean 不对数据延长,曲线

下面积即为所求的平均生存时间。此即为通常教科书上所教授的平均生存时间。

p(#) /* 生存时间的百分位数

level(#) /* 可信区间的可信度

也可用survsum 命令计算中位生存时间的中位数。

继续以例1数据为例,在命令窗口键入:

. stsum, by(treat)

数字化结果如下:

stata生存分析的详细步骤及解释_stata 生存分析

第二组(改进手术组) 较早出现了截尾数据,故该组的中位生存时间无法进行估计,Stata

用缺失值表示。

用stci 命令可以计算平均生存时间及其可信区间:

命令窗口键入命令如下:

. stci, rmean by(treat)

数字化结果如下:

stata生存分析的详细步骤及解释_stata 生存分析

第二组的平均生存时间明显长于第一组。对于观察队列中最后一例为截尾值者,平均生

存时间的估计值偏低(underestimated)。Stata 在相应数值后加“*”表示。

3.3生存率的估计

生存率的估计一般采用乘积极限(product-limit)法,又称Kaplan-Meier 法,其标准误的计

算用Greenwood 近似法。根据生存率及其标准误,可以绘制生存曲线,估计可信区间。

用于输出生存率、生存率的标准误等统计量的命令是:

sts list [if 表达式] [, by(分组变量) strata(分层变量) adjustfor(校正变量选择项]

选择项有:

failure /* 输出死亡函数 (1-S(t))

na /* 输出累积风险函数

level /* 规定所输出可信区间的可信度

这里,by 与strata 选择项的使用有所不同。使用by 选择项时,Stata 对分组变量的不同水平分别计算生存函数和累积风险函数。而在使用strata 选择项时必须同时使用adjustfor 选择项,此时Stata 将计算adjustfor 选择项中校正变量取值为0 时的生存函数、累积风险函数,即计算基线生存函数、基线累积风险函数。

用于绘制Kaplan-Meier 生存(死亡)曲线的命令是:

stsgraph [,by(分组变量) separate 绘图命令选择项/* 绘制Kaplan-Meier 生存(死亡)曲线

stphplot [,by(分组变量绘图命令选择项/* 绘制log(-log(S(t)))log(time)的线图

stcoxkm [,by(分组变量) separate 绘图命令选择项/* 绘制Kaplan-Meier 生存曲线与Cox 预测曲线

在上述的“绘图命令选择项”中,可以选用xlab,ylab,xtick,ytick 以及标题命令t1,t2,b1,

b2,l1,l2,r1,r2 等,但connect, symbol, pen 等选择项不能用,因y 轴是从0 开始,所以ylog 等选择项亦不能用。在stphplot 命令中,可选用connect, symbol, pen 等选择项。

sts graph 命令中的其他常用选择项:

adjustfor(变量) /* 按指定变量进行调整

failure /* 指定绘制死亡曲线,缺失为绘制生存曲线

gwood /* 绘制生存(或死亡)曲线及其Greenwood 可信区间

na /* 绘制Nelson-Aalen 累积风险函数曲线

cna /* 绘制Nelson-Aalen 累积风险函数曲线及其可信区间

lost /* 在曲线上标出该时间点截尾值例数

计算各组的生存率及标准误,命令及结果如下:

. sts list if treat==1

stata生存分析的详细步骤及解释_stata 生存分析

绘制各组的生存曲线,命令及结果如下:

. sts graph,by(treat) lost

stata生存分析的详细步骤及解释_stata 生存分析

两条曲线分别表示两组的生存曲线,曲线上的数字表示在该时刻的截尾值例数。显然,两

组的生存率不同。绘制各组的生存曲线及其可信区间,使用gwood 选择项。如对第1 组,命令及结果如下:

. sts graph if treat==1, gwood

stata生存分析的详细步骤及解释_stata 生存分析

图中,中间一条线是treat=1 组的生存曲线,上、下两条线分别表示生存率的可信区

间的上下限。注意,率的可信区间是不对称的。

3.4 生存率的比较

一、两组或多组生存率的比较

检验两组或多组生存率是否相同一般采用Log-rank(Mantel–Cox)检验、Wilcoxon(Breslow)

检验。命令如下:

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

选择项有:

logrank /* 进行Log-Rank 检验

wilcoxon /* 进行Wilcoxon(也叫Wilcoxon-Gehan 检验或Breslow 检验)

trend /* 检验死亡(生存)率是否随分组变量取值水平的增高而上升或下降

就例1资料,比较两组病人的生存时间有无差别。键入命令如下:

. sts test treat, logrank

数字化结果如下:

stata生存分析的详细步骤及解释_stata 生存分析

这里的检验假设是第一处理组的生存率与第二组的相同。输出结果中给出了两组的实际

数(Events observed)及理论数(Events expected)。本例中改进手术组的实际实际死亡数小于理论数,说明该组病人预后情况较好,经Log-rank检验,χ2= 6.71,自由度υ=1,P=0.0096,按α=0.05的检验水准认为两组病人的生存时间有差别,以改进手术组为优。

3.5 Cox 比例风险模型

恶性肿瘤患者生存时间的长短,不仅与治疗有关,还受病人的年龄、性别、病情、心理、环境、社会等因素的影响,如果要确切地显示治疗措施的效果,所有的病人除了治疗措施不同以外,其他影响因素必须相同(或相近),但这在实际上是不可能做到的。因此,我们最好能采用多因素分析方法,即分析包括治疗措施在内的可能因素对生存时间长短的影响(大小和方向)。

但生存时间的分布往往不服从正态分布(大多为正偏态分布),有时不知道它的分布类型,又存在截尾数据(Censored data)这样,就不能用多元线性回归方法来分析。而传统的方法只能进行单因素分析,又不能利用截尾数据(Censored data)。1972年,英国统计学家 D. R. COX提出了一种比例风险模型(Cox proportional hazard model),简称COX模型。它可以分析多种因素对生存时间的影响,而且允许有“截尾”存在。是生存分析中最重要的模型之一。COX模型主要用于肿瘤和其它慢性病的预后因素分析,也可以用于一般的临床疗效评价和队列的病因探索。Cox 比例风险模型的一般形式是:

stata生存分析的详细步骤及解释_stata 生存分析

h(t)表示研究对象在时点t暴露于协变量(x1, x2, …, xm)之下的风险函数,h0(t)表示所有协变量取值均为0 时的基线风险函数。在Cox模型中h0(t)不能由样本得出,因而不能估计生存率。但这并不妨碍对各协变量相对危险度的估计。

估计Cox比例风险模型的命令格式为:

stcox [协变量] [,nohr strata(分层变量)]

估计含有时依变量的Cox比例风险模型的命令格式为:

cox 时间变量 [协变量] [,hr dead(截尾变量) strata(分层变量选择项]

进行逐步Cox 回归分析的命令为:

sw cox 时间变量协变量 [, cox 命令选择项逐步回归选择项]

这里选择项有:

nohr /* 指定输出回归系数而不是危险比exp(b)

tvc(时依变量) /* 指定时依变量

texp() /* 时依变量取值变化表达式

level(#) /* 可信区间的可信度

maximize-options /* 进行最大似然估计的控制选项

[应用命令cox时无须事先应用stset 对数据进行定义,且进行逐步回归时只能使用cox 命令。

用 sw cox 命令可以进行逐步Cox 回归分析。

就例1资料进行Cox 回归分析。

在应用stset 对数据进行规定后,可直接用stcox 命令进行Cox 回归分析。键入命令如下:

. stcox treat,nohr

或者也可以使用如下命令:

. cox time treat, dead(outcome)

结果如下:

stata生存分析的详细步骤及解释_stata 生存分析

stata生存分析的详细步骤及解释_stata 生存分析

风险函数一般用极大似然估计,用Newton-Raphson 法迭代。结果中给出了每次迭代的似

然函数之对数值(Log Likelihood), 本例经四次迭代得极大似然估计变量treat 的系数

(Coefficient)为-1.371774,其标准误(Std.Error)为.5708971,z 值为-2.40,是对treat 的系数是否为0的检验,P>|z|是相应的概率。本例的P=0.0016,故可认为treat的系数不为0。从而得Cox 比例风险函数:

如果计算HR则可使用如下命令:

. stcox treat, hr

或者使用以下命令:

. cox time treat, dead(outcome) hr

结果如下:

stata生存分析的详细步骤及解释_stata 生存分析

stata生存分析的详细步骤及解释_stata 生存分析

以例2数据为例继续演示Stata软件实现Cox回归

某临床试验比较A,B 两治疗方案对某病的治疗效果,A 组(group=0)12 人,B

组(group=1)13 人。病人分组后检验其肾功能(kidney),功能正常者记0,不正常者记为1;治疗后生存时间为stime(天);数据已存入文件2.dta。问不同治疗方案及肾功能对病人的生存时间是否有影响?

这里,时间变量是stime,终检变量是censor,治疗方案(group)是研究因素,而肾功能

(kidney)是混杂因素。例2数据如下图所示:

stata生存分析的详细步骤及解释_stata 生存分析

键盘键入命令设置数据为生存数据,如下:

. stset stime, failure(censor)

结果如下:

stata生存分析的详细步骤及解释_stata 生存分析

继续键入命令如下:

. stcox group kidney, nohr nolog

得结果如下:

stata生存分析的详细步骤及解释_stata 生存分析

计算HR,则输入如下命令:

. stcox group kidney, hr

stata生存分析的详细步骤及解释_stata 生存分析

3.6 随访生存资料的寿命表法

当样本含量较大或不能准确得知研究结果出现的时间时,可以将各研究对象的生存时间

按年或月进行分组计算其生存率。Stata相应的命令是:

ltable 时间变量 [dead(截尾变量)] [,by(分组变量选择项绘图选择项]

ltable命令中大部分选择项前面已经介绍过,未介绍过的有:

weight /* 指定权重变量

survival /* 指定输出生存率

failure /* 指定输出死亡率

hazard /* 指定输出风险函数

interval /* 指定寿命表中生存时间组距

test /* 应用似然比检验、Log-rank检验对各总体生存率曲线是否相同进行检验

notab /* 不输出寿命表

noconf /* 绘制生存率曲线时不绘制各时间点生存率的可信区间

graph /* 指定绘制生存率曲线

例3 随访某种恶性肿瘤患者生存情况如下图所示,试作统计分析。这是一个分组资料,先将数据整理成下列形式,包括处理变量treat,生存年数year,是否截尾censor,以及频数num。其中,生存年数输入时“0~”输为0.5,“1~”输为1.5,其他依此类推。

stata生存分析的详细步骤及解释_stata 生存分析

计算寿命表并,进行统计学检验,命令如下:

. ltable year censor [weight=num], test by(treat) interval(1)

得数字化结果如下:

stata生存分析的详细步骤及解释_stata 生存分析

Stata依次输出各段生存时间起点及终点、期初人数、期内死亡人数、截尾例数、生存率

及其标准误和相应的95%可信区间。同时给出了两组的齐性检验(Lawlsee,1982)及log-rank检

验。

绘制第一组(group=1)病人的生存率曲线图。命令如下:

. ltable year censor [weight=num], test by(treat) interval(1) graph

得数字化结果如下:

stata生存分析的详细步骤及解释_stata 生存分析

4. 结语

生存分析应用广泛,作为一个临床医生至少应该掌握使用一种统计学软件实现生存分析,本文在参考了《现代医学统计方法与Stata应用(第2版)》基础上给大家演示了Stata软件实现生存分析的过程,希望能对大家的科研工作有所帮助。

转载请注明:数据分析 » stata生存分析的详细步骤及解释_stata 生存分析

喜欢 (10)or分享 (0)