从多变量GLM输出中提取各水平样本数的技术问询

针对多变量二项式GLM模型无法提取各协变量水平样本数、合并结果表的问题，这里提供一套基于broom和tidyverse的简便解决方案，无需手动处理模型矩阵，直接从原始数据和模型结果中整合所需信息。

步骤1：准备环境与模拟数据 #

先加载所需包，并模拟包含多分类协变量的示例数据（贴合你的gender+occupation场景）：

library(tidyverse)
library(broom)

# 模拟示例数据
set.seed(123)
n <- 500
dat <- tibble(
  gender = sample(c("Male", "Female"), n, replace = TRUE),
  occupation = sample(c("Student", "Worker", "Retired"), n, replace = TRUE, prob = c(0.4, 0.5, 0.1)),
  outcome = rbinom(n, 1, plogis(-1 + 0.5*(gender=="Female") + 0.8*(occupation=="Worker") + 1.2*(occupation=="Retired")))
)

# 拟合多变量二项式GLM
model <- glm(outcome ~ gender + occupation, data = dat, family = binomial)

步骤2：提取各协变量水平的样本数 #

从模型对应的原始数据中统计每个协变量各水平的样本量，避免直接操作模型矩阵的麻烦：

sample_counts <- dat %>%
  select(gender, occupation) %>% # 选择所有协变量
  pivot_longer(cols = everything(), names_to = "Covariate", values_to = "Level") %>%
  count(Covariate, Level, name = "Sample_Size")

步骤3：提取模型的系数、P值与置信区间 #

用broom::tidy()一键提取模型结果，并通过字符串处理拆分协变量与水平：

model_results <- tidy(model, conf.int = TRUE) %>%
  mutate(
    # 拆分协变量名称
    Covariate = case_when(
      term == "(Intercept)" ~ "(Intercept)",
      TRUE ~ str_remove(term, "^.*(?=\\()") %>% str_remove("\\(|\\)")
    ),
    # 拆分水平名称
    Level = case_when(
      term == "(Intercept)" ~ "Reference",
      TRUE ~ str_extract(term, "(?<=\\().*(?=\\))")
    )
  ) %>%
  select(Covariate, Level, Beta_Estimate = estimate, P_val = p.value, CI_Lower = conf.low, CI_Upper = conf.high)

步骤4：合并生成最终结果表 #

将样本数与模型结果合并，补充截距项的总样本量：

final_table <- model_results %>%
  left_join(sample_counts, by = c("Covariate", "Level")) %>%
  mutate(Sample_Size = ifelse(Covariate == "(Intercept)", n, Sample_Size)) %>%
  select(Covariate, Level, Sample_Size, Beta_Estimate, P_val, CI_Lower, CI_Upper)

# 查看结果
print(final_table, digits = 3)

最终输出样式（示例）：

# A tibble: 4 × 7
  Covariate    Level     Sample_Size Beta_Estimate P_val CI_Lower CI_Upper
  <chr>        <chr>           <int>         <dbl> <dbl>    <dbl>    <dbl>
1 (Intercept)  Reference         500         -1.36 0.001    -2.17    -0.55
2 gender       Female            252          0.49 0.129    -0.13     1.11
3 occupation   Worker            257          0.79 0.003     0.22     1.36
4 occupation   Retired            47          1.41 0.019     0.23     2.59

循环处理多个模型的适配 #

如果需要批量处理多个模型，可将上述步骤封装为函数，用purrr::map()批量执行：

process_glm <- function(model) {
  dat <- model.frame(model)
  covars <- names(dat)[-1]
  
  sample_counts <- dat %>%
    select(all_of(covars)) %>%
    pivot_longer(cols = everything(), names_to = "Covariate", values_to = "Level") %>%
    count(Covariate, Level, name = "Sample_Size")
  
  model_results <- tidy(model, conf.int = TRUE) %>%
    mutate(
      Covariate = case_when(
        term == "(Intercept)" ~ "(Intercept)",
        TRUE ~ str_remove(term, "^.*(?=\\()") %>% str_remove("\\(|\\)")
      ),
      Level = case_when(
        term == "(Intercept)" ~ "Reference",
        TRUE ~ str_extract(term, "(?<=\\().*(?=\\))")
      )
    ) %>%
    select(Covariate, Level, Beta_Estimate = estimate, P_val = p.value, CI_Lower = conf.low, CI_Upper = conf.high)
  
  final_table <- model_results %>%
    left_join(sample_counts, by = c("Covariate", "Level")) %>%
    mutate(Sample_Size = ifelse(Covariate == "(Intercept)", nrow(dat), Sample_Size)) %>%
    select(Covariate, Level, Sample_Size, Beta_Estimate, P_val, CI_Lower, CI_Upper)
  
  return(final_table)
}

# 假设output是存储多个GLM模型的列表，批量处理
# map(output, process_glm)

核心思路说明 #

• 放弃直接操作模型矩阵（output[[i]]$model）的方式，转而从模型绑定的原始数据中统计样本量，避免多变量场景下模型矩阵的哑变量干扰
• 用broom包标准化提取模型结果，减少手动提取系数、P值的代码量
• 通过字符串处理自动拆分协变量与水平，适配二分类、多分类协变量的统一处理

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