如何在嵌套数据上应用多公式并提取GLM模型结果

我正在尝试在嵌套数据上运行多个不同的公式，已完成大部分代码编写，但因不熟悉函数及嵌套数据处理，不知道如何提取列表中的模型信息。

现有代码如下：

#packages and data
library(tidyverse)
library(broom)
install.packages("AER")
library("AER")
data(Affairs, package="AER")

#creating response variable
Affairs$ynaffair[Affairs$affairs >  0] <- 1
Affairs$ynaffair[Affairs$affairs == 0] <- 0

#factoring some of the variables
Affairs <- Affairs %>% 
  mutate_at(c("affairs", "religiousness", "occupation", "rating", "ynaffair"), as.factor)

#nesting by religious group
by_relig_group <- Affairs %>%
  group_by(religiousness) %>%
  nest()

#Create formulas
gender_formula <- as.formula("ynaffair ~ gender")
age_formula <- as.formula("ynaffair ~ age")
occupation_formula <- as.formula("ynaffair ~ occupation")
formula_list <- list(gender_formula,age_formula, occupation_formula)

# Designate the models and function
model <- function(df, formula_list, index) {
    list(glm(formula = formula_list[[index]] , family=binomial, data = Affairs))
}

#Apply to each data group and formula list by formula number 
df1 <- by_relig_group %>% 
  mutate(model1 = map(data, ~model(., formula_list, 1)), 
         model2 = map(data, ~model(., formula_list, 2))) 

代码看似运行成功，但不知如何提取其中的信息。

在小范围测试中，我成功运行了以下代码：

#create nested data frame
by_relig_group <- Affairs %>%
  group_by(religiousness) %>%
  nest()

#create a function to run the linear models
model <- function(df){
  glm(ynaffair~gender + age + yearsmarried + children + religiousness + education + occupation + rating,
               family = binomial,
             data = Affairs)
}

#run the models
by_relig_group <- by_relig_group %>%
  mutate(model = map(data,model))


#unnest the models using the broom package 
by_relig_group_unnest <- by_relig_group %>%
  mutate(tidy = map(model, broom::tidy)) %>% ## I can't figure out how to apply these two lines of code above
  unnest(tidy)

我的最终目标是获取各水平频数、估计值（将其指数化得到OR）、置信区间（CI）及p值，并导出至Excel。同时希望得到在代码中获取置信区间的建议。

1. 修正模型运行代码 #

原代码存在两个核心问题：一是模型函数误用全局数据集Affairs而非传入的分组数据df，导致所有分组模型都用全量数据计算；二是将glm对象包裹在list中，增加后续处理复杂度。同时，分开创建model1/model2的方式效率较低，建议直接遍历公式列表批量处理。

修改后的代码：

# 加载包和数据
library(tidyverse)
library(broom)
library(AER)
data(Affairs, package="AER")

# 创建响应变量（简化写法）
Affairs$ynaffair <- ifelse(Affairs$affairs > 0, 1, 0) %>% as.factor()

# 转换因子变量
Affairs <- Affairs %>% 
  mutate_at(c("affairs", "religiousness", "occupation", "rating"), as.factor)

# 按宗教分组嵌套
by_relig_group <- Affairs %>%
  group_by(religiousness) %>%
  nest()

# 创建命名公式列表（方便后续识别模型类型）
formula_list <- list(
  gender_model = as.formula("ynaffair ~ gender"),
  age_model = as.formula("ynaffair ~ age"),
  occupation_model = as.formula("ynaffair ~ occupation")
)

# 定义模型函数：传入分组数据和单个公式，返回glm对象
run_single_model <- function(df, formula) {
  glm(formula = formula, family = binomial, data = df)
}

# 为每个分组运行所有公式的模型
by_relig_models <- by_relig_group %>%
  mutate(
    # 遍历公式列表，每个分组生成对应模型的列表
    models = map(data, ~map(formula_list, run_single_model, df = .x)),
    # 把模型列表拆分为单独列（可选，方便查看单个模型）
    models = unnest_wider(models)
  )

2. 提取模型信息（OR、CI、p值） #

使用broom::tidy()可直接提取模型系数，通过参数设置自动计算指数化的优势比（OR）和置信区间：

# 提取模型系数并计算OR、置信区间、p值
model_results <- by_relig_models %>%
  # 把多列模型转为长格式，统一处理
  pivot_longer(cols = starts_with("gender_model"):starts_with("occupation_model"),
               names_to = "model_type", values_to = "model") %>%
  # 用tidy提取结果：exponentiate=TRUE直接生成OR，conf.int=TRUE计算95%CI
  mutate(tidy_output = map(model, ~broom::tidy(.x, exponentiate = TRUE, conf.int = TRUE))) %>%
  unnest(tidy_output) %>%
  # 保留所需字段
  select(religiousness, model_type, term, 
         OR = estimate, CI_low = conf.low, CI_high = conf.high, p_value = p.value)

3. 获取各水平频数 #

在嵌套前计算每个宗教分组内的变量水平频数：

# 计算每个宗教分组内的变量水平频数
group_frequencies <- Affairs %>%
  group_by(religiousness) %>%
  # 按需选择需要统计的变量
  count(gender, occupation, name = "frequency") %>%
  ungroup()

# （可选）合并频数与模型结果
combined_results <- model_results %>%
  left_join(group_frequencies, by = c("religiousness", "term" = "gender")) %>%
  left_join(group_frequencies, by = c("religiousness", "term" = "occupation"), suffix = c("_gender", "_occupation"))

4. 导出至Excel #

使用writexl包可将多个数据框导出到Excel的不同工作表：

# 安装并加载包
install.packages("writexl")
library(writexl)

# 导出到Excel
write_xlsx(
  list(
    model_coefficients = model_results,
    group_counts = group_frequencies,
    combined_data = combined_results
  ),
  path = "affairs_analysis_results.xlsx"
)

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