如何使用 R 构建不同广义线性模型的平均部分效应的漂亮表格？

Question

我正在估计具有多个变量的 logit 模型，并希望以这种方式巧妙地显示模型的平均部分效应 (APE)：

基本上，显示一个类似于 stargazer 命令为任何类型的 lm 或 glm 对象生成的 table，但使用 APE 而不是斜率系数及其标准误差，而不是斜率系数的标准误差。

我的代码是这样的：

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

margins_summary 对象的输出虽然本身是一个 data.frame 对象，但不能仅传递给 stargazer 以产生与 [=15] 相同的漂亮输出=] 对象，就像我之前代码中的 fit1。

> 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) 的系数更改为 APE 及其错误。虽然显示 APE 很简单，但试图显示它们的标准误差是有问题的，因为它无法找到变量的名称以便将它们与它们的系数（在本例中是它们的 APE）相匹配。

请帮忙！由于这个问题，我一直无法提供像样的结果。你可以看到一个 MWE here.

Answer 1

您可以结合使用 modelsummary 和 marginaleffects 包来完成此操作。（重要免责声明：我维护这两个包。）

You can install modelsummary from CRAN:

install.packages("modelsummary")

You can install marginaleffects from Github（警告：这个包还很年轻，仍在开发中）：

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

加载库并拟合两个模型。将这两个模型存储在列表中：

library(marginaleffects)
library(modelsummary)

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

现在，我们要将 marginaleffects 函数应用于两个模型，因此我们使用 lapply 将其应用于列表的每个元素：

mfx <- lapply(mod, marginaleffects)

最后我们调用 modelsummary 并将 output 参数设置为 "markdown" 因为 Markdown 表格在 Stack Overflow 上看起来不错：

modelsummary(mfx, output = "markdown")

	Model 1	Model 2
mpg	0.046	0.056
	(0.017)	(0.035)
factor(cyl)6		0.097
		(0.174)
factor(cyl)8		0.093
		(0.237)
Num.Obs.	32	32
AIC	33.7	37.4
BIC	36.6	43.3
Log.Lik.	-14.838	-14.702

如何使用 R 构建不同广义线性模型的平均部分效应的漂亮表格？

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

r

standard-error

xtable

stargazer

coefficients