如何在R中实现带Bonferroni校正p值的分组多元线性回归？

问题描述 #

用户提供样本数据：

dad <- data.frame(type = c("new", "new", "old", "new", "old", "new", "old", "old", "new", "new"),
                   outcome = c(68,76,57,67,89,98,99,120,99,67),
                  loan1=c(98000 ,56668,67000,87999,87945, 89768, 65738, 87547, 78937, 20000),
                  loan2 =c(56768,45679, 23453, 87673, 56749, 45783, 67836, 54673, 45379, 78483),
                  house=c("semi", "detached", "town", "semi", "town", "town", "semi", "detached", "semi", "town"))

需求：针对每种house类型分别运行多元线性回归模型（实际场景为7组，初始自变量相同，通过向后逐步选择得到最终模型），对所有模型的p值进行Bonferroni校正，询问现有代码是否能实现，以及正确的实现方法。

用户现有示例代码：

#Example 
full <- lm(outcome ~loan1+loan2,data = subset(dad, house=="semi"))
summary(full)
tab_model(full)

解决方案 #

你提供的现有代码只能单独运行某一组的回归、查看未校正的结果，无法直接获取Bonferroni校正后的p值。下面是完整的实现步骤：

1. 分组运行向后逐步回归模型 #

先按house分组，对每组数据运行向后逐步选择的回归模型：

# 加载必要包
library(MASS)
library(broom)

# 按house分组，批量运行向后逐步回归
group_models <- lapply(split(dad, dad$house), function(df) {
  # 初始全模型
  full_model <- lm(outcome ~ loan1 + loan2, data = df)
  # 向后逐步选择（trace=FALSE关闭过程输出）
  step_model <- stepAIC(full_model, direction = "backward", trace = FALSE)
  return(step_model)
})

2. 提取所有原始p值 #

从每个模型中提取自变量的p值，统一收集：

# 提取所有模型的自变量p值（排除截距项）
all_pvals <- lapply(group_models, function(model) {
  tidy(model) %>% 
    filter(term != "(Intercept)") %>% 
    pull(p.value)
})

# 将所有p值合并为一个向量，用于后续校正
flat_pvals <- unlist(all_pvals)

3. 执行Bonferroni校正 #

Bonferroni校正的核心是将原始p值乘以总检验次数（即所有模型中自变量的总数）：

# 计算总检验次数
total_tests <- length(flat_pvals)
# 执行Bonferroni校正
corrected_pvals <- p.adjust(flat_pvals, method = "bonferroni")

4. 整合校正结果并展示 #

将校正后的p值对应回原模型的自变量，生成结构化结果：

# 把校正p值匹配到对应模型的自变量上
result_list <- mapply(function(model, p_corr) {
  tidy(model) %>% 
    filter(term != "(Intercept)") %>% 
    mutate(corrected_p.value = p_corr) %>% 
    select(term, estimate, std.error, statistic, p.value, corrected_p.value)
}, group_models, split(corrected_pvals, rep(names(group_models), sapply(all_pvals, length))), SIMPLIFY = FALSE)

# 查看任意组的结果，比如semi组
print(result_list[["semi"]])

5. 用tab_model展示校正结果（可选） #

如果想用tab_model展示带校正p值的结果，可手动添加校正列后输出：

library(sjPlot)

# 给每个模型添加校正p值属性
for(i in seq_along(group_models)) {
  model <- group_models[[i]]
  corr_p <- corrected_pvals[grepl(names(group_models)[i], names(corrected_pvals))]
  attr(model, "corrected_p") <- corr_p
}

# 自定义输出，包含校正p值
tab_model(group_models, 
          show.p = TRUE,
          show.stat = TRUE,
          add.lines = list(
            "校正p值" = sapply(group_models, function(x) paste(round(attr(x, "corrected_p"), 4), collapse = "; "))
          ))

关键说明 #

• Bonferroni校正的检验次数是所有模型中需要检验的自变量总数，不是模型的组数。比如你有7个模型，每个最终保留2个自变量，总检验次数就是14次。
• 如果每组逐步回归保留的自变量数量不同，上述代码会自动适配，无需手动调整。

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