请求协助计算中介分析中效应指标的Bootstrap 95%置信区间

要计算总效应(TE)、直接效应(DE)、间接效应(IE)及中介比例(PM)的Bootstrap 95%置信区间，可按以下步骤实现：

1. 加载所需R包 #

先确保安装并加载拟合模型和抽样所需的包：

# 安装包（若未安装）
install.packages(c("survival", "boot", "rms"))

# 加载包
library(survival)
library(boot)
library(rms)

2. 定义Bootstrap统计量函数 #

编写自定义函数，用于每次抽样后拟合模型并计算目标指标：

boot_cox_mediation <- function(data, indices) {
  # 获取Bootstrap抽样样本
  boot_data <- data[indices, ]
  
  # 拟合加权Cox模型
  fit <- coxph(Surv(time, event) ~ const(factor(exp)) + const(factor(expstar)) + covariates,
               data = boot_data, weights = boot_data$weightM)
  
  # 提取系数并计算各效应指标
  coefs <- coef(fit)
  te <- exp(sum(coefs[c('const(factor(exp))1', 'const(factor(expstar))1')]))
  de <- exp(unname(coefs['const(factor(exp))1']))
  ie <- exp(coefs['const(factor(expstar))1'])
  pm <- log(ie)/log(te)
  
  # 返回四个指标结果
  return(c(TE = te, DE = de, IE = ie, PM = pm))
}

3. 执行Bootstrap抽样 #

设置随机种子保证结果可重复，执行抽样（建议至少1000次）：

set.seed(123)
boot_results <- boot(data = newmyd, statistic = boot_cox_mediation, R = 1000)

4. 计算95%置信区间 #

用百分位数法提取各指标的置信区间：

# 提取每个指标的置信区间
boot_ci <- lapply(1:4, function(i) {
  boot.ci(boot_results, type = "perc", index = i)
})

# 命名并输出结果
names(boot_ci) <- c("TE", "DE", "IE", "PM")
print(boot_ci)

注意事项 #

• 若抽样过程中出现模型收敛问题，可适当增加抽样次数，或检查原始数据的变量分布与样本量。
• const()函数来自rms包，必须加载该包才能正常拟合模型。

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