Stata生存分析极简指南：从KM曲线到Cox回归

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(status1)意思是设定生存时间变量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(status1) 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(status1)设定时间与事件变量单因素筛选stcox 单个变量快速看每个因素的粗略效应多因素建模stcox 多个变量同时调整混杂因素检验PH假定estat phtest检查比例风险假定输出结果esttab 多模型整理成论文表格单因素做筛选、多因素做最终模型是临床研究里最常见的建模顺序。如果你的变量数量不多且专业上有强关联可以直接做多因素分析不一定非要先单因素筛选。4.2 比例风险假定绕不开的模型体检Cox模型的前提是“比例风险假定PH假定”意思是各组的风险比在整个随访期间保持不变。如果治疗组和对照组的KM曲线明显交叉多半意味着比例风险假定有问题这时Cox模型的平均HR解释起来就很别扭。Stata的检验命令是estat phteststcox 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选项是把输出的系数转换为HRb(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 (ageagemax)更常用的做法是直接用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时设置width2400以上中位生存时间不显示生存概率未降至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命令怎么敲而是数据清洗阶段怎么确保时间变量和事件变量录对。基础数据没问题后面的分析全是水到渠成的事。