如何用R语言的非参数百分位Bootstrap法计算1-RR的置信区间

一、数据预处理：筛选符合条件的匹配对 #

先提取接种第三剂疫苗7天后仍处于风险的1:1匹配对（假设数据框为df，包含match_id匹配对ID、vaccine_group分组标识、days_since_3rd_dose第三剂接种后天数、outcome结局变量、follow_up_days随访天数字段）：

library(dplyr)

# 筛选接种7天后仍在风险的匹配对，确保每个匹配对的两个个体都符合条件
df_included <- df %>%
  filter(days_since_3rd_dose >= 7) %>%
  group_by(match_id) %>%
  filter(n() == 2) %>%
  ungroup()

二、Kaplan-Meier累积发病率曲线与风险估计 #

使用survival包实现曲线绘制与风险值提取：

library(survival)
library(survminer)

# 构建生存对象，转换为累积发病率曲线（1-生存概率）
km_fit <- survfit(Surv(follow_up_days, outcome) ~ vaccine_group, data = df_included)

# 绘制累积发病率曲线
ggsurvplot(km_fit, 
           data = df_included,
           fun = function(x) 1 - x, # 生存曲线转累积发病率
           xlab = "随访天数",
           ylab = "累积发病率",
           legend.title = "分组",
           palette = c("#2E86AB", "#F24236"))

# 提取特定时间点的累积风险估计（示例：30、60、90天）
summary(km_fit, times = c(30, 60, 90))

三、非参数百分位Bootstrap计算1-RR的95%置信区间 #

自定义Bootstrap抽样逻辑，重复1000次：

set.seed(123) # 固定随机种子保证结果可复现
n_boot <- 1000
boot_ve <- numeric(n_boot)

for (i in 1:n_boot) {
  # 按匹配对有放回抽样
  sampled_matches <- sample(unique(df_included$match_id), 
                            size = length(unique(df_included$match_id)), 
                            replace = TRUE)
  boot_sample <- df_included %>% filter(match_id %in% sampled_matches)
  
  # 计算两组的累积发病风险
  vaccine_risk <- mean(boot_sample$outcome[boot_sample$vaccine_group == "vaccinated"])
  control_risk <- mean(boot_sample$outcome[boot_sample$vaccine_group == "control"])
  
  # 计算1-RR（疫苗有效性）
  boot_ve[i] <- 1 - (vaccine_risk / control_risk)
}

# 提取95%百分位置信区间
ve_ci <- quantile(boot_ve, c(0.025, 0.975))
cat("1-RR的95%Bootstrap置信区间：", round(ve_ci[1], 3), "-", round(ve_ci[2], 3), "\n")

四、敏感性分析：Poisson回归估计1-发病率比 #

通过Poisson模型拟合发病率比，调整随访时间偏移：

# 拟合Poisson回归，加入随访天数作为偏移项
poisson_fit <- glm(outcome ~ vaccine_group + offset(log(follow_up_days)),
                   data = df_included,
                   family = poisson(link = "log"))

# 计算1-发病率比（疫苗有效性）及95%置信区间
irr <- exp(coef(poisson_fit)["vaccine_groupvaccinated"])
ve_irr <- 1 - irr
irr_ci <- exp(confint(poisson_fit)["vaccine_groupvaccinated", ])
ve_irr_ci <- c(1 - irr_ci[2], 1 - irr_ci[1])

cat("敏感性分析：1-发病率比的95%置信区间：", round(ve_irr_ci[1], 3), "-", round(ve_irr_ci[2], 3), "\n")

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