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.

Goeselt

Are you opposed to doing it by hand using the add.lines option? Then you could use the coeftest object and add the other stats from the lm object: stargazer(m,type="text", single.row = T,add.lines = list(c("Observations",length(m1$data[,1])),c("Null Deviance" ,round(m1$null.deviance,3)))) – Cadency

Thanks @Cadency but in my real example that's too cumbersome. I'm actually inputting lists with N models inside the list. It would be too much manual work for 5-6 models. – Goeselt

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)}

Ame

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!

Marquita

I don't think that this will work...did you try it? Stargazer produces text (or html or latex), rather than objects that can be manipulated. The assignment to stargazer.values just assigns the text the stargazer command produces. – Cadency

@BrorGiesenbauer, @Cadency is right. If you type, for example, str(m), you can't see the insides of the stargazer object but rather text. – Goeselt

Ah, good to know, I had misread your question and supposed that stargazer returns values of class "stargazer". To reiterate: Does coeftest show all relevant values and stargazer just omits some of them that you'd like to be included in your report? – Marquita

One thing that might help is to include some general values (like n) in the notesection: stargazer(m, type = "text", notes = paste("n =", nrow(mtcars), sep =" "), notes.align = "l", single.row = T) – Marquita

Well, if I append the values to the coeftest object stargazer doesn't recognize it. I really have no clue how stargazer picks estimations. As for the previous comment, in my particular example this is cumbersome since I have lists with several models inside. – Goeselt

Three years later, I also could not figure out how to get stargazer and coeftest to provide . But it looks like estimatr can provide non-vanilla homogonous standard errors similar to coeftest and has a built-in helper function starprep to be compatible with stargazer. I don't think this would work with a glm obejct though. Here's the resource: declaredesign.org/r/estimatr/articles/regression-tables.html – Plutocrat