如何编写函数将GMM检验结果自动添加至stargazer输出？

需求背景 #

用plm包做大量GMM回归分析时，需要报告Sargan检验、AR(1)/AR(2)检验、Wald检验等结果，但stargazer无法直接输出这些统计量。手动编写代码逐个提取并添加效率极低，因此需要一个可复用的函数，传入回归方程对象后自动生成所有检验结果，直接用于stargazer的add.lines参数。

解决方案代码 #

1. 定义自动提取检验结果的函数 #

library(plm)
library(stargazer)
library(stringr)

# 定义函数：输入pgmm对象列表和对应的工具变量描述，输出stargazer可用的add.lines列表
extract_gmm_tests <- function(pgmm_list, inst_descriptions) {
  # 初始化结果列表
  add_lines_list <- list()
  
  # 处理工具变量行
  add_lines_list[["Instruments"]] <- c("Instr", inst_descriptions)
  
  # 遍历每个pgmm对象，提取检验结果
  summaries <- lapply(pgmm_list, function(x) summary(x, robust = TRUE))
  
  # 提取Sargan检验结果
  sargan_stats <- sapply(summaries, function(s) round(s$sargan$statistic, 4))
  sargan_pvals <- sapply(summaries, function(s) round(s$sargan$p.value, 4))
  sargan_conclusion <- ifelse(sargan_pvals > 0.05, "Good", "Bad")
  
  add_lines_list[["Sargan_stat"]] <- c("Sargan", sargan_stats)
  add_lines_list[["Sargan_p"]] <- c("Sargan (p-value)", sargan_pvals)
  add_lines_list[["Sargan_conclusion"]] <- c("Conclusion about Sargan", sargan_conclusion)
  
  # 提取AR(1)检验结果
  ar1_stats <- sapply(summaries, function(s) round(s$m1$statistic, 2))
  ar1_pvals <- sapply(summaries, function(s) round(s$m1$p.value, 4))
  
  add_lines_list[["AR1_stat"]] <- c("AR(1)", ar1_stats)
  add_lines_list[["AR1_p"]] <- c("AR(1) p-value", ar1_pvals)
  
  # 提取AR(2)检验结果
  ar2_stats <- sapply(summaries, function(s) round(s$m2$statistic, 2))
  ar2_pvals <- sapply(summaries, function(s) round(s$m2$p.value, 4))
  
  add_lines_list[["AR2_stat"]] <- c("AR(2)", ar2_stats)
  add_lines_list[["AR2_p"]] <- c("AR(2) p-value", ar2_pvals)
  
  # 提取Wald系数检验结果
  wald_coef_stats <- sapply(summaries, function(s) {
    str_c(round(s$wald.coef$statistic, 2), " (", s$wald.coef$parameter, ")")
  })
  wald_coef_pvals <- sapply(summaries, function(s) round(s$wald.coef$p.value, 4))
  
  add_lines_list[["Wald_coef"]] <- c("Wald Coef (df)", wald_coef_stats)
  add_lines_list[["Wald_coef_p"]] <- c("Wald Coef p-value", wald_coef_pvals)
  
  # 提取Wald时间效应检验结果
  wald_td_stats <- sapply(summaries, function(s) {
    str_c(round(s$wald.td$statistic, 2), " (", s$wald.td$parameter, ")")
  })
  wald_td_pvals <- sapply(summaries, function(s) round(s$wald.td$p.value, 4))
  
  add_lines_list[["Wald_time"]] <- c("Wald Time (df)", wald_td_stats)
  add_lines_list[["Wald_time_p"]] <- c("Wald Time p-value", wald_td_pvals)
  
  # 返回整理好的列表
  return(add_lines_list)
}

2. 使用示例 #

# 加载数据并拟合GMM模型
data("EmplUK")

# Arellano Bond 91 table 4b 
z1 <- pgmm(dynformula(log(emp)~log(wage)+log(capital)+log(output),list(2,1,0,1)),
           data=EmplUK, effect="twoways", model="twosteps",
           gmm.inst=~log(emp),lag.gmm=list(c(2,99)))

# Blundell and Bond tab4 (cf DPD for OX p.12 col.4)
z2 <- pgmm(dynformula(log(emp)~log(wage)+log(capital),list(1,1,1)),
           data=EmplUK, effect="twoways", model="onestep",
           gmm.inst=~log(emp)+log(wage)+log(capital),lag.gmm=c(2,99),
           transformation="ld")

# 准备工具变量描述
inst_descs <- c(
  "~log(emp),lag.gmm=list(c(2,99))",
  "~log(emp)+log(wage)+log(capital),lag.gmm=c(2,99)"
)

# 调用函数生成检验结果
gmm_test_lines <- extract_gmm_tests(list(z1, z2), inst_descs)

# 用stargazer输出结果
stargazer(z1, z2,
          column.labels = c("Eq z1", "Eq z2"),
          title= "2 GMM equations",
          align = TRUE,
          type = "text",
          digits = 4,
          no.space = TRUE,
          add.lines = gmm_test_lines,
          notes="Sargan检验p值大于5%时视为可接受（Good）")

函数说明 #

• 输入参数：pgmm_list是拟合好的pgmm对象组成的列表；inst_descriptions是对应每个模型的工具变量设置描述字符串组成的向量。
• 输出：直接可传入stargazer的add.lines参数的列表，包含所有需要的检验统计量、p值和结论。

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