Преобразование выходных данных нескольких моделей в таблицы регрессии

По разным причинам у меня есть доступ только к объектам model_coeffs и model_summary ниже. Они образуют объекты broom::tidy() и broom::glance() соответственно. Другими словами, у меня нет доступа к базовым объектам lm/glm, но есть доступ к model_coeffs и model_summary в виде файлов .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))

Как я могу преобразовать model_coeffs и model_summary во что-то вроде этого? Это не обязательно должно быть именно так, плюс к нему следует добавить некоторую сводную информацию о модели. Может ли {modelsummary} это сделать? Я пытался следовать инструкциям здесь , но безрезультатно.

Обновлено: Спасибо @Vincent. Используя пакет {modelsummary}, вы можете обмануть его, заставив поверить, что он принадлежит к классу (modelsummary_list)

# Filtering data for outcome 'vs'
model_coeffs_vs <- model_coeffs %>%
  filter(lhs == "vs")
model_summary_vs <- model_summary %>%
  filter(lhs == "vs")

# Creating a list of unique model numbers
model_numbers <- unique(model_coeffs_vs$model_number)

# Creating a list of modelsummary_list objects
model_list <- map(model_numbers, function(model_num) {

  # Extract tidy data
  tidy_data <- model_coeffs_vs %>%
    filter(model_number == !!model_num) %>%
    mutate("OR" = estimate,
           estimate = log(estimate))
  
  # Extract glance data
  glance_data <- model_summary_vs %>%
    filter(model_number == !!model_num) %>%
    select(-c(lhs, rhs, model_number, model_name)) %>%
    rename(
      `Sample size` = nobs,
      "Pseduo R<sup>2</sup> (Nagelkerke)" = r_squared_nagelkerke, 
      "Pseudo R<sup>2</sup> (Tjur)" = r_squared_tjur)
  
  # Create the modelsummary_list object
  mod <- list(
    tidy = tidy_data,
    glance = glance_data
  )
  
  class(mod) <- "modelsummary_list"
  
  return(mod)
})

# Naming models with model numbers
names(model_list) <- set_names(model_numbers)

# Use modelsummary to summarize the models

modelsummary(model_list,
             shape = term ~ model + statistic,
             coef_omit = "Intercept",
             statistic = c("Odds ratio" = "OR", 
                           "S.E." = "std.error"),
             coef_rename = c(
               "drat" = "Rear axle ratio", 
                             "hp" = "Horsepower",
                             "wt" = "Weight",
                             "qsec" = "1/4 mile time"),
             notes = list('Text of the first note.', 
                          'Text of the second note.')
)