R中GLM相对风险、风险差多重检验的Bonferroni校正方法咨询

1. 适配RR/RD场景的多重检验校正R包说明 #

不存在专门仅适配RR/RD效应量的多重检验校正R包，因为多重检验校正的核心是调整检验水准/置信区间覆盖度，和效应量类型无关。目前常用的通用R包都可以支持你的场景：

• emmeans：可先计算GLM模型的边际预测值，得到RR/RD估计值后直接调用内置的adjust参数完成多种校正（含Bonferroni、BH法等）
• multcomp：支持自定义效应量对比，可直接输出校正后的p值和置信区间
• 基础R的p.adjust函数也可以直接对批量计算得到的原始p值做校正，适配所有效应量场景。

2. Bonferroni校正拓宽bootstrap置信区间的实现方法 #

你目前用的是自助法（bootstrap）百分位置信区间，Bonferroni校正的逻辑非常简单：
你总共有7个结局指标 × 3个随访时间点 = 21次独立检验，原95%置信区间对应的一类错误概率α=0.05，Bonferroni校正需要将α除以总检验次数得到校正后的检验水准：
校正后α = 0.05 / 总检验次数m
对应的置信区间分位数就从原来的c(0.025, 0.975)替换为c(α/2, 1-α/2)即可。

代码修改示例 #

假设你总共有21次检验，修改你原来结果汇总的代码段即可：

# 定义总检验次数
n_test <- 7 * 3
# 计算校正后的分位点
alpha_adj <- 0.05 / n_test
lower_q <- alpha_adj / 2
upper_q <- 1 - alpha_adj / 2

# 生成校正后的结果
m_adj = rbind(
  c(mean(rr), quantile(rr, c(lower_q, upper_q)), twosidep(rr, test = 1) ),
  c(mean(rd), quantile(rd, c(lower_q, upper_q)), twosidep(rd) )
)
rownames(m_adj) = c('rr', 'rd')
colnames(m_adj) = c('Estimate', 'Adjusted Lower CI', 'Adjusted Upper CI', 'Raw p-value')
# 如果需要输出校正后的p值，收集所有检验的原始p值后做如下处理
# adjusted_p <- p.adjust(all_raw_p, method = "bonferroni")
m_adj

注意事项 #

• Bonferroni校正属于比较保守的校正方法，若你认为检验次数多导致统计效力太低，可替换为FDR校正（Benjamini-Hochberg法），p值校正只需把p.adjust的method参数改为"BH"即可。
• 若不同结局/时间点的检验存在相关性，也可以考虑使用置换检验（Permutation test）做校正，会比Bonferroni更高效。

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