生存分析预测：如何优化R代码生成合理多模拟路径？

问题背景 #

基于survival包的lung数据集，构建仅包含前500期数据的lung1子集，拟合Weibull参数模型（已得到shape=1.804891、scale=306.320693），需要生成501-1000期的单调递减模拟生存路径，并与实际501-1000期的生存数据对比，完成假设情景分析。

核心问题分析 #

原代码无法生成合理单调路径，通常是因为未正确应用条件生存概率逻辑：Weibull分布的生存概率随时间单调递减，但直接生成独立生存时间可能忽略“从t=500开始的条件生存”约束，导致模拟路径出现波动甚至上升。

优化代码方案 #

library(survival)
library(ggplot2)

# 1. 数据预处理
# 构建lung1：仅保留时间<=500的样本
lung1 <- subset(lung, time <= 500)
# 全数据集用于后续对比实际后期生存曲线
full_lung <- lung

# 2. 拟合Weibull模型（验证用户提供的参数）
weibull_fit <- survreg(Surv(time, status) ~ 1, data = lung1, dist = "weibull")
# 提取Weibull参数：survreg的scale是1/shape，location是log(scale)
shape <- 1/weibull_fit$scale
scale <- exp(weibull_fit$coefficients[1])
cat("拟合得到的Weibull参数：shape =", round(shape, 6), "scale =", round(scale, 6), "\n")

# 3. 计算lung1的K-M生存曲线（t=1到500）
km_lung1 <- survfit(Surv(time, status) ~ 1, data = lung1)
# 提取t=500时的生存概率（作为模拟的起始点）
surv_500 <- stepfun(km_lung1$time, c(1, km_lung1$surv))(500)

# 4. 模拟501-1000期的生存路径
n_sim <- 100  # 模拟路径数量
time_pred <- 501:1000

# 初始化存储模拟生存概率的矩阵
sim_surv <- matrix(NA, nrow = n_sim, ncol = length(time_pred))

set.seed(123)  # 固定随机种子保证可复现
for (i in 1:n_sim) {
  # 基于Weibull条件生存公式计算每个时间点的生存概率
  # S(t | t0) = S(t)/S(t0)，确保从t=500的生存概率开始单调递减
  sim_surv[i, ] <- surv_500 * (1 - pweibull(time_pred, shape, scale)) / (1 - pweibull(500, shape, scale))
}

# 5. 计算实际全数据集501-1000期的K-M曲线
km_full <- survfit(Surv(time, status) ~ 1, data = full_lung)
# 提取实际后期生存概率
actual_surv <- stepfun(km_full$time, c(1, km_full$surv))(time_pred)

# 6. 可视化
# 整理模拟数据为数据框
sim_df <- data.frame(
  time = rep(time_pred, n_sim),
  surv = as.vector(sim_surv),
  sim_id = factor(rep(1:n_sim, each = length(time_pred)))
)

# 整理lung1的K-M数据到t=500
lung1_km_df <- data.frame(
  time = c(km_lung1$time[km_lung1$time <= 500], 500),
  surv = c(km_lung1$surv[km_lung1$time <= 500], surv_500)
) %>% unique()  # 去重避免重复点

# 绘图
ggplot() +
  # 绘制lung1的K-M曲线（1-500期）
  geom_step(data = lung1_km_df, aes(x = time, y = surv), color = "blue", size = 1) +
  # 绘制模拟路径（501-1000期）
  geom_line(data = sim_df, aes(x = time, y = surv, group = sim_id), color = "gray", alpha = 0.3) +
  # 绘制模拟路径的均值
  geom_line(data = data.frame(time = time_pred, surv = colMeans(sim_surv)),
            aes(x = time, y = surv), color = "red", size = 1) +
  # 绘制实际后期K-M曲线（501-1000期）
  geom_step(data = data.frame(time = time_pred, surv = actual_surv),
            aes(x = time, y = surv), color = "green", size = 1, linetype = "dashed") +
  labs(x = "时间", y = "生存概率", title = "肺癌生存概率模拟与实际对比") +
  scale_y_continuous(limits = c(0, 1)) +
  theme_minimal()

关键优化点 #

• 条件生存概率计算：模拟时基于t=500的生存概率，使用Weibull的条件生存公式S(t | t0) = S(t)/S(t0)（t>t0），严格保证生存概率单调递减。
• 可复现性：设置随机种子set.seed()，确保模拟结果可重复。
• 可视化分层：区分展示历史K-M曲线、多模拟路径、模拟均值路径和实际后期曲线，清晰对比假设情景与真实情况。
• 参数验证：代码中包含Weibull模型拟合步骤，可直接验证用户提供的参数是否与lung1数据匹配。

内容的提问来源于stack exchange，提问作者Village.Idyot