什么时候用它
生存分析处理的是**"事件发生了没有"和"什么时候发生"合在一起**的数据,而且必然带删失——有人到研究结束还没发生事件,你只知道"至少活到了这一天"。
普通的均值比较处理不了这种数据:把删失的人当成"没事件"会低估风险,直接扔掉会引入偏倚。Kaplan-Meier 曲线是专门为此设计的:每次有事件发生才下降一级,删失只是把这个人从风险集里移走,不让曲线掉。
怎么读
- 起点是 1。 时间 0 时所有人都还没发生事件。
- 每级台阶是一次事件。 曲线是阶梯状的,不是平滑的——它只在观察到事件的那一刻下降。
- 看两条线分不分得开。 分得越开,两组的生存经验差别越大。
- 看置信区间。 后期变宽是正常的:风险集里人少了。
- 看风险表。 这是最容易被跳过、但最要紧的一行——某个时间点之后只剩几个人,那之后的曲线就别当真了。
- 最后看 p 值。 log-rank 概括的是整个随访期的差异,不说明差多少。
常见陷阱
- 事件编码是最容易错的一步。
survival::lung 里 status 是 1 = 删失、2 = 死亡,但大多数函数期待的是 0 = 删失、1 = 事件。Tessera 的 lung_survival 已经转换为明确的 event = 0/1,recipe 仍会在作图前验证它——换自己的数据时也要保留这一步。
- 删失不让曲线下降。 曲线平着走一段不代表"没人退出",可能只是那段时间里退出的人都是删失。
- 尾部不可信。 风险集只剩个位数时,一次事件就能让曲线掉一大截。
- 中位生存可能估不出来。 曲线没降到 0.5 就没有中位数,报告时要写"未达到"而不是留空。
- log-rank 假设风险比恒定。 两条曲线交叉时它的检出力会很差——那种情况该看的是别的方法。
- p 值不是效应量。 要说"差多少"得给风险比和置信区间。
- 别在图上量某个时间点的生存率。 曲线是阶梯状的、还压着置信带,肉眼取值误差很大。要 1 年、2 年生存率就直接算:
summary(fit_km, times = c(365.25, 730.5)),标准误和置信区间一并给出。
ggsurvplot() 的返回值不能直接 ggsave()。 它是个复合对象而不是 ggplot,得存 p$plot;而一旦开了 risk.table = TRUE,风险表在 $table 里、不在 $plot 里,直接存主图会把风险表丢掉——要合起来得走 arrange_ggsurvplots()。方法 A 没这个问题,ggsurvfit() 返回的就是普通 ggplot 对象。
配色
两组使用 Tessera 的 gene_red:黑 #000000 和绯红 #B11522。这是很克制的二元对照,红色只负责把第二组从中性基线中拎出来;置信区间沿用曲线本色并降低透明度,不再额外引入一套浅色。
曲线图的配色要点是曲线比填充重:置信区间是背景信息,深浅必须拉开,否则两条带子叠在一起时谁也看不清。
配方
数据准备
示例直接读取 Tessera Toy lung_survival。原始 lung 的编码转换已经固化在数据生成脚本中;recipe 再验证一次事件列,避免换数据后静默画错。
library(survival)
library(ggsurvfit)
library(dplyr)
data_km <- read.csv("content/tessera/data/csv/lung_survival.csv") |>
mutate(
event = as.integer(event),
sex = factor(sex, levels = c("male", "female"))
) |>
filter(!is.na(time_days), !is.na(event), !is.na(sex))
stopifnot(all(data_km$event %in% c(0L, 1L)))
# survfit2() 来自 ggsurvfit,比基础的 survfit() 更好接后面的绘图层
fit_km <- survfit2(Surv(time_days, event) ~ sex, data = data_km)
# log-rank:rho = 0 是标准 log-rank;rho = 1 是 Peto 权重,更看重早期事件
sd_logrank <- survdiff(Surv(time_days, event) ~ sex, data = data_km, rho = 0)
p_logrank <- pchisq(sd_logrank$chisq, df = length(sd_logrank$n) - 1,
lower.tail = FALSE)
方法 A · ggsurvfit()
一张完整的 KM 图有五件东西:曲线、置信区间、风险表、中位线、检验 p 值。
#| fig: km
#| fig-width: 8
#| fig-height: 7
pal_sex <- c(male = "#000000", female = "#B11522")
ggsurvfit(fit_km) +
labs(
title = "Overall Survival by Sex",
subtitle = "Tessera Toy: lung_survival · Kaplan-Meier with 95% CI",
x = "Time (days)", y = "Survival probability"
) +
scale_color_manual(values = pal_sex, breaks = names(pal_sex),
labels = c("Male", "Female")) +
# 置信区间沿用曲线颜色,只降低存在感;不再引入 recipe 外的颜色
scale_fill_manual(values = pal_sex, guide = "none") +
add_confidence_interval(type = "ribbon", alpha = 0.14) +
add_risktable() + # 风险表 —— 判断尾部可不可信全靠它
add_quantile(y_value = 0.5, color = "gray50", linewidth = 0.75) + # 中位生存参考线
scale_ggsurvfit() + # 把 y 轴调成 0–1 的百分比刻度
ggplot2::theme_classic() +
theme(
legend.title = element_blank(),
legend.position = "bottom",
plot.title = element_text(size = 16, face = "bold")
) +
geom_text(
aes(x = 500, y = 0.8, label = sprintf("log-rank p = %.3f", p_logrank)),
inherit.aes = FALSE, # 不继承曲线的分组映射,否则每组各画一次
hjust = 0, size = 4.2
)
survival_curve-km
方法 B · 累计风险
同一批数据的另一个视角:不看"还没发生事件的概率",看"风险累积了多少"。
#| fig: cumhaz
#| fig-width: 8
#| fig-height: 6
library(survminer)
fit_km0 <- survfit(Surv(time_days, event) ~ sex, data = data_km)
ggsurvplot(
fit_km0,
data = data_km,
fun = "cumhaz", # 关键参数:换成 "event" 是累计发生率,默认是生存概率
conf.int = TRUE,
conf.int.style = "ribbon",
conf.int.alpha = 0.20,
palette = unname(pal_sex),
legend = "bottom",
legend.title = "",
legend.labs = c("Male", "Female"),
xlab = "Time (days)", ylab = "Cumulative hazard",
title = "Cumulative Hazard by Sex",
ggtheme = ggplot2::theme_classic(base_size = 12),
break.time.by = 200 # x 轴每 200 天一个刻度
)
survival_curve-cumhaz