如何基于broom输出的model_coeffs与model_summary生成美观回归表格？

问题背景 #

由于限制，仅能获取model_coeffs和model_summary两个对象（分别对应broom::tidy()和broom::glance()的输出），无法访问底层的lm/glm模型，且这两个对象仅能以CSV文件形式获取。

可复现代码（生成目标对象） #

library(tidyverse)

pairings <- tidyr::expand_grid(
  lhs = c('vs', 'am'),
  rhs = c('hp', 'hp + drat', 'hp + drat + wt', 'hp + drat + wt + qsec')
)

pairings <- pairings %>%
  dplyr::mutate(
    formula = map2(lhs, rhs, ~ as.formula(paste(.x, "~", .y))),
    model = map(formula, ~ glm(.x, family = binomial(link = "logit"), data = mtcars)),
    tidied_model = map(model, \(x) broom::tidy(x, conf.int = T, exp = T)),
    model_summaries = map(model, broom::glance),
    r_squared_nagelkerke = map(model, performance::r2_nagelkerke),
    r_squared_tjur = map(model, performance::r2_tjur),
    model_name = map_chr(formula, deparse),
    model_number = case_when(
      str_count(model_name, "\\+") == 0 ~ 'Model 1',
      str_count(model_name, "\\+") == 1 ~ 'Model 2',
      str_count(model_name, "\\+") == 2 ~ 'Model 3',
      str_count(model_name, "\\+") == 3 ~ 'Model 4'
    )
  ) %>% 
  dplyr::select(-formula, -model)

model_coeffs <- pairings %>% 
  dplyr::select(-model_summaries, -r_squared_nagelkerke, -r_squared_tjur) %>%
  tidyr::unnest(c(tidied_model),   names_repair = "unique") 

model_summary <- pairings %>% 
  dplyr::select(-tidied_model) %>% 
  tidyr::unnest(c(model_summaries, r_squared_nagelkerke, r_squared_tjur))

需求与解决方案 #

希望将上述两个对象转换为附带模型汇总信息的美观回归表格，尝试{modelsummary}包教程后未成功，以下是可行的解决方案代码：

解决方案代码 #

# 筛选因变量为'vs'的数据
model_coeffs_vs <- model_coeffs %>%
  filter(lhs == "vs")
model_summary_vs <- model_summary %>%
  filter(lhs == "vs")

# 获取唯一的模型编号列表
model_numbers <- unique(model_coeffs_vs$model_number)

# 创建modelsummary_list对象列表
model_list <- map(model_numbers, function(model_num) {

  # 提取tidy格式数据
  tidy_data <- model_coeffs_vs %>%
    filter(model_number == !!model_num) %>%
    mutate("OR" = estimate,
           estimate = log(estimate))
  
  # 提取glance格式数据并重命名
  glance_data <- model_summary_vs %>%
    filter(model_number == !!model_num) %>%
    select(-c(lhs, rhs, model_number, model_name)) %>%
    rename(
      `样本量` = nobs,
      "伪R²（Nagelkerke）" = r_squared_nagelkerke, 
      "伪R²（Tjur）" = r_squared_tjur)
  
  # 构建modelsummary_list对象
  mod <- list(
    tidy = tidy_data,
    glance = glance_data
  )
  
  class(mod) <- "modelsummary_list"
  
  return(mod)
})

# 为模型列表命名
names(model_list) <- set_names(model_numbers)

# 使用modelsummary生成汇总表格
modelsummary(model_list,
             shape = term ~ model + statistic,
             coef_omit = "Intercept",
             statistic = c("优势比" = "OR", 
                           "标准误" = "std.error"),
             coef_rename = c(
               "drat" = "后轴比", 
               "hp" = "马力",
               "wt" = "重量",
               "qsec" = "1/4英里耗时"),
             notes = list('第一条注释文本.', 
                          '第二条注释文本.')
)

内容的提问来源于stack exchange，提问作者allen.joseph