如何在R中将Kaplan-Meier生存曲线转换为以连续预测变量为X轴、生存概率为Y轴并按指定时间点着色的可视化图

我完全懂你想要的这种可视化——把传统KM图的横轴换成你的连续预测变量，纵轴保留生存概率，用颜色区分不同时间点的生存趋势，这样能直观看到整个连续变量范围内的生存概率变化，代价是不再展示时间维度的整体分布。确实survfit()和ggsurvplot()默认做不了这种转换，但我们可以手动提取KM模型的数据，用ggplot来实现，完全基于KM的结果，不会涉及Cox模型。

下面是具体的实现步骤，我用示例数据来演示：

1. 加载必要的包 #

首先我们需要用到survival包拟合KM模型，tidyverse做数据整理和绘图，broom用来把survfit对象转换成整洁的数据框：

library(survival)
library(tidyverse)
library(broom)

2. 拟合Kaplan-Meier模型 #

这里我用R自带的lung数据集做示例，你可以换成自己的数据集。注意这里我们把连续预测变量（比如age）作为分层变量传入survfit()：

# 预处理示例数据：把status转换成0/1的删失变量（1=事件发生，0=删失）
df <- lung %>% mutate(censor = ifelse(status == 1, 1, 0))

# 拟合KM模型
km <- survfit(Surv(time, censor) ~ age, data = df)

3. 提取并整理KM数据 #

我们需要从survfit对象中提取每个分层（即每个连续变量值）的时间点和对应的生存概率，然后整理成适合绘图的格式：

# 把survfit对象转换成整洁数据框
km_tidy <- tidy(km) %>%
  # 从stratum列中提取连续预测变量的数值（比如这里的age）
  mutate(continuous_predictor = as.numeric(str_remove(stratum, "age="))) %>%
  # 保留我们需要的核心列
  select(continuous_predictor, time, estimate)

4. 获取目标时间点的生存概率 #

KM曲线的时间点是事件实际发生的时间，不一定正好是你关心的1年、5年这类时间点，所以我们用插值来得到精确时间点的生存概率：

# 定义你关心的时间点，比如1年（365天）、5年（1825天）
target_times <- c(365, 1825)

# 对每个连续变量值，插值得到目标时间点的生存概率，并转换成长格式
km_target <- km_tidy %>%
  group_by(continuous_predictor) %>%
  summarise(
    # 对每个目标时间点做线性插值
    across(target_times, ~ approx(time, estimate, xout = .x)$y)
  ) %>%
  # 转换成长格式，方便按时间点分组绘图
  pivot_longer(
    cols = -continuous_predictor,
    names_to = "time_point",
    values_to = "survival_prob"
  ) %>%
  # 给时间点添加更友好的标签
  mutate(time_point = paste0(time_point, " 天"))

5. 绘制可视化图 #

现在我们可以用ggplot来绘制你想要的图，横轴是连续预测变量，纵轴是生存概率，颜色区分不同时间点：

ggplot(km_target, aes(x = continuous_predictor, y = survival_prob, color = time_point)) +
  # 绘制折线和点，也可以换成平滑曲线
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  # 添加标签
  labs(
    x = "你的连续预测变量名称",
    y = "生存概率 P(Survival)",
    color = "时间点",
    title = "生存概率随连续预测变量的变化"
  ) +
  # 用简洁的主题
  theme_minimal()

如果你的连续变量取值很多，折线图看起来太杂乱，也可以换成平滑曲线来展示趋势：

ggplot(km_target, aes(x = continuous_predictor, y = survival_prob, color = time_point)) +
  geom_smooth(method = "loess", se = TRUE, linewidth = 1) +
  labs(
    x = "你的连续预测变量名称",
    y = "生存概率 P(Survival)",
    color = "时间点",
    title = "生存概率随连续预测变量的变化（平滑趋势）"
  ) +
  theme_minimal()

注意事项 #

• 如果某个连续变量值对应的KM曲线在到达目标时间点前就已经全部发生事件，插值会得到NA，你可以用filter(!is.na(survival_prob))去掉这些行。
• 如果你的连续变量有极高的基数（比如上千个唯一值），可以考虑先对变量进行分箱，或者用更平滑的拟合方法（比如GAM）来展示趋势。

备注：内容来源于stack exchange，提问作者purpleblade98