Converting broom outputs of multiple models into regression tables

For various reasons, I only have access to objects model_coeffs and model_summary below. These form broom::tidy() and broom::glance() objects respectively. In other words, I do not have access to the underlying lm/glm objects, but have access to model_coeffs and model_summary as .csv files. For a reproducible example, I have provided code so that others can get to the same point

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))

How can I convert model_coeffs and model_summary into something pretty like this? It does not have to be exactly that, plus it should append some model summary stuff. Can {modelsummary} do this? I tried to follow instructions here but to no avail.

EDIT: Thank you to @Vincent. Using the {modelsummary} package, you can trick it to believe that it's of class(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.')
)

The easiest way to do this might be to store your data in a format that modelsummary will automatically recognize: a modelsummary_list object. This strategy is described here:

https://modelsummary.com/vignettes/modelsummary_extension.html

Here’s a minimal example example:

library(modelsummary)

models <- list(
    glm(am ~ hp, data = mtcars, family = binomial),
    glm(am ~ hp + wt, data = mtcars, family = binomial)
)

results <- list()
for (i in seq_along(models)) {
    results[[i]] <- list()
    class(results[[i]]) <- "modelsummary_list"
    results[[i]][["tidy"]] <- data.frame(
        "term" = names(coef(models[[i]])),
        "estimate" = coef(models[[i]]),
        "OR" = exp(coef(models[[i]])),
        "std.error" = sqrt(diag(vcov(models[[i]]))))
    results[[i]][["glance"]] <- data.frame(nobs = nobs(models[[i]]))
}

modelsummary(results,
    statistic = c("OR", "std.error"),
    shape = term ~ model + statistic)

+-------------+--------+-------+-------+--------+--------------+-------+
|             | (1)                    | (2)                           |
+-------------+--------+-------+-------+--------+--------------+-------+
|             | Est.   | OR    | S.E.  | Est.   | OR           | S.E.  |
+=============+========+=======+=======+========+==============+=======+
| (Intercept) | 0.777  | 2.174 | 0.915 | 18.866 | 1.561455e+08 | 7.444 |
+-------------+--------+-------+-------+--------+--------------+-------+
| hp          | -0.008 | 0.992 | 0.006 | 0.036  | 1.037000e+00 | 0.018 |
+-------------+--------+-------+-------+--------+--------------+-------+
| wt          |        |       |       | -8.083 | 0.000000e+00 | 3.069 |
+-------------+--------+-------+-------+--------+--------------+-------+
| Num.Obs.    | 32     |       |       | 32     |              |       |
+-------------+--------+-------+-------+--------+--------------+-------+