How can I build nice tables of average partial effects of different generalized linear models with R?

I am estimating logit models with more than a few variables, and would like to neatly show average partial effects (APEs) for the model in this way:

Basically, show a table like the one that the stargazer command would produce for any kind of lm or glm object, but with APEs instead of slope coefficients and their standard errors rather than the ones for the slope coefficients.

My code goes something like this:

# Estimate the models 

fit1<-glm(ctol ~ y16 + polscore + age,
          data = df46,
          family = quasibinomial(link = 'logit'))

fit2<-glm(ctol ~ y16*polscore + age,
          data = df46,
          family = quasibinomial(link = 'probit'))

fit3<-glm(ctol ~ y16 + polscore + age + ed,
          data = df46,
          family = quasibinomial(link = 'logit'))

# Calculate marginal effects

me_fit1<-margins_summary(fit1)
me_fit2<-margins_summary(fit2)
me_fit3<-margins_summary(fit3)

The output of a margins_summary object, while itself a data.frame object, cannot just be passed to stargazer to produce the nice looking output it would do with a glm object, like fit1 in my code before.

> me_fit1

   factor     AME     SE       z      p   lower   upper
      age -0.0031 0.0005 -5.8426 0.0000 -0.0041 -0.0020
 polscore  0.0033 0.0031  1.0646 0.2871 -0.0028  0.0093
      y16  0.1184 0.0166  7.1271 0.0000  0.0859  0.1510

Trying to pass me_fit1 to stargazer simply prints the data.frame summary stats, as stargazer would normally do with objects of this type.

> stargazer(me_fit1, type = 'text')

=========================================================
Statistic N Mean  St. Dev.  Min   Pctl(25) Pctl(75)  Max 
---------------------------------------------------------
AME       3 0.040  0.068   -0.003  0.0001   0.061   0.118
SE        3 0.007  0.009   0.001   0.002    0.010   0.017
z         3 0.783  6.489   -5.843  -2.389   4.096   7.127
p         3 0.096  0.166     0       0       0.1      0  
lower     3 0.026  0.052   -0.004  -0.003   0.042   0.086
upper     3 0.053  0.085   -0.002  0.004    0.080   0.151
---------------------------------------------------------

I've tried using the coef and se options from stargazer to change the coefficients presented of stargazer(fit1) to APEs and their errors. While its simple to show APEs, trying to show their standard errors is problematic because it cannot find the names of the variables in order to match them with their coefficients (in this case, their APEs).

Please help! I haven't been able to present decent results because of this problem. You can see an MWE here.

You can do this by using a combination of the modelsummary and marginaleffects packages. (Massive Disclaimer: I maintain both packages.)

You can install modelsummary from CRAN:

install.packages("modelsummary")

You can install marginaleffects from Github (warning: this package is young and still in development):

library(remotes)
install_github("vincentarelbundock/marginaleffects")

Load libraries and fit two models. Store those two models in a list:

library(marginaleffects)
library(modelsummary)

mod <- list(
    glm(am ~ mpg, data = mtcars, family = binomial),
    glm(am ~ mpg + factor(cyl), data = mtcars, family = binomial))

Now, we want to apply the marginaleffects function to both models, so we use lapply to apply it to each element of the list:

mfx <- lapply(mod, marginaleffects)

Finally we call modelsummary with the output argument set to "markdown" because Markdown tables look good on Stack Overflow:

modelsummary(mfx, output = "markdown")