将统计信息添加到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也会识别模型并正确打印出来。为此，我尝试将coeftest模型中的系数、SE、z统计量和p值代入m1模型中，但其中一些统计量是用summary.glm计算的，不包括在m1输出中。我可以轻松地将这些系数代入summary输出中，但stargazer无法识别summary类型类。我尝试向m对象添加具有特定统计数据的属性，但它们不会显示在输出中，stargazer也无法识别它。注意：我知道stargazer可以计算稳健的SE，但我还要进行其他计算，因此示例需要包括coeftest输出。感谢任何帮助。

- cimentadaj

---

你反对手动使用 add.lines 选项吗？那么你可以使用 coeftest 对象，并从 lm 对象中添加其他统计信息：stargazer(m,type="text", single.row = T,add.lines = list(c("Observations",length(m1$data[,1])),c("Null Deviance" ,round(m1$null.deviance,3)))) - paqmo

谢谢@paqmo，但在我的实际例子中，这太麻烦了。我实际上正在输入包含N个模型的列表。对于5-6个模型来说，这将是太多的手动工作。 - cimentadaj

Answer 1

将原始模型传入stargazer可能是最简单的方法，然后使用coeftest传递自定义值来进行标准误差（se = ），置信区间（ci.custom = ）和/或p值（p = ）。请参见下面如何轻松处理包含多个模型的列表。

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 包 (v0.3.0) 在2020年11月9日创建。}

- JWilliman

Answer 2

如果我理解正确，您可以尝试以下操作：

首先，将您的stargazer分析分配给一个对象，如下所示：

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

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

stargazers.values$null.deviance

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

stargazers.values$df$null.deviance

也许像这样的代码会有帮助。

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

希望这有所帮助！

- Bror Giesenbauer

---

我不认为这会起作用...你试过吗？Stargazer生成文本（或HTML或LaTeX），而不是可以操作的对象。对stargazer.values的赋值只是将stargazer命令生成的文本分配给它。 - paqmo

@BrorGiesenbauer，@paqmo是正确的。例如，如果您键入str(m)，则无法看到stargazer对象内部而是文本。 - cimentadaj

啊，好知道了，我误读了你的问题，并假设stargazer返回"class stargazer"的值。再重申一遍：coeftest是否显示所有相关的值，而stargazer只是省略了一些你想在报告中包括的值？ - Bror Giesenbauer

一个可能有帮助的方法是在notesection中包含一些通用值（如n）： stargazer(m, type ="text", notes = paste ("n =", nrow(mtcars), sep =" "), notes.align = "l", single.row = T) - Bror Giesenbauer

好吧，如果我将值附加到“coeftest”对象中，“stargazer”无法识别它。我真的不知道stargazer如何选择估计。至于之前的评论，在我的特定示例中，这很麻烦，因为我有多个模型列表。 - cimentadaj

三年后，我仍然无法弄清楚如何获得stargazer和coeftest提供的非香草同质标准误差。但是看起来estimatr可以提供类似于coeftest的非香草同质标准误差，并且具有内置的帮助函数starprep以与stargazer兼容。不过我认为这可能无法用于glm对象。以下是资源链接：https://declaredesign.org/r/estimatr/articles/regression-tables.html - Ari Anisfeld