In a previous post, we discussed how to obtain clustered standard errors in R. While the previous post described how one can easily calculate cluster robust standard errors in R, this post shows how one can include cluster robust standard errors in stargazer and create nice tables including clustered standard errors.
I prepared a short tutorial to explain how to include clustered standard errors in stargazer. The following R code does the following. First, it loads the function that is necessary to compute clustered standard errors. Second, it downloads an example data set from this blog that is used for the OLS estimation and thirdly, it calculates a simple linear model using OLS. Finally, the script uses the summary.lm() function, the one that we loaded at the beginning, to calculate and recover STATA like clustered standard errors and passes them on to the stargazer function. In the example, I print the stargazer output as text, however, one replace can the argument type to “tex” or “html” in order to obtain perfectly formatted “tex” or “html” tables.
# start with an empty workspace
rm(list=ls())
# load necessary packages for importing the function
library(RCurl)
# load necessary packages for the example
library(gdata)
library(zoo)
# import the function
url_robust <- "https://raw.githubusercontent.com/IsidoreBeautrelet/economictheoryblog/master/robust_summary.R"
eval(parse(text = getURL(url_robust, ssl.verifypeer = FALSE)),
envir=.GlobalEnv)
# download data set for example
url_data <- "https://economictheoryblog.files.wordpress.com/2016/12/data.xls"
data <- read.xls(gsub("s:",":",url_data))
# estimate simple linear model
reg <- lm(id_score ~ class_size,
data=data)
# use new summary function to obtain clustered standard errors
summary(reg,cluster = c("class_id"))
# create stargazer output with cluster robust standard errors
require("stargazer")
# save cluster robust standard errors
cluster_se <- as.vector(summary(reg,cluster = c("class_id"))$coefficients[,"Std. Error"])
# print stargazer output with robust standard errors
stargazer(reg,type = "text",se = list(cluster_se))
# the last command prints the stargazer output (in this case as text)
# with cluster robust standard errors.
Economic Theory Blog 
Thank you so much for the code! Do you eventually also know how to report the clustered errors in something like mtable? I’m working with Markdown and knit everything Word and, hence, cannot use stargazer.
Hey Bea, unfortunately, I have not an already prepared way of using clustered standard errors and mtable. But you can, of course, use the clustered standard errors and pass them on to mtable.
Best, ad
Hey Ad, thank you so much for the quick reply! Best, Bea
I hope you manged to include the clustered standard errors in your table. Best, ad
Could you help me figure out how to modify the following line of code if there are more than one regressions? Thank you! cluster_se <- as.vector(summary(reg,cluster = c("class_id"))$coefficients[,"Std. Error"])
Hi, thanks for you comment. If there are one than one regression in you stargazer, you can simply add the standard errors to the list, i.e. list(cluster_se1,cluster_se2). Hope this works for you. Cheers, ad
Thanks for this, and thank you for this great public good!
Hi,
thanks a lot for your clustering function and this explaination.
I just have an additional question concerning the F-stat in stargazer. To my knowledge, it is not possible to update the clustered F-stat in stargazer. Do you have any idea how I can update the output?
Hi, I am not sure if there is an argument in the stargazer function through which you can directly pass the f test value. In case that there is not, why don’t you try to adjust the F-statistic in you regression object. Stargazer will then take the updated value. Best, ad
This is a great function, thanks for sharing it! However, I seemed to have broken it somehow. It worked fine until the ~20th time that I ran it again. Now the line you provide for calculating clustered SE does not work anymore:
cluster_se1 <- as.vector(summary(reg,cluster = c("classid"))$coefficients[,"Std. Error"])
Instead, now when I run the model, I get NaN values instead. I swear I didn't change anything at all in the model. Everything else runs fine.
Hey there. Thanks for you comment. Unfortunately, I need more information, i.e. a reproducible example, to help you out. I use the function and stargazer regularly, however, I never encountered similar problems. Best, ad
I had this problem a few times as well. Maybe we experienced the same kind of problem
This was because my dataset was not a dataset any more but of type “tbl_df” and in this case.