R 与 RStudio 生存分析实战:Kaplan–Meier 曲线和 Cox 回归
安装 R 与 RStudio,使用 survival 包示例数据绘制生存曲线、拟合 Cox 模型并检查比例风险假设。
生存分析处理的是“从起点到某事件”的时间资料,还需要记录观察结束时是否已经发生事件。Kaplan–Meier 曲线适合展示不同组的时间分布,Cox 模型可在指定协变量下估计相对风险。下面的代码使用 survival 包公开的教学数据,只练习计算,不给出治疗建议。
先安装 R,再安装 RStudio
从 R 项目官网安装 R,再从 Posit 官网安装 RStudio Desktop。R 是计算环境,RStudio 是编辑和运行代码的界面;如果只装了 RStudio 却找不到 R,先核对 R 是否安装及其路径。RStudio 的入门文档介绍了脚本、控制台与输出窗格。
打开 RStudio,新建一个 .R 脚本,在控制台安装并加载包:
install.packages("survival", repos = "https://cloud.r-project.org")
library(survival)
packageVersion("survival")
install.packages() 只需在初次安装或更新时运行;每次新会话使用 library() 加载即可。正式分析应记录 R、RStudio 和软件包版本,以便以后复算。
认识事件与删失编码
survival 包自带的 ovarian 数据含有随访时间 futime、事件状态 fustat、组别 rx 等字段。执行 ?ovarian 和 table(ovarian$fustat),先核对字段含义与编码,再写模型。数据集文档说明了变量来源;这些公开示例资料也不能替代一个新研究自己的数据字典。
图 1:虚构对象的随访时间线;实心点是观察到事件,空心点表示随访结束时未观察到事件。
一个真实数据集至少应明确三列:起点到末次观察的时间、是否发生目标事件、分组或协变量。时间起点若有人从入组算、有人从治疗开始算,曲线可能失去可比性。状态字段也常见 0/1 与 1/2 两种编码;Surv() 对常见编码有处理规则,但研究者必须自己检查原始数据字典,不能靠软件猜测。Surv 文档列出了右删失资料的状态编码。
library(survival)
data(cancer, package = "survival") # 加载包含 ovarian 的示例数据
str(ovarian)
table(ovarian$rx, ovarian$fustat)
km <- survfit(Surv(futime, fustat) ~ rx, data = ovarian)
plot(km, col = c("#4f8070", "#c38d69"), lwd = 2,
xlab = "随访时间", ylab = "估计生存概率")
legend("topright", legend = levels(factor(ovarian$rx)),
col = c("#4f8070", "#c38d69"), lwd = 2, title = "rx")
Surv() 组合随访时间与事件状态,survfit() 估计各组曲线。曲线在后期样本很少时会变得不稳定;报告图形时应同时给出样本数、事件数和合理的风险人数表。此处只展示最基本的绘图步骤。survfit 文档说明了公式和输出。
要检查曲线在各时间点有多少人仍处于风险集中,可以运行:
km_table <- summary(km, data.frame = TRUE)
head(km_table[, c("strata", "time", "n.risk", "n.event", "surv")])
n.risk 是时间点前仍处于观察且尚未发生目标事件的人数;n.event 是该点发生事件的人数。删除一条记录、改变缺失值处理或改动随访起点,都会改变这些数。曲线尾部若风险人数很少,纵轴看似明显分开的两条线也可能具有很大不确定性。summary.survfit 文档说明了这些字段。
拟合 Cox 模型并检查前提
fit <- coxph(Surv(futime, fustat) ~ rx + age, data = ovarian, x = TRUE)
summary(fit)
ph_check <- cox.zph(fit)
print(ph_check)
plot(ph_check)
summary(fit) 会显示进入模型的例数、事件数、回归系数、风险比及区间。风险比的含义取决于变量编码和参考组,不能仅看数字大于或小于 1。cox.zph() 与图形用于检查比例风险假设;若假设不合适,应结合研究问题考虑不同模型,而不是忽略诊断。coxph 文档和比例风险检验文档给出函数定义。
具体解释前还要检查 summary(fit) 中的有效例数与缺失情况。若 age 有缺失,Cox 模型可能使用的对象少于 Kaplan–Meier 曲线,此时两个输出并非基于完全相同的人群。风险比也不是某个时间点的绝对风险,更不等于生存时间“增加了多少百分比”。要报告绝对时间结局,应明确选择了哪一个时间点或其他合适的指标。
一个可复查的结果段可按顺序写:研究对象与起点、事件定义、随访长度、各组 N 与事件数、曲线及风险人数、Cox 模型中的协变量与参考组、风险比和 95% 区间、比例风险检查与限制。若研究中存在竞争事件、延迟入组或重复事件,上面的单事件右删失示例就不够,需要使用更适合的数据结构和方法。
真实研究还要处理时间起点、失访、缺失、共变量预设和竞争事件。尤其是小样本教学数据,估计区间可能很宽;不要把练习输出写成临床效果结论。
最后保存脚本与软件环境信息,而不是只留一个 PNG:
sessionInfo()
saveRDS(list(km = km, cox = fit), "survival-practice.rds")
sessionInfo() 记录 R 与已加载包的版本;.rds 可以保存当前拟合对象,脚本则说明对象如何产生。公开发布代码前,确认输入数据是否允许共享;真实患者数据应按机构要求处理,不要把其原始表格连同脚本上传公开仓库。
运行模型前加一段数据核查
下面的检查不会代替研究设计,但能及早发现明显的录入问题:
stopifnot(all(ovarian$futime >= 0, na.rm = TRUE))
table(ovarian$fustat, useNA = "ifany")
table(ovarian$rx, useNA = "ifany")
colSums(is.na(ovarian[c("futime", "fustat", "rx", "age")]))
正式数据还应检查日期先后、同一对象是否有多行、事件日期与末次随访日期是否矛盾,以及是否有人在“研究起点”之前已发生目标事件。若观察对象的时间单位有的用天、有的用月,模型仍可能运行,但所得曲线和风险比没有可解释的统一时间尺度。
曲线与 Cox 模型分别回答什么
Kaplan–Meier 曲线以组别描述观察到的事件时间分布;Cox 模型在指定协变量条件下估计相对瞬时风险,两者并非同一指标。曲线分开也不等于因果作用;模型加入年龄等变量后仍可能有未测混杂。若研究是随机试验,还应依照试验方案报告分析人群;若是观察性研究,更要解释组间基线差异与混杂控制。
删失意味着在其观察结束前没有看到目标事件,可能是研究随访结束,也可能是失访。对非信息性删失等前提不能只靠软件默认值来保证。若失访与病情严重程度有关,结果可能偏倚。小型教学数据可以帮助掌握命令,但正式论文需要事先计划结局定义、随访窗口、缺失处理及敏感性分析。保存每次分析对应的代码提交或文件版本,才能说明不同稿件图表为何不一致。
作为练习的最后一步,关闭 RStudio 后重新打开脚本,在一个新会话中从第一行运行到最后一行。若代码依赖了上次会话中手工创建、却没有写入脚本的对象,会在这一步暴露。能从干净会话再现图和模型,是“可复现分析”的基本检查;如果包版本变化造成结果差异,可考虑使用项目级环境记录工具,并在论文方法中写清环境版本。不要把 .RData 中偶然保存的工作空间当成分析流程本身。
保存图时也应给纵轴、横轴写出定义和单位,并区分“随访到该时间仍未发生事件的估计概率”与“事件发生概率”。若图中未展示风险人数,读者可能误读尾部少数对象造成的陡降;正式报告宜配风险人数表或在图注解释尾部样本量限制。