如何使用R生成不同广义线性模型平均偏效应(APEs)的美观汇总表

我正在估计包含多个变量的logit模型，希望按如下形式清晰展示模型的平均偏效应（APEs）：

大体而言，我希望得到类似stargazer命令对lm或glm对象输出的汇总表格，但表格中要展示APEs而非斜率系数，对应的标准误也要替换为APEs的标准误而非原斜率系数的标准误。

我的代码示例如下：

# 估计模型
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'))

# 计算边际效应
me_fit1<-margins_summary(fit1)
me_fit2<-margins_summary(fit2)
me_fit3<-margins_summary(fit3)

margins_summary对象的输出本身是data.frame类型，但无法直接传入stargazer生成和之前代码中fit1这类glm对象一样的美观输出。

> 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

尝试将me_fit1传入stargazer只会输出data.frame的描述统计结果，这是stargazer处理此类对象的默认行为。

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

我试过使用stargazer的coef和se参数，将stargazer(fit1)输出的系数替换为APEs及对应标准误。展示APEs的操作很简单，但展示对应标准误时会出错，因为stargazer无法匹配变量名与对应的系数（此处为APEs）。

解决方案 #

方法1：适配stargazer的手动传参方案 #

通过构造带变量名的系数、标准误、p值列表，关闭stargazer默认的统计量自动计算逻辑即可实现需求，示例代码如下：

library(stargazer)
library(margins)

# 构造带变量名的系数、标准误、p值列表
coef_list <- list(
  setNames(me_fit1$AME, me_fit1$factor),
  setNames(me_fit2$AME, me_fit2$factor),
  setNames(me_fit3$AME, me_fit3$factor)
)
se_list <- list(
  setNames(me_fit1$SE, me_fit1$factor),
  setNames(me_fit2$SE, me_fit2$factor),
  setNames(me_fit3$SE, me_fit3$factor)
)
p_list <- list(
  setNames(me_fit1$p, me_fit1$factor),
  setNames(me_fit2$p, me_fit2$factor),
  setNames(me_fit3$p, me_fit3$factor)
)

# 调用stargazer输出表格
stargazer(fit1, fit2, fit3,
          type = "text", # 需要输出LaTeX/HTML格式可直接修改该参数
          coef = coef_list,
          se = se_list,
          p = p_list,
          # 自定义表格说明
          dep.var.caption = "平均边际效应(APE)",
          dep.var.labels = "ctol",
          column.labels = c("Logit主效应", "Probit交互项", "Logit加教育变量"),
          # 关闭默认统计量计算，使用传入的参数
          t.auto = FALSE,
          p.auto = FALSE,
          # 按需移除不需要的原始模型统计量
          omit.stat = c("aic", "ll", "resdev")
)

如果有交互项需要修改显示名，可补充covariate.labels参数自定义变量的展示名称。

方法2：使用原生支持边际效应的modelsummary包 #

如果可以更换输出包，modelsummary原生支持margins对象的输出，不需要手动做变量匹配，语法更简洁：

library(modelsummary)
# 直接传入边际效应对象即可生成规范表格
modelsummary(list(
  "Logit主效应" = me_fit1,
  "Probit交互项" = me_fit2,
  "Logit加教育变量" = me_fit3
), output = "markdown")

内容的提问来源于stack exchange，提问作者Daniel Sánchez