Appending statistics to coeftest output to include in stargazer tables

Question

I have a glm model for which I use coeftest from the lmtest package to estimate robust standard errors. When I use stargazer to produce regression tables I get the correct results but without the number of observations and other relevant statistics like the null deviance and the model deviance.

Here's an example:

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.

I'm specifically interested in the number of observations, the null deviance and the residual deviance.

I thought that If I replaced the old coefficient matrix with the new one, I'd get the correct estimates with the correct statistics and stargazer would recognize the model and print it correctly. For that, I've tried substituting the coefficients, SE's, z statistic and p values from the coeftest model in the m1 model but some of these statistics are computed with summary.glm and are not included in the m1 output. I could easily substitute these coefficients in the summary output but stargazer doesn't recognize summary type class. I've tried adding attributes to the m object with the specific statistics but they don't show up in the output and stargazer doesn't recognize it.

Note: I know stargazer can compute robust SE's but I'm also doing other computations, so the example needs to include the coeftest output.

Any help is appreciated.

Answer 1

It may be easiest to pass the original models into stargazer, and then use coeftest to pass in custom values for standard errors (se = ), confidence intervals (ci.custom = ) and/or p values (p = ). See below for how to easily handle a list containing multiple models.

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

^{Created on 2020-11-09 by the reprex package (v0.3.0)}

Answer 2

If I get you right, you could try the following:

First, assign your stargazer analysis to an object like this

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

then check the code of the stargazer command with body(stargazer). Hopefully you can find objects for values that stargazers uses but does not report. You can then address them like this (if there is, for example, an object named "null.deviance"

stargazers.values$null.deviance

Or, if it is part of another data frame, say df, it could go like this

stargazers.values$df$null.deviance

maybe a code like this could be helpful

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

Hope this helps!