拟合生存曲线的Weibull分布后，如何正确估算生存值并绘图？

你遇到的问题核心是用错了拟合工具的场景——fitdist函数是用来拟合随机变量的观测值的分布（比如你有一批失效时间数据，用它来拟合这些时间的Weibull分布），但你直接把生存概率survival的值传给它，相当于让它去拟合“生存概率”这个变量的分布，这和你想要的“拟合生存曲线对应的Weibull模型”完全是两回事，所以得到的参数和计算结果肯定不对。

下面给你两种正确的拟合思路和实现代码：

方法一：用生存分析专用工具拟合（推荐） #

我们可以把现有的月度生存概率数据，转换为生存分析标准的“风险集-事件数”格式，然后用flexsurv包的flexsurvreg函数（专门为生存模型设计）来拟合Weibull分布。

代码实现 #

library(flexsurv)
library(tidyverse)

# 加载你的数据
data <- tibble(month = 0:52, survival = c(1, 1, 1, 1, 1, 1, 0.999, 0.998, 0.997, 0.993, 0.984, 0.976, 0.973, 0.971, 0.969, 0.969, 0.969, 0.969, 0.968, 0.968, 0.968, 0.968, 0.968, 0.968, 0.967, 0.966, 0.966, 0.962, 0.957, 0.952, 0.948, 0.944, 0.942, 0.941, 0.941, 0.941, 0.941, 0.940, 0.939, 0.938, 0.938, 0.938, 0.935, 0.934, 0.930, 0.920, 0.910, 0.895, 0.884, 0.881, 0.879, 0.878, 0.878))

# 转换数据：计算每个时间点的失效概率、事件数（假设初始风险集为1000，比例不影响拟合结果）
data <- data %>%
  mutate(
    prev_survival = lag(survival, default = 1),
    fail_prob = 1 - survival / prev_survival,
    n_risk = 1000,
    n_events = round(n_risk * fail_prob),
    # 计算删失人数（风险集减去事件数再减去下一个时间点的风险集）
    n_censor = n_risk - n_events - lead(n_risk, default = 0)
  )

# 拟合Weibull生存模型
fit_weibull <- flexsurvreg(
  formula = Surv(month, event = (n_events > 0)) ~ 1,
  data = data,
  dist = "weibull"
)

# 查看拟合得到的参数
print(fit_weibull)

# 计算拟合的生存值（Weibull生存函数：S(t) = exp(-(t/scale)^shape)）
shape <- fit_weibull$coefficients["shape"]
scale <- fit_weibull$coefficients["scale"]
data <- data %>%
  mutate(fitted_survival = exp(-(month / scale)^shape))

# 绘制原始曲线和拟合曲线
ggplot(data, aes(x = month)) +
  geom_line(aes(y = survival, color = "原始生存曲线"), linewidth = 1) +
  geom_line(aes(y = fitted_survival, color = "拟合Weibull曲线"), linewidth = 1, linetype = "dashed") +
  labs(x = "月份", y = "生存概率", color = "曲线类型") +
  theme_minimal()

方法二：直接用非线性最小二乘拟合生存函数 #

如果你只是想找到Weibull生存函数的参数，让它尽可能贴合你的观测生存值，可以直接用非线性最小二乘（nls）拟合Weibull生存函数的表达式，这种方法更直接。

代码实现 #

library(tidyverse)

# 加载数据（和上面一致）
data <- tibble(month = 0:52, survival = c(1, 1, 1, 1, 1, 1, 0.999, 0.998, 0.997, 0.993, 0.984, 0.976, 0.973, 0.971, 0.969, 0.969, 0.969, 0.969, 0.968, 0.968, 0.968, 0.968, 0.968, 0.968, 0.967, 0.966, 0.966, 0.962, 0.957, 0.952, 0.948, 0.944, 0.942, 0.941, 0.941, 0.941, 0.941, 0.940, 0.939, 0.938, 0.938, 0.938, 0.935, 0.934, 0.930, 0.920, 0.910, 0.895, 0.884, 0.881, 0.879, 0.878, 0.878))

# 拟合Weibull生存函数：S(t) = exp(-(t/scale)^shape)
nls_fit <- nls(
  survival ~ exp(-(month/scale)^shape),
  data = data,
  start = list(shape = 1, scale = 20) # 初始参数值，可根据调整
)

# 查看拟合结果
summary(nls_fit)

# 计算拟合的生存值
data <- data %>%
  mutate(fitted_survival = predict(nls_fit))

# 绘图对比
ggplot(data, aes(x = month)) +
  geom_line(aes(y = survival, color = "原始生存曲线"), linewidth = 1) +
  geom_line(aes(y = fitted_survival, color = "拟合Weibull曲线"), linewidth = 1, linetype = "dashed") +
  labs(x = "月份", y = "生存概率", color = "曲线类型") +
  theme_minimal()

为什么之前的方法错误？ #

你之前用fitdist(data[['survival']], 'weibull')，相当于告诉R：“这里有53个独立的观测值，它们服从Weibull分布，请拟合这个分布的参数”。但你的survival是不同时间点的生存概率，不是独立的随机变量观测值，所以拟合出来的参数完全不符合生存曲线的需求，自然计算出的生存值也完全不对。

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