如何合并R中adjustedCurves包生成的50个校正生存风险表？

一、提取adjustedsurv()中的风险表 #

首先要确保调用adjustedsurv()时开启风险表返回，修改你的分析代码：

library(mice)
library(survival)
library(adjustedCurves)
library(data.table)
library(gt)

# 假设imp是你生成的多重插补对象
adj_surv_list <- lapply(1:50, function(m) {
  dat <- complete(imp, m)
  # 计算IPtW权重
  ps_model <- glm(trt ~ cov1 + cov2 + cov3, data = dat, family = binomial)
  dat$iptw <- weights(ps_model, type = "ipw")
  # 拟合加权Cox模型
  cox_model <- coxph(Surv(time, event) ~ trt, data = dat, weights = iptw)
  # 生成校正曲线时，必须设置return_risktable=TRUE
  adj_surv <- adjustedsurv(data = dat,
                           variable = "trt",
                           ev_time = "time",
                           event = "event",
                           method = "iptw",
                           outcome_model = cox_model,
                           return_risktable = TRUE)
  adj_surv
})

# 提取所有插补数据集的风险表
risktables_list <- lapply(adj_surv_list, function(x) {
  rt <- x$risktable
  rt$imp_num <- m # 标记插补编号
  rt
})

# 合并为一个大数据表
all_risktables <- rbindlist(risktables_list)

二、用Rubin法则合并风险表 #

针对风险表中的n_at_risk（风险人数）和n_events（事件数），按时间点、治疗组分组合并：

m <- 50 # 插补数据集数量

merged_risktable <- all_risktables[, .(
  # 点估计取各插补数据集的均值
  n_at_risk = mean(n_at_risk),
  n_events = mean(n_events),
  # 计算标准误与95%置信区间（基于t分布）
  se_at_risk = sd(n_at_risk)/sqrt(m),
  se_events = sd(n_events)/sqrt(m),
  ci_lower_at_risk = n_at_risk - qt(0.975, df = m-1)*se_at_risk,
  ci_upper_at_risk = n_at_risk + qt(0.975, df = m-1)*se_at_risk,
  ci_lower_events = n_events - qt(0.975, df = m-1)*se_events,
  ci_upper_events = n_events + qt(0.975, df = m-1)*se_events
), by = .(time, strata)]

# 格式化输出列
merged_risktable <- merged_risktable[, .(
  time = round(time, 1),
  treatment = strata,
  n_at_risk = round(n_at_risk),
  n_at_risk_ci = paste0("[", round(ci_lower_at_risk), ", ", round(ci_upper_at_risk), "]"),
  n_events = round(n_events),
  n_events_ci = paste0("[", round(ci_lower_events), ", ", round(ci_upper_events), "]")
)]

三、用data.table自行计算风险表 #

如果需要手动计算，基于校正生存概率推导：

# 提取所有插补数据集的校正生存数据
adj_data_list <- lapply(adj_surv_list, function(x) x$adj)
all_adj_data <- rbindlist(adj_data_list, idcol = "imp_num")

# 计算每个插补数据集的初始加权样本量
initial_n <- all_adj_data[, .(initial_weighted_n = sum(weights)), by = .(imp_num, strata)]
all_adj_data <- merge(all_adj_data, initial_n, by = c("imp_num", "strata"))

# 推导风险人数与事件数
all_adj_data <- all_adj_data[order(imp_num, strata, time)]
all_adj_data[, surv_prev := shift(surv, fill = 1), by = .(imp_num, strata)]
all_adj_data[, n_at_risk := initial_weighted_n * surv]
all_adj_data[, n_events := initial_weighted_n * (surv_prev - surv)]
all_adj_data[, cum_events := cumsum(n_events), by = .(imp_num, strata)]

# 后续合并步骤同第二部分

四、用gt包可视化合并后的风险表 #

gt可以轻松生成美观的分层风险表：

gt_table <- merged_risktable %>%
  gt(groupname_col = "treatment") %>%
  tab_header(title = "校正后风险表（多重插补合并）",
             subtitle = "基于IPtW Cox模型，50个插补数据集合并结果") %>%
  cols_label(
    time = "时间",
    n_at_risk = "风险人数",
    n_at_risk_ci = "95%置信区间",
    n_events = "事件数",
    n_events_ci = "95%置信区间"
  ) %>%
  fmt_number(columns = c(n_at_risk, n_events), decimals = 0) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_row_groups()
  ) %>%
  tab_options(
    table.width = pct(90),
    row_group.border.top = border(width = 2),
    row_group.border.bottom = border(width = 1)
  )

# 输出表格
gt_table

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