使用bshazard包绘制风险函数、获取生命表风险比及置信区间问题

问题描述 #

我在博士阶段的生存分析研究中需要绘制风险函数，对比两种不同条件下的风险率（hazard rate）。
我参考了相关文献第7页图4的实现方案，想要获取预测变量两个水平对应的平滑风险线置信区间，但代码无法达成预期效果。
我的参考代码如下：

fitt<-bshazard(Surv(time,event) ~ session.type,data=data,lambda=10,nbin=60) 
plot(fitt,overall=FALSE, col=1, conf.int = TRUE)

设置overall=FALSE参数后我得到了两条平滑风险曲线，但未显示我用于推导绘图结果所需的置信区间。
如果有方法可以导出包含各时间区间风险率及上下置信区间的时间表格，也可满足我的需求。

解决方案 #

1. 手动绘制带置信区间的分组风险曲线

bshazard包的默认绘图方法在overall=FALSE时不会自动渲染分组的置信区间，你可以直接提取拟合结果里的分层数据手动绘图，示例代码如下：

# 加载所需包
library(bshazard)
library(ggplot2)

# 原有拟合逻辑保留
fitt <- bshazard(Surv(time, event) ~ session.type, data = data, lambda = 10, nbin = 60)

# 提取两个分组的拟合数据，session.type的两个水平对应strata列表的两个元素
group1_data <- data.frame(
  time = fitt$strata[[1]]$time,
  hazard = fitt$strata[[1]]$hazard,
  lower_ci = fitt$strata[[1]]$lower,
  upper_ci = fitt$strata[[1]]$upper,
  group = names(fitt$strata)[1]
)
group2_data <- data.frame(
  time = fitt$strata[[2]]$time,
  hazard = fitt$strata[[2]]$hazard,
  lower_ci = fitt$strata[[2]]$lower,
  upper_ci = fitt$strata[[2]]$upper,
  group = names(fitt$strata)[2]
)
plot_data <- rbind(group1_data, group2_data)

# 绘制带置信区间的风险曲线
ggplot(plot_data, aes(x = time, color = group, fill = group)) +
  geom_line(aes(y = hazard), linewidth = 1) +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), alpha = 0.2, color = NA) +
  labs(x = "时间", y = "风险率 (Hazard Rate)", color = "分组", fill = "分组") +
  theme_bw()

2. 导出风险率与置信区间表格

直接将整理好的全量数据导出为csv文件即可：

write.csv(plot_data, "分组风险率置信区间明细.csv", row.names = FALSE, fileEncoding = "UTF-8")

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