how to create one output with multiple linear regression models including Cluster-Robust Standard Errors in Stargazer

Question

I know how I can create a linear regression model lm and how to use the summary function to obtain clustered standard errors and add them to the stargazer output:

# 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")

Does anyone know how the code supposed to look like if I want to create one stargazer output with multiple regression models and clustered standard errors?

The logic of the code is as follows:

1 step: create lm models

ols1 <- lm(y ~ x) 
ols2 <- lm(y ~ x + z) 
ols3 <- lm(y ~ x + z + a) 
ols2 <- lm(y ~ x + z + a + b)

2 step: include standard errors

summary(ols1, cluster="cluster_id")
summary(ols2, cluster="cluster_id")
summary(ols3, cluster="cluster_id")
summary(ols4, cluster="cluster_id")

3 step: create one output with 4 different models

stargazer(ols1,ols2,ols3,ols4, type="html", dep.var.labels=c("ROA"), intercept.bottom = FALSE,
          out="OLS1")

I think step 1 and step 2 are not critical but I do not know how to set up the code for step 3.

I do not know how to implement the following code in step 3:

# create table in stargazer
stargazer(ols1, se=list(coef(summary(ols1,cluster = c("cluster_id")))[, 2]), type = "text")

Can anyone help?

Thank you so much!!!

Marco Marco · Accepted Answer · 2020-01-14T16:09:08

You already saw that the option se can handle a list. So just save the robust SE somewhere else and easily accessible, such as:

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

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

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

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

# Adjust standard errors
cov1         <- vcovHC(model1, type = "HC1")
cov2         <- vcovHC(model2, type = "HC1")
robust_se1    <- sqrt(diag(cov1))
robust_se2    <- sqrt(diag(cov2))

# Stargazer output (with and without RSE)
stargazer(model1, model2, type = "text",
          se = list(robust_se1, robust_se2))

Now you can have as any models next to each other with robust SE in stargazer.

how to create one output with multiple linear regression models including Cluster-Robust Standard Errors in Stargazer

1 Answers