1. 为什么说Stata做生存分析可以“极简”
干临床科研这几年,Stata给我留下的最深印象是:手里有数据、心里有分析思路,就能靠那几行命令把生存分析快速跑完。谈不上花哨,但胜在极简。和R语言需要写一堆tidyverse管道、SPSS需要在菜单里来回点选相比,Stata的命令行方式特别适合随访数据的处理节奏——你不需要中断思路去翻菜单,只要记住几个核心动词,就能把数据清洗、描述统计、KM曲线、Cox回归和亚组分析串成一条流水线。
很多人一听“生存分析”就觉得门槛高,其实拆开看就两件事:第一,描述每个时间点的存活情况;第二,找出影响存活的因素。Stata恰好把这两件事都压缩成了几个命令。从stset声明数据结构开始,到sts list输出生存表,再到stcox跑Cox模型,正常情况下半小时内就能完成整套分析。如果你正在做随访队列、临床试验数据,或者毕业论文需要处理“复发时间”“死亡时间”“并发症发生时间”这类数据,这套极简流程可以直接抄作业。
这篇文章我不打算铺开讲生存分析的理论公式,而是聚焦在“怎么用Stata一步步跑通生存分析”这件事上。后面所有内容都依托Stata实操,适合医学生、公卫方向的研究者,以及想快速上手生存分析但不想啃大部头教材的社科研究者。学过一点Stata基础操作就够,没学过也没关系,跟着每一步走就行。
2. 生存分析必会的基础概念与数据结构
2.1 时间变量与事件变量:数据的两个基本盘
任何生存分析数据集,核心都离不开两个变量:生存时间(time)和结局事件(status/event)。
生存时间的定义看似简单,实际容易出错。随访研究中,每个人从某个起点开始被观察,一直到出现目标事件、失访、或者研究结束。我们记录的是“从起点到终点的时间长度”。这个终点可能是术后复发、死亡、出院、或某种并发症发生。起点是手术日期、入组日期、治疗开始日期,需要一开始就明确。
我见过不少数据在时间变量上翻车的情况。比如有人把“随访日期”当成时间变量,还有人把“确诊年份”直接拿来算生存时间,却没有减去起点的日期。Stata本身不负责判断你的时间计算方式对不对,它只按你指定的变量去跑,所以时间变量一定要在数据整理阶段就计算清楚。稳妥的做法是在原始数据里保留三个日期:起点日期、终点日期、以及一个状态变量(0=未发生事件/截尾,1=发生目标事件)。用双日期相减生成生存时间,而不是手工填写一个“随访月数”。
2.2 截尾数据:不能删,也别当普通观察
生存分析里最特别的概念是截尾(censoring)。简单说,一个人到研究结束时还没发生事件、中途失访、或者死于其他原因(和目标事件无关),那他的实际事件时间是不完整的,这类数据就叫截尾数据。
很多初学者会犯一个错误:直接把这些截尾病例删掉。这绝对不行。截尾数据虽然不知道确切的事件发生时间,但至少知道“他撑过了这段时间没出事”,这个信息本身就是有价值的。比如一个病人在术后随访到第10个月时失访,我们至少知道他在前10个月没有复发,这对估计生存函数是有贡献的。全部删掉会导致生存率被严重低估。
Stata在处理截尾时用的是状态变量编码习惯:通常用0表示截尾,1表示发生事件。这个编码不是唯一标准,但绝大多数命令默认按这个逻辑理解。所以清洗数据时,建议统一把事件变量的取值设为0和1,并加上值标签,这样跑出来的结果和图表标注都清晰。
2.3 stset:生存分析的第一步,声明数据结构
在Stata里做任何生存分析,第一件事就是告诉Stata:“我们这个数据集里的时间变量是哪个,事件变量是哪个”。这个动作就是stset。
用法非常直接:
stset time, failure(status==1)意思是:设定生存时间变量time,并定义事件发生条件为status等于1。这条命令执行后,Stata会自动建立一个名为_st的变量记录观察时间、_d记录事件状态,并在输出窗口告诉你总样本量、事件数、截尾数、观察人时等信息。这些都是后续所有生存分析命令的基础。
stset还有一些常用选项。如果数据是带时间区间的计数类型,比如随访是区间数据,可以用id()标识个体;如果时间是以日期存储的,还可以直接做日期转换后计算。但极简操作用法里,最核心的就是上面这一行。
值得注意的是,stset可以反复执行。比如你换了子样本分析,或改变了事件定义,重新跑一次stset并不会产生破坏性,Stata会用新的设定覆盖旧的。这一点很实用,亚组分析时不需要反复创建临时数据集。
3. 手把手跑通Kaplan-Meier分析与生存曲线
3.1 sts list:图表之前先看数字
生存分析的描述性阶段,最常用的就是Kaplan-Meier方法,简称KM分析。它能计算每个时间点的累积生存概率,并且能自动处理截尾数据。
极简流程的第一步是sts list。这个命令输出一张生存表,包含每个时间点的生存人数、事件数、失访数、生存概率和标准误。刚开始做分析时,我会强烈建议你先跑sts list,而不是直接画图。原因很简单:曲线图容易让人产生“很直观”的错觉,但数字才是判断数据是否合理的基础。
看表的时候重点看几个信息:总人数有多少、事件数是多少、最后一个时间点还剩下多少人。如果最后只剩几个人,那生存曲线尾巴部分其实非常不稳定,报告时需要谨慎。比如一个500人的队列,随访到60个月时只剩下5个人,那60个月的生存概率就算算出来了,也没有太大参考价值。
stset time, failure(status==1) sts list输出结果里关注N_subjects、N_failed、N_censored这几列,配合生存概率和置信区间,基本就能对总体生存情况有一个准确判断。
3.2 sts graph:出图不等于做完了分析
数字看过了,接下来用sts graph出KM曲线。
sts graph, survival这一行命令直接出默认的生存曲线图。如果你希望按某个分组变量分别画曲线,比如按治疗组和对照组分组,可以加上by()选项:
sts graph, by(treatment) failure这里如果用failure选项,画的就是累积风险曲线,也就是1减生存概率的曲线,在期刊里也很常见。按需选择即可。
出图之后别急着截图保存。先检查曲线形态是否合理。几条线有没有交叉?曲线下降的坡度是否和临床预期一致?置信区间的宽带是否过宽?这些问题是审稿人关注的焦点,也是判断整个分析是否可靠的直观线索。
图形的美化可以放在最后做,用graph export导出高分辨率图片。Stata默认图形颜色和字体可能不合期刊要求,但修改起来也不难:
sts graph, by(treatment) name(km, replace) /// title("Kaplan-Meier 生存曲线") ytitle("生存概率") graph export "km.png", width(2400) replace导出时设置宽度2400像素,基本能满足多数期刊的300dpi要求。
3.3 sts test:组间差异到底有没有意义
图形上两条KM曲线看起来分开了,但到底是抽样误差还是真实差异,就要用检验统计量来判断。Stata里最简单的组间比较命令是sts test。
sts test treatment默认状态下sts test跑的是log-rank检验,这也是生存分析文献里最常用的组间比较方法。如果生存曲线的组间差异在早期比较明显、后期逐渐消失,log-rank检验的权重均等可能不太敏感,这时可以考虑用wilcoxon选项给早期事件更高的权重:
sts test treatment, wilcoxon一般情况下,log-rank检验已经足够应付大多数临床研究会。结果看Probability值(p值),以0.05为界判断差异是否显著。
需要提醒的是,sts test只适合分组变量的比较。如果自变量是连续变量,或者需要调整混杂因素,就应该进入Cox回归的环节。
3.4 生存表的中位生存时间
KM分析还有一个常用产出:中位生存时间。这在肿瘤研究、器械临床报告中几乎是必报的指标。
Stata里用sts list的时候,通常会顺带输出中位生存时间估计值。如果你只要中位生存时间,可以用:
sts list, med surv中位生存时间指的是生存概率下降到50%所对应的时间点。注意,如果随访结束时生存概率仍然高于50%,就报告“未达到”,而不是硬填一个随访终点的数值。我在实际审稿里见过一些报告把未达到的中位生存时间直接写成随访终点时间,这是不准确的。
4. Cox比例风险模型:极简多因素分析
4.1 stcox的基本用法和HR解读
KM分析解决的是“两组或多组的生存时间有没有差异”,但临床研究里往往需要同时考虑年龄、性别、分期、治疗方式等多个因素,这时就需要Cox比例风险模型。Stata里跑Cox模型极为简单:
stcox age treatment stage这行命令就完成了多因素Cox回归。结果里每个变量会给出风险比(Hazard Ratio, HR),标准误,z值和p值,以及95%置信区间。HR大于1表示该因素增加事件风险,小于1则表示保护因素。
解读时有个常见的坑:连续变量的HR解释。比如年龄作为连续变量纳入模型,HR如果是1.03,意思是年龄每增加1岁,事件风险增加3%。但如果你想报告“年龄每增加5岁”的风险变化,需要在建模时对变量做相应转换,或者用margin等命令计算,直接拿年龄的HR说“50岁比40岁风险高一截”是不严谨的。
Cox模型的极简流程适合用表格归纳:
| 步骤 | 命令 | 目的 |
|---|---|---|
| 声明生存数据 | stset time, failure(status==1) | 设定时间与事件变量 |
| 单因素筛选 | stcox 单个变量 | 快速看每个因素的粗略效应 |
| 多因素建模 | stcox 多个变量 | 同时调整混杂因素 |
| 检验PH假定 | estat phtest | 检查比例风险假定 |
| 输出结果 | esttab 多模型 | 整理成论文表格 |
单因素做筛选、多因素做最终模型,是临床研究里最常见的建模顺序。如果你的变量数量不多且专业上有强关联,可以直接做多因素分析,不一定非要先单因素筛选。
4.2 比例风险假定:绕不开的模型体检
Cox模型的前提是“比例风险假定(PH假定)”,意思是各组的风险比在整个随访期间保持不变。如果治疗组和对照组的KM曲线明显交叉,多半意味着比例风险假定有问题,这时Cox模型的平均HR解释起来就很别扭。
Stata的检验命令是estat phtest:
stcox age treatment stage estat phtest输出结果看全局检验的p值。p值不显著(通常>0.05)说明没有充分证据推翻PH假定,模型基本可接受。如果p值显著,说明比例风险假定被违反,需要处理。
处理PH假定问题有个相对简单的办法:做分层Cox模型,对违规的变量进行分层控制,或者把时间依存变量纳入模型。极简情况下,做完estat phtest后如果p值小于0.05,我会先把该变量作为分层变量重新建模:
stcox age treatment, strata(stage)strata()选项允许不同层的基线风险不同,但仍估计其他变量的共同效应。这比直接放弃Cox模型要优雅得多。
4.3 多个模型的表格输出
实际操作中往往要建好几个模型:模型1是单因素结果,模型2是调整了部分变量,模型3是全部调整。传统的做法是逐个跑完后手动抄写HR和置信区间到论文里,又慢又容易抄错。
Stata里的esttab命令可以帮你把多个模型并排列在同一个表格里。先用est store把模型存起来,再统一输出:
quietly stcox treatment est store Model1 quietly stcox treatment age stage est store Model2 quietly stcox treatment age stage grade est store Model3 esttab Model1 Model2 Model3, eform b(2) se(2) star(* 0.05 ** 0.01)eform选项是把输出的系数转换为HR,b(2)保留两位小数,star标记显著性。这样直接就能得到一张三个模型并列的表格,复制到Word里稍微调格式即可。
5. 亚组分析与分层分析的正确打开方式
5.1 做亚组分析前先想清楚的问题
“亚组分析”是很多临床研究里被要求做的高频操作。收到审稿意见经常能看到类似的建议:“请补充年龄亚组的分析结果”。
Stata里做亚组分析,最稳妥的思路是分别在不同的亚组内运行Cox回归,然后比较各亚组的HR大小方向和置信区间。举个例子,如果你想看治疗效果在不同年龄段里是否一致,先按年龄分组跑模型:
stcox treatment, by(agegroup)这个by()选项也可以加在stcox命令里,Stata会分别输出每个亚组的结果。还有一种做法是先stsplit再交互,但极简场景下用by()分组每条命令都清楚,结果也容易单独导出。
做亚组分析时最需要警惕的是过度解读。亚组分析本质上是探索性的,亚组样本量通常比全组小,容易出现假阳性或假阴性。我自己的习惯是:亚组结果只能作为补充证据,不能因为某个亚组显著、另一个不显著就直接断言两亚组之间存在差异。
5.2 用交互项判断组间差异:p for interaction
严谨的亚组分析,只对比各亚组P值是否小于0.05是不够的。正确做法是检验分组变量和分析变量之间是否存在交互作用,也就是所谓的p for interaction。
Stata里实现交互项很简单:
gen treatment_age = treatment * agegroup stcox treatment agegroup treatment_age或者用更简洁的因子变量语法:
stcox treatment##i.agegroup看交互项的p值。如果交互项不显著,即使不同亚组的HR看起来有差异,也没有充分证据说明治疗效应随年龄改变;如果交互项显著,才可以说存在真正的效应差异。很多论文里的“亚组分析森林图”,本质就是在展示多个亚组内的HR和置信区间,同时报告总的p for interaction。
5.3 分层分析和亚组分析不是一回事
分析时容易混淆的两个概念:分层分析(stratified analysis)和亚组分析。Cox模型中的分层分析是用strata()选项控制一个分层变量,比如不同中心、不同癌症分期,在各层内部拥有不同的基线风险函数,但核心变量的效应是全样本共同估计的。
亚组分析则是把所有样本分成几个子样本分别分析。Stata里用by()选项或循环命令都能实现。
两者的适用场景不同。当某个变量只是“nuisance variable”也就是你不想研究它,但它会影响风险基线时,用分层更合适;当你明确想知道某个变量是否改变核心效应的方向或大小时,用亚组分析更合理。
实际报告里,分层分析的结果喝亚组分析可以一起出现:先用多层模型调整个重要分层变量,再做亚组分析展示效应修饰作用。
6. 高频实战命令:外部命令、最大值最小值与meta分析前缀
6.1 外部命令从哪里找
Stata自带命令覆盖了基础统计功能,但一些特殊分析需要安装外部命令才能用。比如某些高级图表、特定类型的meta分析工具包,都不是默认自带的。
遇到不认识的命令,先不要慌。Stata提供了一个官方社区的外部命令仓库,常见的命令都能在这里找到。安装方法是在Stata命令窗口输入:
ssc install 命令名比如你可能在搜“网状meta分析Stata操作”时看到过ftool之类的命令,这类第三方命令能不能用、版本是否支持,最好先用findit确认一下官方说明再安装。我一般不推荐从不明来源下载.ado文件直接丢进Stata目录,一旦文件冲突或版本不匹配,整个软件都会出问题。
判断一个外部命令是否可用的方法很简单:输入ssc hot,查看最近流行和更新的命令列表;输入findit 关键词,查看相关命令的介绍与安装说明。极简思路是,能跑内置命令就用内置命令,不要为了炫技盲目安装一堆外部包。核心的生存分析流程用默认命令就完全可以覆盖。
6.2 最大值最小值与描述统计命令
有时候做数据清洗时需要快速看变量的取值范围,比如检查生存时间是否有异常值、年龄是否有录入错误。这里高频使用的就是summarize命令:
summarize time age, detail加了detail选项后会输出最小值、最大值、中位数、百分位数等详细信息,通过min和max基本能一眼发现异常。
如果你需要生成一个新的变量,比如把年龄大于80岁的患者标记出来,就会用到egen命令:
egen agemax = max(age) generate old = (age==agemax)更常用的做法是直接用max()函数在generate里完成。但egen的优势是可以按组生成统计量,比如按治疗组生成各组最大年龄。极简场景下,记住这两条命令就够用了:
summarize 变量, detail egen 新变量 = min(变量) / max(变量)6.3 meta分析前缀是另一个方向
很多人在搜索Stata生存分析时,会同时搜到meta分析相关的命令。确实,Stata的meta分析模块也很成熟,尤其是network meta分析,在Stata里有一整套配套命令和绘图功能。和生存分析的问题类型不同,meta分析是用来合并多项独立研究的效应量,回答的是“多个研究综合起来总体效应怎样”的问题,和单个队列的生存分析不是一回事。
如果你的研究路线更偏向meta分析,需要系统学的是meta set、meta summarize、network graph等一套流程;但如果你当前的任务就是分析一个随访队列的生存数据,那专注于sts和stcox这一套就够了。不要被搜索时同时出现的热词带偏,不同分析目的用不同的命令体系。
7. 常见报错与排查技巧实录
7.1 “no observations”或“variable not found”防不胜防
新手跑生存分析时最容易出现的报错之一是variable not found。核对命令,逐一检查变量名是否输错、大小写是否一致。Stata是区分大小写的,例如Time和time是两个不同的变量。
另一种常见情况是“no observations”,多见于stset后做亚组分析。比如你写了by(treatment),但treatment变量里存在缺失值,缺失值所在的样本会被自动剔除,样本量突然变小。排查思路是先用tab treatment查看分组变量的缺失情况,再用count确认实际有效样本数。
7.2 stset后遗症的困扰
stset声明数据结构后,Stata会自动生成_st、_d、_t等变量。如果后续你手动删除或覆盖了这些变量,Stata很可能弹出奇怪报错,或者生存分析命令无法运行。
解决办法很简单:重新stset一次,把数据结构恢复回来。如果你在stset之后又执行了sort、merge或其他数据处理操作,也建议重新stset确认数据状态。就算数据没有变化,重复执行stset的成本也极低,比因底层变量被改动导致报错盲查半天更省时间。
7.3 日期变量直接放进模型的问题
部分初学者会把日期变量本身当作生存时间纳入模型,比如直接把“手术日期”和“终点日期”放进模型,而不是先算出时间差。这样跑出来的结果几乎肯定是错的,因为日期变量在Stata里是一串连续的日期代码,不是以天、月为单位的随访时间。
正确的做法是先用时间差生成生存时间变量。Stata里最简单的日期差计算:
gen double time_days = end_date - start_date gen time_months = time_days / 30.44建议把时间单位统一好再进入分析,比如统一用月或统一用天。不要在同一个模型里混用不同时间单位,模型结果的解释会非常别扭。
7.4 surv曲线画出来了但图形“不交叉”
KM曲线画出来如果几条线紧紧贴在一起,甚至完全重合,不一定是没有效应,可能是分组变量没有正确载入或图形选项写错了。先检查by()的变量类型是否数值型,如果分组变量是字符串,Stata会提示类型不符,需要先encode转换。
还有一种可能性是真没有差异。这时候不要强行通过选项调整y轴范围来制造“视觉差异”。p值就在那里,不显著就是不显著。改变刻度范围放大差异感是不诚实的数据展示方式,审稿人很容易识破。
常见报错与解决思路整理成一个速查表:
| 现象 | 可能原因 | 解决方式 |
|---|---|---|
| stset后无样本 | 时间变量存在缺失 | 用missing()检查并填补或剔除 |
| 变量找不到 | 变量名拼写或大小写不当 | tab命令确认准确变量名 |
| 交互项不显著 | 亚组差异是随机波动 | 用p for interaction判断再下结论 |
| 图形导出不清晰 | 分辨率设置不够 | export时设置width=2400以上 |
| 中位生存时间不显示 | 生存概率未降至50% | 报告为“未达到”并说明随访时间 |
| HR结果很异常 | 时间单位不一致或编码反转 | 检查事件变量0/1与stset的failure定义 |
7.5 极简流程的三个省时习惯
最后分享三个我实际操作中形成的习惯。第一个,脚本化管理。所有分析都写成.do文件,不要直接在命令窗口手输。生存分析经常要反复调整变量和模型,脚本化之后改个变量名重跑一遍就行,还方便留档应对审稿人的原始数据核查要求。
第二个,分析过程记录。跑完每个模型以后顺手执行describe和count,确认样本量没有莫名减少,然后把结果用log或outfile保存起来。等写论文时你会发现,这些零散的输出要比“当时我记得结果是多少”可靠得多。
第三个,报错先读英文。Stata的报错信息大部分是英文的,很多问题直接读一遍报错就能知道是什么原因。不要遇到报错就盲目搜索引擎,你先理解这句报错在说什么,往往省下的时间按小时计。
我自己的习惯是,从头到尾的极简流程固定写成几行do文件:stset、sts list、sts graph、stcox、estat phtest、esttab。整个生存分析的主线就这几步,数据量再大,跑的也很快。熟练之后你会发现,真正花时间的不是Stata命令怎么敲,而是数据清洗阶段怎么确保时间变量和事件变量录对。基础数据没问题,后面的分析全是水到渠成的事。