Robust Standard Errors in lm() using stargazer()

Question

I have read a lot about the pain of replicate the easy robust option from STATA to R to use robust standard errors. I replicated following approaches: StackExchange and Economic Theory Blog. They work but the problem I face is, if I want to print my results using the stargazer function (this prints the .tex code for Latex files).

Here is the illustration to my problem:

reg1 <-lm(rev~id + source + listed + country , data=data2_rev)
stargazer(reg1)

This prints the R output as .tex code (non-robust SE) If i want to use robust SE, i can do it with the sandwich package as follow:

vcov <- vcovHC(reg1, "HC1")

if I now use stargazer(vcov) only the output of the vcovHC function is printed and not the regression output itself.

With the package lmtest() it is possible to print at least the estimator, but not the observations, R2, adj. R2, Residual, Residual St.Error and the F-Statistics.

lmtest::coeftest(reg1, vcov. = sandwich::vcovHC(reg1, type = 'HC1'))

This gives the following output:

t test of coefficients:

            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -2.54923    6.85521 -0.3719 0.710611   
id           0.39634    0.12376  3.2026 0.001722 **
source       1.48164    4.20183  0.3526 0.724960   
country     -4.00398    4.00256 -1.0004 0.319041   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

How can I add or get an output with the following parameters as well?

Residual standard error: 17.43 on 127 degrees of freedom
Multiple R-squared:  0.09676,   Adjusted R-squared:  0.07543 
F-statistic: 4.535 on 3 and 127 DF,  p-value: 0.00469

Did anybody face the same problem and can help me out? How can I use robust standard errors in the lm function and apply the stargazer function?

Answer 1

You already calculated robust standard errors, and there's an easy way to include it in the stargazeroutput:

library("sandwich")
library("plm")
library("stargazer")

data("Produc", package = "plm")

# Regression    
model <- plm(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp,
             data = Produc, 
             index = c("state","year"),
             method="pooling")

# Adjust standard errors
cov1         <- vcovHC(model, type = "HC1")
robust_se    <- sqrt(diag(cov1))

# Stargazer output (with and without RSE)
stargazer(model, model, type = "text",
          se = list(NULL, robust_se))

Solution found here: https://www.jakeruss.com/cheatsheets/stargazer/#robust-standard-errors-replicating-statas-robust-option

Update I'm not so much into F-Tests. People are discussing those issues, e.g. https://stats.stackexchange.com/questions/93787/f-test-formula-under-robust-standard-error

When you follow http://www3.grips.ac.jp/~yamanota/Lecture_Note_9_Heteroskedasticity

"A heteroskedasticity-robust t statistic can be obtained by dividing an OSL estimator by its robust standard error (for zero null hypotheses). The usual F-statistic, however, is invalid. Instead, we need to use the heteroskedasticity-robust Wald statistic."

and use a Wald statistic here?

Answer 2

This is a fairly simple solution using coeftest:

reg1 <-lm(rev~id + source + listed + country , data=data2_rev)
cl_robust <- coeftest(reg1, vcov = vcovCL, type = "HC1", cluster = ~ 
   country)
se_robust <- cl_robust[, 2]
stargazer(reg1, reg1, cl_robust, se = list(NULL, se_robust, NULL))

Note that I only included cl_robust in the output as a verification that the results are identical.