13

Does anyone know how to get stargazer to display clustered SEs for lm models? (And the corresponding F-test?) If possible, I'd like to follow an approach similar to computing heteroskedasticity-robust SEs with sandwich and popping them into stargazer as in http://jakeruss.com/cheatsheets/stargazer.html#robust-standard-errors-replicating-statas-robust-option.

I'm using lm to get my regression models, and I'm clustering by firm (a factor variable that I'm not including in the regression models). I also have a bunch of NA values, which makes me think multiwayvcov is going to be the best package (see the bottom of landroni's answer here - Double clustered standard errors for panel data - and also https://sites.google.com/site/npgraham1/research/code)? Note that I do not want to use plm.

Edit: I think I found a solution using the multiwayvcov package...

library(lmtest) # load packages
library(multiwayvcov)

data(petersen) # load data
petersen$z <- petersen$y + 0.35  # create new variable

ols1 <- lm(y ~ x, data = petersen) # create models
ols2 <- lm(y ~ x + z, data = petersen)

cl.cov1 <- cluster.vcov(ols1, data$firmid) # cluster-robust SEs for ols1
cl.robust.se.1 <- sqrt(diag(cl.cov1))
cl.wald1 <- waldtest(ols1, vcov = cl.cov1)

cl.cov2 <- cluster.vcov(ols2, data$ticker) # cluster-robust SEs for ols2
cl.robust.se.2 <- sqrt(diag(cl.cov2))
cl.wald2 <- waldtest(ols2, vcov = cl.cov2)

stargazer(ols1, ols2, se=list(cl.robust.se.1, cl.robust.se.2), type = "text") # create table in stargazer

Only downside of this approach is you have to manually re-enter the F-stats from the waldtest() output for each model.

Community
  • 1
  • 1
RSS
  • 163
  • 1
  • 2
  • 11

2 Answers2

17

Using the packages lmtest and multiwayvcov causes a lot of unnecessary overhead. The easiest way to compute clustered standard errors in R is the modified summary() function. This function allows you to add an additional parameter, called cluster, to the conventional summary() function. The following post describes how to use this function to compute clustered standard errors in R:

https://economictheoryblog.com/2016/12/13/clustered-standard-errors-in-r/

You can easily the summary function to obtain clustered standard errors and add them to the stargazer output. Based on your example you could simply use the following code:

# estimate models
ols1 <- lm(y ~ x) 

# summary with cluster-robust SEs
summary(ols1, cluster="cluster_id") 

# create table in stargazer
stargazer(ols1, se=list(coef(summary(ols1,cluster = c("cluster_id")))[, 2]), type = "text")
Flavian Barton
  • 186
  • 2
  • 4
  • 1
    I want to create one output with 4 regression models. Each model should have clustered standard errors. How would you set up the code? This is how it looks like for the model without clustered standard errors: stargazer(m1,m2,m3,m4, type="html", dep.var.labels=c("ROA"), intercept.bottom = FALSE, out="OLS1") – Yaron Nolan Dec 27 '19 at 17:24
  • 1
    @YaronNolan I'd do the following: clust <- function(m, cluster="cluster_id") { return(coef(summary(m, cluster = c(cluster)))[, 2]) } stargazer(m1, m2, m3, m4, se=list(clust(m1), clust(m2), clust(m3), clust(m4)), ...) – Max Ghenis May 03 '20 at 15:42
  • I used this function, but it gave me the same standard errors as the one without cluster option. Also felm gave another error (very small), and stata gave ne a slightly larger error. No idea why that happens, but might be a good idea to first compare several tools for consistency. – Maxim Moloshenko Apr 13 '21 at 11:25
  • Note that the code there doesn't work anymore. – Ken Lee Jun 01 '21 at 20:45
5

I would recommend lfe package, which is much more powerful package than lm package. You can easily specify the cluster in the regression model:

ols1 <- felm(y ~ x + z|0|0|firmid, data = petersen)
summary(ols1)

stargazer(OLS1, type="html")

The clustered standard errors will be automatically produced. And stargazer will report the clustered-standard error accordingly.

By the way (allow me to do more marketing), for micro-econometric analysis, felm is highly recommended. You can specify fixed effects and IV easily using felm. The grammar is like:

ols1 <- felm(y ~ x + z|FixedEffect1 + FixedEffect2 | IV | Cluster, data = Data)
Xing Zhang
  • 175
  • 1
  • 5
  • when I am using the "lm" command my industry and year dummies are displayed. When I m using the "felm" command my industry and yer dummies are not beeing displayed. Does r still take these dummies into account when using your code? This is how my code look like: – Yaron Nolan Dec 27 '19 at 19:47
  • m1_1 <- felm(ROA ~ fam_ownership + lag_investment + dual_class + age + crisis |0 | 0 | gvkey, data) m1_2 <- felm(ROA ~ fam_ownership + fam_ownership_squared + lag_investment + dual_class + age+crisis |0 | 0 | gvkey, data) m1_3 <- felm(ROA ~ fam_ownership + fam_ownership_squared + lag_investment + dual_class + age +crisis |as.factor(industry)+as.factor(year)|0| gvkey, data) stargazer(m1_1,m1_2,m1_3, type="html", dep.var.labels=c("ROA"), intercept.bottom = FALSE, out="OLS1") – Yaron Nolan Dec 27 '19 at 19:48
  • I just love the fact that the `felm` model output works with `stargazer`, which saved me much time. I also checked that the estimated robust std. errors are the same as STATA output. – MECoskun Mar 30 '23 at 03:20