如何在R的rpart生存树中计算估计率的置信区间？

背景说明 #

rpart生存树中每个节点的相对死亡率是节点风险率与整体样本风险率的比值，rpart本身不提供该值的置信区间。由于Cox模型的HR是组间对比（以参考组为基准），与rpart的相对死亡率定义不同，需要通过以下方法计算置信区间。

方法一：基于Cox截距模型的近似CI（简单但存在样本重叠偏倚） #

该方法通过拟合无协变量的Cox模型，分别估计整体和各节点的风险率，再计算相对比值的CI。但因节点样本是整体的子集，估计量不独立，CI可能存在偏倚。

代码实现 #

library(rpart)
library(survival)
library(rpart.plot)
library(dplyr)

# 加载数据并构建剪枝后的生存树
data(stagec)
pfit <- rpart(Surv(pgtime, pgstat) ~ age + eet + g2 + grade + gleason + ploidy, data = stagec)
pfit2 <- prune(pfit, cp = 0.016)

# 计算整体样本的风险率（无协变量Cox模型）
cox_overall <- coxph(Surv(pgtime, pgstat) ~ 1, data = stagec)
overall_risk <- exp(coef(cox_overall))

# 遍历每个节点计算相对死亡率及CI
node_ids <- unique(pfit2$where)
node_results <- list()

for (node in node_ids) {
  node_data <- stagec[pfit2$where == node, ]
  # 节点的无协变量Cox模型
  cox_node <- coxph(Surv(pgtime, pgstat) ~ 1, data = node_data)
  node_risk <- exp(coef(cox_node))
  # 相对死亡率 = 节点风险率 / 整体风险率
  rel_risk <- node_risk / overall_risk
  # 转换节点风险率的CI为相对死亡率的CI
  node_risk_ci <- exp(confint(cox_node))
  rel_risk_ci <- node_risk_ci / overall_risk
  
  node_results[[as.character(node)]] <- list(
    相对死亡率 = rel_risk,
    95%CI下限 = rel_risk_ci[1],
    95%CI上限 = rel_risk_ci[2],
    节点样本量 = nrow(node_data)
  )
}

# 输出结果
node_df <- bind_rows(node_results, .id = "节点编号")
print(node_df)

方法二：Bootstrap方法（推荐，更可靠） #

通过重复抽样构建生存树，模拟树结构和风险估计的变异，从而得到相对死亡率的置信区间，能更好反映真实的统计变异。

代码实现 #

set.seed(123) # 固定种子保证可重复性
n_boot <- 1000 # 抽样次数，次数越多结果越稳定

# 存储bootstrap结果
boot_rel_risks <- list()

for (b in 1:n_boot) {
  # 有放回抽样生成bootstrap数据集
  boot_data <- stagec[sample(nrow(stagec), replace = TRUE), ]
  # 构建并剪枝生存树
  boot_pfit <- rpart(Surv(pgtime, pgstat) ~ age + eet + g2 + grade + gleason + ploidy, data = boot_data)
  boot_pfit2 <- prune(boot_pfit, cp = 0.016)
  
  # 提取每个节点的相对死亡率（rpart的predict(type="risk")返回相对于根节点的风险）
  boot_node_data <- data.frame(
    节点编号 = boot_pfit2$where,
    相对死亡率 = predict(boot_pfit2, type = "risk")
  ) %>%
    group_by(节点编号) %>%
    summarise(平均相对死亡率 = mean(相对死亡率), .groups = "drop") %>%
    mutate(bootstrap_id = b)
  
  boot_rel_risks[[b]] <- boot_node_data
}

# 合并bootstrap结果并计算CI
boot_df <- bind_rows(boot_rel_risks)
node_ci <- boot_df %>%
  group_by(节点编号) %>%
  summarise(
    中位相对死亡率 = median(平均相对死亡率),
    95%CI下限 = quantile(平均相对死亡率, 0.025),
    95%CI上限 = quantile(平均相对死亡率, 0.975),
    有效bootstrap次数 = n()
  ) %>%
  arrange(节点编号)

print(node_ci)

# 将CI添加到生存树可视化中
node_labels <- node_ci %>%
  mutate(标签 = sprintf("%.2f (%.2f-%.2f)", 中位相对死亡率, 95%CI下限, 95%CI上限)) %>%
  pull(标签, name = 节点编号)

rpart.plot(pfit2, extra = 100,
           node.fun = function(x, labs, digits, varlen) {
             node_labels[as.character(x$node)]
           })

关键说明 #

• rpart的相对死亡率定义为节点风险率与整体样本风险率的比值，与Cox模型中基于参考组的HR不同，不能直接用Cox的HR CI替代。
• Bootstrap方法考虑了树构建过程中的随机变异，是更可靠的CI估计方式，尤其适合小样本节点。
• 若节点样本量极小，可适当增加bootstrap抽样次数（如1000次）以提高CI稳定性。

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