R语言Weibull参数生存模型5年预测的95%置信区间求解

针对你用survival包拟合Weibull模型后无法计算生存概率95%置信区间的问题，这里提供三种可行方法：

方法1：基于Delta方法的标准误计算 #

survreg拟合的是对数时间的极值分布模型，我们可以先获取位置参数（mu）的预测值和标准误，再通过Delta方法推导生存概率的方差，进而得到置信区间。

首先修正你代码中的语法错误（补全括号），并获取预测的标准误：

library(survival)
model <- survreg(Surv(time, status) ~ 1, data = lung)
time_days <- c(365.25*1, 365.25*3, 365.25*5)

# 提取link尺度的预测值及标准误
pred <- predict(model, newdata = data.frame(time = time_days), type = "link", se.fit = TRUE)
md_scale <- pred$fit
se_md <- pred$se.fit

# 计算基础生存预测值
sigma <- model$scale
shape <- 1/sigma
scale_weibull <- exp(md_scale)
surv_pred <- 1 - pweibull(time_days, shape = shape, scale = scale_weibull)

接着用Delta方法计算生存概率的标准误，再得到95%置信区间：

# 计算导数项（Delta方法核心）
weibull_term <- (time_days / scale_weibull)^shape
dS_dmu <- surv_pred * weibull_term * shape
se_surv <- abs(dS_dmu) * se_md

# 计算95%置信区间，同时限制在[0,1]范围内
lower_ci <- pmax(surv_pred - qnorm(0.975) * se_surv, 0)
upper_ci <- pmin(surv_pred + qnorm(0.975) * se_surv, 1)

# 整理结果
result <- data.frame(
  随访年限 = c(1,3,5),
  生存概率 = surv_pred,
  95%CI下限 = lower_ci,
  95%CI上限 = upper_ci
)
print(result)

方法2：基于模型参数的置信区间 #

直接利用模型位置参数mu的置信区间，代入Weibull生存函数计算生存概率的上下限（此方法默认sigma固定，若需考虑其不确定性，可结合轮廓似然，一般场景下足够用）：

# 提取模型参数及协方差矩阵
mu_est <- coef(model)
mu_se <- sqrt(vcov(model)[1,1])

# 计算mu的95%置信区间
mu_lower <- mu_est - qnorm(0.975) * mu_se
mu_upper <- mu_est + qnorm(0.975) * mu_se

# 代入生存函数计算生存概率的置信区间
surv_lower <- 1 - pweibull(time_days, shape = shape, scale = exp(mu_lower))
surv_upper <- 1 - pweibull(time_days, shape = shape, scale = exp(mu_upper))

# 整理结果
param_result <- data.frame(
  随访年限 = c(1,3,5),
  生存概率 = surv_pred,
  95%CI下限 = surv_lower,
  95%CI上限 = surv_upper
)
print(param_result)

方法3：Bootstrap法（更稳健） #

通过重采样数据集多次拟合模型，用预测值的分位数作为置信区间，能同时考虑所有参数的不确定性，结果更稳健：

set.seed(123) # 设置随机种子保证结果可重复
n_boot <- 1000 # 重采样次数
boot_surv <- matrix(NA, nrow = n_boot, ncol = length(time_days))

for(i in 1:n_boot){
  # 重采样生成bootstrap数据集
  boot_data <- lung[sample(nrow(lung), replace = TRUE), ]
  # 拟合模型
  boot_model <- survreg(Surv(time, status) ~ 1, data = boot_data)
  # 计算当前bootstrap样本的生存预测
  boot_mu <- predict(boot_model, newdata = data.frame(time = time_days), type = "link")
  boot_scale <- exp(boot_mu)
  boot_surv[i,] <- 1 - pweibull(time_days, shape = 1/boot_model$scale, scale = boot_scale)
}

# 取2.5%和97.5%分位数作为95%置信区间
boot_lower <- apply(boot_surv, 2, quantile, 0.025)
boot_upper <- apply(boot_surv, 2, quantile, 0.975)

# 整理结果
boot_result <- data.frame(
  随访年限 = c(1,3,5),
  生存概率 = surv_pred,
  95%CI下限_bootstrap = boot_lower,
  95%CI上限_bootstrap = boot_upper
)
print(boot_result)

内容的提问来源于stack exchange，提问作者user51962