R：稳健的 SE 和 stargazer 中的模型诊断 table

Question

我尝试使用 stargazer 包将 AER 包中通过 ivreg() 生成的一些 2SLS 回归输出放入 Latex 文档中。我有几个问题，但我自己似乎无法解决。

我不知道如何插入 ivreg() 摘要中提供的模型诊断。即弱仪器测试、Wu-Hausmann 和 Sargan 测试。我想让它们包含通常在 table 下方报告的统计数据，例如观察数、R 平方和 Resid。 SE。 stargazer 函数似乎没有参数，您可以在其中提供包含额外诊断的列表。我没有把它放到我的示例中，因为老实说，我不知道从哪里开始。
我想用稳健标准误差交换正常标准误差，我发现唯一的方法是生成具有稳健标准误差的对象，并将它们添加到 stargazer() 函数中 se=list().我将其放入下面的最小工作示例中。是否有更优雅的方式对此进行编码，或者重新估计模型并使用稳健的标准误差保存它？

library(AER)
library(stargazer)

y <- rnorm(100, 5, 10)
x <- rnorm(100, 3, 15)
z <- rnorm(100, 3, 7)
a <- rnorm(100, 1, 7)
b <- rnorm(100, 3, 5)

# Fitting IV models
fit1 <- ivreg(y ~ x + a  |
             a + z,
             model = TRUE)
fit2 <- ivreg(y ~ x + a  |
             a + b + z,
             model = TRUE)

# Here are the se's and the diagnostics i want
summary(fit1, vcov = sandwich, diagnostics=T)
summary(fit2, vcov = sandwich, diagnostics=T)

# Getting robust se's, i think HC0 is the standard
# used with "vcov=sandwich" from the  above summary
cov1        <- vcovHC(fit1, type = "HC0")
robust1     <- sqrt(diag(cov1))
cov2        <- vcovHC(fit2, type = "HC0")
robust2     <- sqrt(diag(cov1))

# Create latex table
stargazer(fit1, fit2, type = "latex", se=list(robust1, robust2))

Answer 1

这里有一种方法可以做你想做的事：

require(lmtest)

rob.fit1        <- coeftest(fit1, function(x) vcovHC(x, type="HC0"))
rob.fit2        <- coeftest(fit2, function(x) vcovHC(x, type="HC0"))
summ.fit1 <- summary(fit1, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
summ.fit2 <- summary(fit2, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)

stargazer(fit1, fit2, type = "text", 
          se = list(rob.fit1[,"Std. Error"], rob.fit2[,"Std. Error"]), 
          add.lines = list(c(rownames(summ.fit1$diagnostics)[1], 
                             round(summ.fit1$diagnostics[1, "p-value"], 2), 
                             round(summ.fit2$diagnostics[1, "p-value"], 2)), 
                           c(rownames(summ.fit1$diagnostics)[2], 
                             round(summ.fit1$diagnostics[2, "p-value"], 2), 
                             round(summ.fit2$diagnostics[2, "p-value"], 2)) ))

这将产生：

==========================================================
                                  Dependent variable:     
                              ----------------------------
                                           y              
                                   (1)            (2)     
----------------------------------------------------------
x                                 -1.222        -0.912    
                                 (1.672)        (1.002)   

a                                 -0.240        -0.208    
                                 (0.301)        (0.243)   

Constant                          9.662         8.450**   
                                 (6.912)        (4.222)   

----------------------------------------------------------
Weak instruments                   0.45          0.56     
Wu-Hausman                         0.11          0.18     
Observations                       100            100     
R2                                -4.414        -2.458    
Adjusted R2                       -4.526        -2.529    
Residual Std. Error (df = 97)     22.075        17.641    
==========================================================
Note:                          *p<0.1; **p<0.05; ***p<0.01

如您所见，这允许手动将诊断包含在各个模型中。

您可以通过创建一个函数来自动执行此方法，该函数接收模型列表（例如 list(summ.fit1, summ.fit2)）并输出 se 或 add.lines 参数所需的对象。

gaze.coeft <- function(x, col="Std. Error"){
    stopifnot(is.list(x))
    out <- lapply(x, function(y){
        y[ , col]
    })
    return(out)
}
gaze.coeft(list(rob.fit1, rob.fit2))
gaze.coeft(list(rob.fit1, rob.fit2), col=2)

将同时接收 list 个 coeftest 对象，并按照 se 的预期生成 SEs 向量：

[[1]]
(Intercept)           x           a 
  6.9124587   1.6716076   0.3011226 

[[2]]
(Intercept)           x           a 
  4.2221491   1.0016012   0.2434801

同样可以进行诊断：

gaze.lines.ivreg.diagn <- function(x, col="p-value", row=1:3, digits=2){
    stopifnot(is.list(x))
    out <- lapply(x, function(y){
        stopifnot(class(y)=="summary.ivreg")
        y$diagnostics[row, col, drop=FALSE]
    })
    out <- as.list(data.frame(t(as.data.frame(out)), check.names = FALSE))
    for(i in 1:length(out)){
        out[[i]] <- c(names(out)[i], round(out[[i]], digits=digits))
    }
    return(out)
}
gaze.lines.ivreg.diagn(list(summ.fit1, summ.fit2), row=1:2)
gaze.lines.ivreg.diagn(list(summ.fit1, summ.fit2), col=4, row=1:2, digits=2)

两次调用都会产生：

$`Weak instruments`
[1] "Weak instruments" "0.45"             "0.56"            

$`Wu-Hausman`
[1] "Wu-Hausman" "0.11"       "0.18"

现在 stargazer() 调用变得如此简单，产生与上面相同的输出：

stargazer(fit1, fit2, type = "text", 
      se = gaze.coeft(list(rob.fit1, rob.fit2)), 
      add.lines = gaze.lines.ivreg.diagn(list(summ.fit1, summ.fit2), row=1:2))

R：稳健的 SE 和 stargazer 中的模型诊断 table

R: Robust SE's and model diagnostics in stargazer table

latex

r

stargazer