make R report adjusted R squared and F-test in output with robust standard errors

Question

I have estimated a linear regression model using lm(x~y1 + y1 + ... + yn) and to counter the present heteroscedasticity I had R estimate the robust standard errors with

coeftest(model, vcov = vcovHC(model, type = "HC0"))

I know that (robust) R squared and F statistic from the "normal" model are still valid, but how do I get R to report them in the output? I want to fuse several regression output from different specifications together with stargazer and it would become very chaotic if I had to enter the non-robust model along just to get these statistics. Ideally I want to enter a regression output into stargazer that contains these statistics, thus importing it to their framework.

Thanks in advance for all answers

Answer 1

This is how to go about it. You need to use model object that is supported by stargazer as a template and then you can provide a list with standard errors to be used:

library(dplyr)
library(lmtest)
library(stargazer)


# Basic Model ---------------------------------------------------------------------------------

model1 <- lm(hp ~ factor(gear) + qsec + cyl + factor(am), data = mtcars)
summary(model1)


# Robust standard Errors ----------------------------------------------------------------------

model_robust <- coeftest(model1, vcov = vcovHC(model1, type = "HC0"))

# Get robust standard Errors (sqrt of diagonal element of variance-covariance matrix)
se = vcovHC(model1, type = "HC0") %>% diag() %>% sqrt()

stargazer(model1, model1, 
          se = list(NULL, se), type = 'text')

Using this approach you can use stargazer even for model objects that are not supported. You only need coefficients, standard errors and p-values as vectors. Then you can 'mechanically insert' even unsupported models.

One last Note. You are correct that once heteroskedasticity is present, Rsquared can still be used. However, overall F-test as well as t-tests are NOT valid anymore.

Answer 2

I don't have a solution with stargarzer, but I do have a couple of viable alternatives for regression tables with robust standard errors:

Option 1

Use the modelsummary package to make your tables.

it has a statistic_override argument which allows you to supply a function that calculates a robust variance covariance matrix (e.g., sandwich::vcovHC.

library(modelsummary)
library(sandwich)

mod1 <- lm(drat ~ mpg, mtcars)
mod2 <- lm(drat ~ mpg + vs, mtcars)
mod3 <- lm(drat ~ mpg + vs + hp, mtcars)
models <- list(mod1, mod2, mod3)

modelsummary(models, statistic_override = vcovHC)

Note 1: The screenshot above is from an HTML table, but the modelsummary package can also save Word, LaTeX or markdown tables.

Note 2: I am the author of this package, so please treat this as a potentially biased view.

Option 2

Use the estimatr::lm_robust function, which automatically includes robust standard errors. I believe that estimatr is supported by stargazer, but I know that it is supported by modelsummary.

library(estimatr)

mod1 <- lm_robust(drat ~ mpg, mtcars)
mod2 <- lm_robust(drat ~ mpg + vs, mtcars)
mod3 <- lm_robust(drat ~ mpg + vs + hp, mtcars)
models <- list(mod1, mod2, mod3)

modelsummary(models)