分层分析中如何对回归系数子集的p值进行校正？

两种方法结果不一致的核心原因是校正的对象范围错误：

• 方法1的问题在于：如果是在管道内直接对全量p.value执行p.adjust再用if_else筛选替换，相当于把截距项的p值也纳入了校正池，干扰了非截距项的校正结果；若按分组单独校正，每个组仅1个非截距项，校正后p值会和原p值一致，显然不符合你的需求。
• 方法2的思路是对的，但需要把校正后的p值合并回原数据集，才能保留截距项的完整信息。

我们需要先提取所有非截距项的p值，统一进行多重比较校正，再将校正结果映射回原数据集，同时保留截距项的原始估计值和未校正p值（或标记为无需校正）。

完整代码如下：

library(purrr)
library(dplyr, warn.conflicts = FALSE)
library(broom)
library(tidyr)

# 1. 拟合分层回归，得到原始系数结果
mtcars_fit <- mtcars %>%
    group_by(cyl) %>%
    nest() %>%
    mutate(
        model = map(data, ~ lm(mpg ~ wt, data = .)),
        coeff = map(model, tidy, conf.int = FALSE)
    ) %>%
    unnest(coeff) %>%
    select(-statistic)

# 2. 提取所有非截距项，计算校正p值
non_intercept_p <- mtcars_fit %>%
    filter(term != "(Intercept)") %>%
    mutate(p.adj = p.adjust(p.value, method = "fdr")) # 可替换为你需要的校正方法，比如"bonferroni"

# 3. 合并回原数据集，保留截距项信息
mtcars_final <- mtcars_fit %>%
    left_join(non_intercept_p %>% select(cyl, term, p.adj), by = c("cyl", "term")) %>%
    # 截距项的校正p值设为NA，也可保留原p值，按需调整
    mutate(p.adj = if_else(term == "(Intercept)", NA_real_, p.adj))

执行上述代码后，非截距项的校正p值会和方法2的结果一致：0.0918、0.0206、0.0206。这是因为我们正确地将所有需要校正的非截距项p值作为一个整体进行校正，符合多重比较校正的逻辑——仅针对同一类假设检验（所有组的wt系数显著性检验）控制一类错误，排除了无关的截距项检验干扰。

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