将统计信息附加到 coeftest 输出以包含在 stargazer 表中

Question

我有一个 glm 模型，我使用来自 lmtest 包的 coeftest 来估计稳健的标准误差。当我使用 stargazer 生成回归表时，我得到了正确的结果，但没有观察次数和其他相关统计数据，如零偏差和模型偏差。

这是一个例子：

library(lmtest)
library(stargazer)

m1 <- glm(am ~ mpg + cyl + disp, mtcars, family = binomial)
# Simple binomial regression

# For whatever reason, let's say I want to use coeftest to estimate something
m <- coeftest(m1)

stargazer(m, type = "text", single.row = T) # This is fine, but I want to also include the number of observations
                                            # the null deviance and the model deviance.

我对观察次数、零偏差和残余偏差特别感兴趣。

我认为如果我用新的系数矩阵替换旧的系数矩阵，我会得到具有正确统计数据的正确估计值，并且 stargazer 会识别模型并正确打印它。为此，我尝试在 m1 模型中替换 coeftest 模型中的系数、SE、z 统计量和 p 值，但其中一些统计量是用 summary.glm 计算的，不包括在 m1 中输出。我可以在summary 输出中轻松替换这些系数，但 stargazer 无法识别摘要类型类。我尝试使用特定统计信息向 m 对象添加属性，但它们没有显示在输出中，并且 stargazer 无法识别它。

注意：我知道 stargazer 可以计算稳健的 SE，但我也在进行其他计算，因此示例需要包含 coeftest 输出。

感谢任何帮助。

Answer 1

将原始模型传递给 stargazer 可能是最简单的方法，然后使用 coeftest 传递标准误差 (se = )、置信区间 (ci.custom = ) 和/或 p 值 (@987654325) 的自定义值@)。请参阅下文，了解如何轻松处理包含多个模型的列表。

suppressPackageStartupMessages(library(lmtest))
suppressPackageStartupMessages(library(stargazer))

mdls <- list(
  m1 = glm(am ~ mpg, mtcars, family = poisson),
  
  m2 = glm(am ~ mpg + cyl + disp, mtcars, family = poisson)
)

# Calculate robust confidence intervals
se_robust <- function(x)
  coeftest(x, vcov. = sandwich::sandwich)[, "Std. Error"]


# Original SE
stargazer(mdls, type = "text", single.row = T, report = "vcsp")
#> 
#> ===============================================
#>                        Dependent variable:     
#>                   -----------------------------
#>                                am              
#>                        (1)            (2)      
#> -----------------------------------------------
#> mpg               0.106 (0.042)  0.028 (0.083) 
#>                     p = 0.012      p = 0.742   
#> cyl                              0.435 (0.496) 
#>                                    p = 0.381   
#> disp                             -0.014 (0.009)
#>                                    p = 0.151   
#> Constant          -3.247 (1.064) -1.488 (3.411)
#>                     p = 0.003      p = 0.663   
#> -----------------------------------------------
#> Observations            32             32      
#> Log Likelihood       -21.647        -20.299    
#> Akaike Inf. Crit.     47.293         48.598    
#> ===============================================
#> Note:               *p<0.1; **p<0.05; ***p<0.01

# With robust SE
stargazer(
  mdls, type = "text", single.row = TRUE, report = "vcsp", 
  se = lapply(mdls, se_robust))
#> 
#> ===============================================
#>                        Dependent variable:     
#>                   -----------------------------
#>                                am              
#>                        (1)            (2)      
#> -----------------------------------------------
#> mpg               0.106 (0.025)  0.028 (0.047) 
#>                    p = 0.00002     p = 0.560   
#> cyl                              0.435 (0.292) 
#>                                    p = 0.137   
#> disp                             -0.014 (0.007)
#>                                    p = 0.042   
#> Constant          -3.247 (0.737) -1.488 (2.162)
#>                    p = 0.00002     p = 0.492   
#> -----------------------------------------------
#> Observations            32             32      
#> Log Likelihood       -21.647        -20.299    
#> Akaike Inf. Crit.     47.293         48.598    
#> ===============================================
#> Note:               *p<0.1; **p<0.05; ***p<0.01

^{由reprex package (v0.3.0) 于 2020 年 11 月 9 日创建}

Answer 2

如果我猜对了，您可以尝试以下方法：

首先，将您的观星分析分配给这样的对象

stargazer.values <- stargazer(m, type = "text", single.row = T) 

然后用body(stargazer)检查stargazer命令的代码。 希望您能找到观星者使用但未报告的值的对象。然后，您可以像这样处理它们（例如，如果有一个名为“null.deviance”的对象

stargazers.values$null.deviance

或者，如果它是另一个数据框的一部分，比如 df，它可以是这样的

stargazers.values$df$null.deviance

也许这样的代码会有所帮助

print(null.deviance <- stargazers.values$null.deviance)

希望这会有所帮助！