ARTICLE DETAIL

资讯详情

深耕编程入门与网站建设的一线实战洞察。

Stata生存分析极简教程:从KM曲线到Cox回归实战

Stata生存分析极简教程:从KM曲线到Cox回归实战 做生存分析这件事很多人的第一反应是“这是临床研究才用的东西”或者“得先学懂一堆数学公式再说”。实际上只要手里有Stata哪怕你连风险函数、累计风险曲线这些名词都没听过也能在十分钟内跑出像样的生存分析结果。我这些年用Stata做过不少生存数据建模从最简单的单组生存曲线到多因素Cox回归再到亚组分析核心步骤其实就那么几条命令真正的难点反倒全在数据准备和结果解读上。这篇文章的目标很明确用极简的思路带你完整走一遍Stata生存分析流程。你不需要提前啃统计教材也不需要装一堆乱七八糟的外部命令包我用的全部是Stata内置命令版本从Stata 14到Stata 18都能跑。文章里我会先讲清楚生存分析到底在算什么然后把数据格式、stset声明、KM曲线、Log-rank检验、Cox回归这些关键环节一个个拆开中间穿插我实际踩过的坑和报错排查记录最后再补几个高频实用技巧比如亚组分析、宽表转长表以及外部命令的安装逻辑。无论你是医学生、公卫从业者还是社科、经管方向的研究生这套流程拿过去就能直接用。1. 捅破窗户纸生存分析到底在干什么1.1 生存分析的本质其实就是三个问题别被“生存分析”这四个字吓住它的底层逻辑特别朴素。所谓生存数据本质上包含三个信息一个人或一个对象从某个起点开始到某个事件发生为止经历了多长时间到观察结束时这个事件到底发生了没有以及这个人或对象本身有哪些特征。回想一下你手上的数据不管来自临床随访、设备寿命测试还是用户流失分析只要你能回答三件事——观察了多久、结局是什么、有哪些影响因素——你就能做生存分析。我习惯用一个例子给学生讲这个问题假设你想研究某种新药对术后复发的影响每个病人从手术那天开始观察记录他复发的时间有些人观察期结束后还没复发有些人中途失访了还有些人因为其他原因去世了。对于没复发的、失访的、死于其他原因的我们都只知道“至少到某个时间点还没复发”但并不确定他未来会不会复发。这种数据普通回归处理不了因为你连“因变量”都不完整但生存分析天生就是干这个的它能把这种“不完整的信息”全部利用起来。这也是生存分析区别于普通统计方法的根本原因它不是只看事件有没有发生而是把“事件发生的时间”也纳入建模。时间越长还没出事本身就是一个信号。看设备故障也一样一台机器运行了800天才坏和一台运行了20天就坏了信息量完全不同如果只用一个0/1变量“坏没坏”来建模等于把最有价值的时间信息丢掉了。1.2 Stata做生存分析真的不用额外装什么不少人在网上搜“stata下载”“stata安装包”第一反应是到处找资源其实Stata的生存分析模块是内置的根本不需要额外安装。核心命令就五个stset、sts list、sts graph、sts test、stcox。前两个是数据声明和查看中间两个做描述性分析和组间比较最后一个做多因素回归。整个流程下来你甚至不需要会写任何复杂的程序只要会点鼠标或者敲几行命令就够了。当然有些场景确实需要外部命令比如以后你想画更漂亮的森林图或者做PSM倾向性评分匹配或者像热词里提到的ftool命令这些才需要通过ssc install或者findit去安装。但那是后话对一篇极简入门来说先把内置命令用熟比什么都重要。我用过太多版本了从Stata 14到Stata 18这几个生存分析命令的语法几乎没有变化这意味着你学会一套以后换电脑换版本都不用重新学。2. 数据是爹生存分析的数据准备与stset声明2.1 三个关键变量缺一个都不行做生存分析之前先检查你的数据里有没有这三个要素时间变量、结局变量、分组或解释变量。时间变量是一个连续变量表示从起点到事件发生或观察结束所经历的时间单位可以是天、月、年但整个数据集必须统一。结局变量通常用0和1编码1代表事件发生了比如复发、死亡、故障、流失0代表删失也就是到观察结束时事件还没发生。第三个是你要比较或研究的因素可以是分组变量也可以是年龄、性别、血压这类连续或分类变量。这里我特别想强调一个我在实际中反复遇到的坑很多人把结局变量的编码搞反或者把删失和事件混在一起。比如有人图省事把“未复发”编码成0把“复发”编码成1这没错但有人把“失访”和“未复发”都编码成0结果模型把失访的人当成“观察全程都没复发”来处理时间信息会失真。正确的做法是只要事件没发生不管是因为观察期结束、失访还是其他原因结局都记为0同时保留它的观察时间。Stata并不知道你为什么删失它只需要知道“这个人在这个时间点之前没出事”剩下的交给统计模型处理。2.2 用最大值最小值命令摸清数据底细在声明数据之前我强烈建议先花一分钟摸一摸数据的底细。最常用的就是summarize命令配合detail选项看最大值和最小值。比如你的时间变量叫time直接跑summarize time, detail输出里会给出时间变量的最小值、最大值、分位数、均值等一堆信息。这一步至少能帮你发现三类问题时间变量有没有负数或0时间变量的单位是不是混着天和月有没有极端大值可能是记录错误。我曾经处理过一份数据time变量最大值显示是3650结果一查发现有人的随访时间被录成了3650天实际应该是365天多按了一个0。如果没先看最大值就直接建模整条曲线都会被这个异常值带偏。另外summarize配合return list还能把最大值、最小值提取出来备用。比如你想在数据清洗时把时间大于某个阈值的样本筛掉summarize time, detail local max_time r(max) drop if time 365*5这就用到了极简数据处理里的一个小技巧——先看范围再决定清洗规则。生存数据里时间变量的质量直接决定分析质量这一步千万别跳。2.3 stset声明让Stata听懂你的数据Stata里做生存分析的第一步永远是stset它的作用就是告诉Stata哪个变量是时间哪个变量是结局结局的编码是什么。基本语法长这样stset time, failure(status1)这条命令的意思是时间变量是time结局变量是status其中status1代表事件发生。Stata收到声明之后会自动帮你在后台生成一组生存分析专用变量并用_st开头存起来。你可以通过stdescribe看看数据的基本情况或者通过stsum看看不同组的观察人时、事件数、发病率。stset这一步最核心的就是failure()这个选项它决定了“什么才算事件”。有几种常见写法* 终点事件编码为1 stset time, failure(status1) * 终点事件编码为某个具体数值比如2 stset time, failure(status2) * 不写failure选项默认status1为事件 stset time我个人的习惯是永远显式地写出failure()条件哪怕默认值正好是1也写出来因为代码可读性更好过三个月你再回来看脚本一眼就知道当时的事件定义是什么。还有一个必填项要留意如果数据里包含id变量比如同一个病人有多次随访记录建议加上id(idvar)选项这能帮助Stata正确处理重复记录。多行数据一定要先搞清楚每条记录是一个样本还是一次随访这决定了要不要指定id()很多新手在这里翻车。3. 极简三连KM曲线、Log-rank检验与Cox回归3.1 KM生存曲线五步跑通三步看懂数据声明好之后最直观的第一步就是画KM生存曲线。KM曲线全称Kaplan-Meier生存曲线本质上是把每个事件发生时间点上的存活概率连成一条阶梯线每一步下降都对应一次事件发生。画图命令极其简单sts graph如果只有一个整体这条命令会输出一条阶梯下降的曲线。想看分组对比加上by()选项sts graph, by(group)出来的图就是两条或多条阶梯线一眼能看出哪组“活得好”。但我一般会再补两个选项让它更好看、更实用。一个是failure把纵轴从生存概率换成累积失败概率也就是把曲线反过来看另一个是ci给曲线加上置信区间带。sts graph, by(group) failure ci曲线输出之后怎么读我总结成三步第一步看整体趋势曲线是快速下降还是平缓下降下降越陡说明事件发生得越早越集中第二步看组间差距两条曲线是早早分开、到后期又合拢还是一直保持距离前者提示早期效应后者提示持续效应第三步看阶梯的跳跃位置在某个时间点出现一个大台阶说明那个时间点附近有大量事件集中发生比如临床上常见的术后30天、90天、1年复查节点。这里有个细节值得单独说sts graph默认会在删失点上画小短竖线这对判断数据质量很有帮助。如果竖线特别多说明删失比例高你要有心理准备后续统计检验的把握度可能不足。如果你不想显示删失标记可以加noshowno之类的选项但我不建议这么做信息多一点没坏处。3.2 Log-rank检验两组差异到底显不显著曲线分开看了还需要一个统计量来回答“这个差异是不是真的还是随机波动”。最常用的就是Log-rank检验一句话命令sts test group输出会给出卡方值和P值。当P小于0.05我们通常认为组间生存曲线差异有统计学意义。Log-rank检验的核心思想是在每个事件发生时间点计算各组的期望事件数再和实际事件数对比最后汇总成一个卡方统计量。它给所有时间点相同的权重所以对后期差异比较敏感。那有没有“不极简”但更灵活的替代如果你发现两条曲线在早期就分开后期又缠在一起这时候Log-rank检验可能测不出来因为它对早期差异不敏感。Stata提供了wilcoxon选项用的检验对早期事件加更大的权重sts test group, wilcoxon我的经验是常规场景优先看Log-rank结果如果曲线明显早期分化而Log-rank不显著再用Wilcoxon做敏感性分析然后在一句话里报告两种结果。这样审稿人和导师都会觉得你想得比较周全。3.3 Cox回归多因素下看风险比KM曲线和Log-rank检验说到底都是单因素分析只能比较一个分组变量。要想同时校正年龄、性别、合并症等多个因素必须上Cox比例风险回归。极简版命令stcox age gender group输出里最关键的是每行变量对应的HR值风险比、95%置信区间、P值。HR大于1表示该变量增加事件风险小于1表示降低风险。比如group的HR是0.595%CI为0.3~0.8P0.004意思是在校正了年龄和性别后实验组的事件风险是对照组的一半。这里我要多写几句关于Cox模型的理解。很多新手拿到结果只会念“HR0.5P0.05”但问一句“这个0.5是怎么来的”就答不上来。Cox模型不直接估计生存时间而是估计风险函数它默认不同个体在任何时间点的风险成固定比例这就是“比例风险假定”。Stata提供estat phtest命令来检验这个假定是否成立stcox age gender group estat phtest, detail如果检验P值小于0.05说明比例风险假定不成立这时候结果需要谨慎解读可以考虑加时变协变量或者分层Cox模型。不过对极简入门来说先把基本结果跑对再去考虑这些进阶问题顺序不能乱。4. 我踩过的坑报错排查与曲线救国的现场实录4.1 高频报错与排查速查表用Stata做生存分析最让人头疼的不是统计方法本身而是莫名其妙的报错。我把自己踩过和帮别人解决的报错整理成了下面这张表给正在被报错折磨的朋友一个快速定位入口。报错信息出现场景排查思路invalid syntaxstset或sts命令检查逗号位置、变量名是否含特殊字符常见的是命令里漏了逗号no observationssts list/graph数据是否为空stset是否已执行是否用了if条件但条件筛选后样本为0variable time not foundstset时间变量的名字写错了注意大小写Stata变量名区分大小写failure variable must be coded 0/1stset结局变量不是0/1或者存在缺失值用tab status, missing检查r(198) invalid syntax任意命令通常是命令拼写错误或安装了不兼容的外部命令先查语法time variable has nonpositive valuesstset时间变量里有小于等于0的值生存时间必须为正数第4条真的特别常见。很多数据里的status变量用的是1和2比如1存活、2死亡而不是0和1这时候Stata会报错。解决办法是重新编码* 把1和2编码成0和1 recode status (10) (21), gen(event)或者更简单的在stset里直接用failure(status2)根本不用改数据。这个小技巧帮我省了无数次重新编码的时间。4.2 三个很容易被忽略的实操细节第一个细节是时间单位的统一。很多数据集里的时间是东拼西凑出来的有的人录天有的人录月还有人录年。如果混着用KM曲线的时间轴会乱得没法看。我一般拿到数据先做一次单位校准比如全部转换成天数月数乘以30.44年数乘以365.25。虽然不精确但只要全数据统一转换结果就一致可比。第二个细节是删失比例过高的问题。如果一份数据里删失比例超过80%KM曲线往往长期徘徊在高位组间差异很难测出来。这不是代码问题是数据本身的信息量不够。Stata的sts list可以输出各时间点的风险集人数我建议在报告KM曲线时附上关键时间点的风险集人数number at risk这能直接反映后续曲线的可信度。第三个细节是结果保存。很多人跑完stcox就直接把结果截图存进论文草稿后面换数据重跑图表全部要手动重做。强烈建议用est store和esttab把多个模型整理成一张表格est store m1 stcox age gender group est store m2 esttab m1 m2, b(%9.2f) ci(%9.2f) star(* 0.05 ** 0.01)这样模型结果就能一键导出成规范的三线表样式避免手动誊抄出错。5. 进阶小技巧亚组分析、数据转换与命令扩展5.1 亚组分析按人群拆开看更清楚“亚组分析”是近期热词公共卫生和临床研究里尤其常用核心思路是想知道某个效应在特定人群里是否更明显或更弱。在Stata里面极简做法是用if条件限定样本范围。比如只分析男性stcox age group if sex1或者只分析60岁以上人群stcox age group if age60这种写法的好处是极简、直观、不容易出错坏处是如果亚组很多比如按性别、年龄段、疾病分期拆出六七个组就得写六七遍代码效率低还容易复制粘贴出错。这时候我习惯用循环forvalues i 1/3 { stcox age group if stagei est store stage_i } esttab stage_1 stage_2 stage_3, b(%9.2f) ci(%9.2f)三行代码把三个分期亚组的模型全部跑完并汇总成一张表。这个思路也适用于KM曲线和Log-rank检验先sts graph, by(group)看整体再用if条件分组画亚组曲线。不过要提醒一句亚组分析的样本量通常更小跑出来的置信区间会很宽解释时要特别克制别把一个不显著的亚组结果讲成“有效趋势”。5.2 宽表转长表准备数据时长用的reshape生存分析的数据结构经常是“一人一行”的长表每个样本一行包含时间、结局、协变量。但有时候你拿到的原始数据是宽表比如每次随访记录成一列这种情况下需要先reshape long转成一行一个观测才能做生存分析。举个例子宽表里每个病人有3次随访记录* id 是病人编号f1 f2 f3 是三次随访的复发状态 reshape long f, i(id) j(followup)转成长表后每个病人有3行数据。这时候就能用stset声明并利用时间变量followup构造生存数据。宽表转长表这块儿很多人卡在不知道j()该怎么写。记住一个口诀j()后面跟的是“把多列变成长表后那一列新变量叫什么名字”比如j(visit)就是把f1、f2、f3的序号1、2、3存成变量visit。如果你在做数据管理时还需要判断每个样本的最长随访时间最大值最小值命令还能派上用场bysort id: egen max_followup max(followup)这一行就给每个样本标记出他到底被随访到第几期后面做时间相关变量构造非常方便。5.3 外部命令与ftool少走弯路地安装扩展命令Stata内置命令覆盖90%以上需求但总有偶尔要用的扩展命令。比如热词里的ftool以及psmatch2、meta相关的网络meta分析命令这些都需要自己安装。极简安装逻辑只有两条路ssc install或findit。知道命令名字的时候直接用ssc install ftool不确定命令名字只想搜索相关功能时用findit survival plotStata会打开一个搜索结果窗口点蓝色的“click here to install”就能装上。这里有个经验装外部命令前先看一眼它的依赖包很多命令安装时会提示还需要装别的包比如ftool可能依赖estout或moremata。遇到这种情况别慌Stata的ssc install有时候会自动装依赖如果没自动装就手动把提示里出现的其他包也ssc install一遍基本都能解决。不过我要多说一句极简生存分析的场景下外部命令不是必需品。我见过不少人还没搞懂stcox就到处找美化KM曲线的外部命令结果图是画好看了底层的统计判断反倒说不清楚。先把内置命令吃透需要进阶时再装扩展包这个顺序才不会跑偏。做生存分析这些年我最大的一个体会是统计软件只是工具真正决定分析质量的永远是你对自己数据的理解。Stata的极简之处在于它把复杂的生存分析压缩成了几个语义清晰的命令但每个命令背后都有它的适用边界和前提假设。你不需要立刻弄懂所有数学推导但至少要知道每一步在算什么、结果怎么解读、什么情况下结果可能不靠谱。把这几点做到位生存分析这门手艺就算入门了。最后再分享一个小习惯每次跑完分析把stset声明、模型命令、结果表格完整地存档成一个.do文件加好注释。三个月后回来看你会感谢当时的自己。
返回列表