如何用Stargazer整合probitmfx多组边际效应结果至同一表格？

1. 示例数据与模型拟合 #

直接用R内置的mtcars数据集，拟合三组不同自变量组合的probit边际效应模型：

# 加载必要包
library(mfx)
library(stargazer)

# 处理数据：将am（自动/手动变速箱）设为因变量
data(mtcars)
mtcars$am <- as.factor(mtcars$am)

# 拟合三个不同的probit边际效应模型
model1 <- probitmfx(am ~ mpg + wt, data = mtcars, atmean = TRUE)
model2 <- probitmfx(am ~ mpg + wt + hp, data = mtcars, atmean = TRUE)
model3 <- probitmfx(am ~ mpg + wt + hp + qsec, data = mtcars, atmean = TRUE)

2. 提取边际效应核心结果 #

probitmfx的输出结构特殊，Stargazer无法直接识别，所以手动提取每个模型的边际效应系数、标准误、p值，以及模型统计量：

# 写个小函数批量提取结果
get_mfx_results <- function(model) {
  coefs <- model$mfxest[, "Marginal Effect"]
  ses <- model$mfxest[, "Std. Error"]
  pvals <- model$mfxest[, "Pr(>|z|)"]
  n <- model$fit$nobs
  pseudo_r2 <- model$fit$pseudo.r2
  return(list(coefs = coefs, ses = ses, pvals = pvals, n = n, pseudo_r2 = pseudo_r2))
}

# 提取三个模型的结果
res1 <- get_mfx_results(model1)
res2 <- get_mfx_results(model2)
res3 <- get_mfx_results(model3)

3. 用Stargazer生成汇总表格 #

利用Stargazer的custom.*系列参数，把提取的结果手动传入，即可生成统一的汇总表格：

stargazer(
  type = "text", # 需要LaTeX格式时改为"latex"
  title = "三组Probit模型边际效应汇总",
  column.labels = c("模型1", "模型2", "模型3"),
  custom.coef = list(res1$coefs, res2$coefs, res3$coefs),
  custom.se = list(res1$ses, res2$ses, res3$ses),
  custom.p = list(res1$pvals, res2$pvals, res3$pvals),
  custom.n = c(res1$n, res2$n, res3$n),
  custom.loglik = c(model1$fit$loglik, model2$fit$loglik, model3$fit$loglik),
  custom.stat = list(c("伪R²" = res1$pseudo_r2), 
                     c("伪R²" = res2$pseudo_r2), 
                     c("伪R²" = res3$pseudo_r2)),
  star.cutoffs = c(0.05, 0.01, 0.001),
  digits = 3,
  omit.table.layout = "n"
)

4. 目标表格样式参考（文本输出示例） #

三组Probit模型边际效应汇总
=====================================================
                    模型1      模型2      模型3     
-----------------------------------------------------
mpg                0.051**    0.048**    0.060***   
                  (0.021)    (0.021)    (0.022)    
wt                -0.361***  -0.347***  -0.409***  
                  (0.116)    (0.118)    (0.129)    
hp                           -0.001     -0.002*    
                             (0.001)    (0.001)    
qsec                                      0.116**   
                                         (0.047)    
-----------------------------------------------------
观测数              32         32         32        
伪R²               0.397      0.402      0.495    
对数似然值        -12.832    -12.727    -11.144    
=====================================================
注: *p<0.05; **p<0.01; ***p<0.001

补充说明 #

• 若需嵌入LaTeX文档，将type = "text"改为type = "latex"即可，输出代码可直接使用。
• 如需调整变量顺序或添加自定义注释，可修改custom.coef的变量名，或给stargazer添加notes参数。
• 也可使用broom包的tidy()函数整理结果，但手动提取的方式更直接可控。

内容的提问来源于stack exchange，提问作者Es003