如何为数据中不同地理区域快速批量运行多个glm模型？

错误原因说明 #

你之前的代码报错有两个核心原因：

1. subset参数用法错误：subset需要传入逻辑判断条件（比如region == "north"），直接传入变量名会被解析为数值向量，无法正确筛选目标分组数据。
2. 分组内因子水平不足：报错提示“对比只能应用于有2个及以上水平的因子”，说明你的数据中存在某个region分组，其中grade变量只有1种取值（比如某区域全是"A"），这种情况下GLM无法拟合因子的对比项。

预处理：过滤无效分组 #

先过滤掉grade水平不足2个的区域（如果需要，还可以额外过滤因变量无变异的分组）：

library(tidyverse)

# 过滤grade水平≥2的区域
valid_regions <- df %>%
  group_by(region) %>%
  summarise(n_grade = n_distinct(grade)) %>%
  filter(n_grade >= 2) %>%
  pull(region)

df_filtered <- df %>% filter(region %in% valid_regions)

# （可选）额外过滤因变量无变异的分组（以var1为例）
valid_regions_full <- df_filtered %>%
  group_by(region) %>%
  summarise(var_var1 = var(var1, na.rm = TRUE)) %>%
  filter(var_var1 > 0) %>%
  pull(region)

df_filtered <- df_filtered %>% filter(region %in% valid_regions_full)

批量拟合模型方法 #

方法1：单个因变量，按Region批量跑

以因变量var1为例，用tidyverse批量拟合并整理结构化结果：

library(broom)

# 按region分组拟合GLM，自动整理为数据框
single_var_results <- df_filtered %>%
  group_by(region) %>%
  group_map(~ glm(var1 ~ grade, data = .x), .keep = TRUE) %>%
  map_dfr(tidy, .id = "region")

# 查看结果，包含每个区域的系数、标准误、p值等关键指标
print(single_var_results)

方法2：多个因变量，批量遍历所有组合

如果有多个因变量（比如var1到var10），可以一次性处理所有因变量+区域的组合：

# 定义所有因变量列表
dv_list <- paste0("var", 1:10)

# 遍历每个因变量，按region拟合模型并合并结果
multi_var_results <- map_dfr(dv_list, function(dv) {
  df_filtered %>%
    group_by(region) %>%
    group_map(~ glm(as.formula(paste(dv, "~ grade")), data = .x), .keep = TRUE) %>%
    map_dfr(tidy, .id = "region") %>%
    mutate(dependent_var = dv)
})

# 查看所有结果，包含因变量名、区域、模型参数等信息
print(multi_var_results)

方法3：Base R实现（无需tidyverse）

如果不习惯用tidyverse生态，用Base R也能完成批量处理：

# 拆分数据为各个region的子集
region_data <- split(df_filtered, df_filtered$region)

# 批量拟合模型
models <- lapply(region_data, function(x) glm(var1 ~ grade, data = x))

# 整理结果为结构化数据框
base_results <- do.call(rbind, lapply(names(models), function(reg) {
  broom::tidy(models[[reg]]) %>% mutate(region = reg)
}))

扩展到State层级 #

如果要按state层级拟合模型，只需要把代码中的group_by(region)替换为group_by(state)即可；如果需要同时按region+state二级地理层级分组，改为group_by(region, state)即可。

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