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?

Maidservant

I'm pretty sure none of those statistics depend on the variance-covariance matrix, just the residuals and the variance of y, sample size, df, etc. So the output would be the same. see stats.stackexchange.com/questions/5135/… – Breastwork

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?

Witty

Thanks @marco! But this only explains how to manually add the F-statistics. But furst i have to know my f statistics of the overall regression with the robust SE. Does it remain the same as without RSE? – Maidservant

Hi @Maidservant in the sense of stackoverflow, please accept my answer, when it solves the initial problem (as it says in the heading, robust standard errors in stargazer). I think you can try the wald test, instead of F statistic. If you need more information about some econometric background, I suggest to open a new thread on cross-validated. Best regards – Witty

Hi, I'm looking for a way to include robust SEs in stargazer from link model output. It has robust and cluster option, but I can't get it to print the robust SEs and pvalues in stargazer. – Snowshoe

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.

Hejaz

I found this to be a better solution as it provides p-values and t-distribution from the robust standard errors. The answer above only changes standard errors whilst still displaying p-values based on the standard t-distribution. Be careful! – Catena