R doesn’t have a built-in method for calculating heteroskedastically-robust and clustered standard errors. Since this is a common need for us, I’ve put together a little library to provide these functions.<\/p>\n
Download the library<\/a>.<\/p>\n The file contains three functions:<\/p>\n The arguments for the functions are below, and then an example is at the end.<\/p>\n get.coef.robust(mod, var=c(), estimator=”HC”)<\/b> get.coef.clust(mod, cluster, var=c())<\/b> get.conf.clust(mod, cluster, xx, alpha, var=c())<\/b> The following example compare the results from Stata and R for an example from section 3.3.2 of “An Introduction to Categorical Data Analysis” by Alan Agresti. The data describes 173 female horseshoe crabs, and the number of male “satellites” that live near them.<\/p>\n Download the data as a CSV file<\/a>. The columns are the crabs color (2-5), spine condition (1-3), carapace width (cm), number of satellite crabs, and weight (kg).<\/p>\n The code below calculates a simple model with crab color fixed-effects, relating width (which is thought to be an explanatory variable) and number of satellites. The first model is a straight OLS regression; then we add robust errors, and finally crab color clustered errors.<\/p>\n The analysis in Stata<\/b>.<\/p>\n First, let’s do the analysis in Stata, to make sure that we can reproduce it in R.<\/p>\n import delim data.csv<\/p>\n xi: reg satell width i.color Here are excerpts of the results.<\/p>\n The OLS model:<\/p>\n The robust errors:<\/p>\n Note that we’ve used HC2 robust errors, to provide a direct comparison with the R results.<\/p>\n The clustered errors.<\/p>\n The robust and clustered errors are a little tighter for this example.<\/p>\n The analysis in R<\/b>.<\/p>\n Here’s the equivalent code in R.<\/p>\n And excerpts of the results. The OLS model:<\/p>\n The coefficient estimates will be the same for all three models. Stata gives us The robust errors:<\/p>\n Again, we use HC2 robust errors. Stata provides 2.939015 and .1028225 for the errors, matching these exactly.<\/p>\n The clustered errors:<\/p>\n Stata gives the exact same values as R.<\/p>\n","protected":false},"excerpt":{"rendered":" R doesn’t have a built-in method for calculating heteroskedastically-robust and clustered standard errors. Since this is a common need for us, I’ve put together a little library to provide these functions. Download the library. The file contains three functions: get.coef.robust: Calculate heteroskedastically-robust standard errors get.coef.clust: Calculate clustered standard errors get.conf.clust: Calculate confidence intervals for clustered … Continue reading Simple Robust Standard Errors Library for R<\/span> \n
get.coef.robust<\/code>: Calculate heteroskedastically-robust standard errors<\/li>\nget.coef.clust<\/code>: Calculate clustered standard errors<\/li>\nget.conf.clust<\/code>: Calculate confidence intervals for clustered standard errors<\/li>\n<\/ul>\n
\nCalculate heteroskedastically-robust standard errors<\/i><\/p>\nmod: the result of a call to lm()\n e.g., lm(satell ~ width + factor(color))\nvar: a list of indexes to extract only certain coefficients (default: all)\n e.g., 1:2 [to drop FE from above]\nestimator: an estimator type passed to vcovHC (default: White's estimator)\nReturns an object of type coeftest, with estimate, std. err, t and p values\n<\/pre>\n
\nCalculate clustered standard errors<\/i><\/p>\nmod: the result of a call to lm()\n e.g., lm(satell ~ width + factor(color))\ncluster: a cluster variable, with a length equal to the observation count\n e.g., color\nvar: a list of indexes to extract only certain coefficients (default: all)\n e.g., 1:2 [to drop FE from above]\nReturns an object of type coeftest, with estimate, std. err, t and p values\n<\/pre>\n
\nCalculate confidence intervals for clustered standard errors<\/i><\/p>\nmod: the result of a call to lm()\n e.g., lm(satell ~ width + factor(color))\ncluster: a cluster variable, with a length equal to the observation count\n e.g., color\nxx: new values for each of the variables (only needs to be as long as 'var')\n e.g., seq(0, 50, length.out=100)\nalpha: the level of confidence, as an error rate\n e.g., .05 [for 95% confidence intervals]\nvar: a list of indexes to extract only certain coefficients (default: all)\n e.g., 1:2 [to drop FE from above]\nReturns a data.frame of yhat, lo, and hi\n<\/pre>\n
Example in Stata and R<\/h3>\n
\nxi: reg satell width i.color, vce(hc2)
\nxi: reg satell width i.color, vce(cluster color)<\/p>\n\n------------------------------------------------------------------------------\n satell | Coef. Std. Err. t P>|t| [95% Conf. Interval]\n-------------+----------------------------------------------------------------\n width | .4633951 .1117381 4.15 0.000 .2428033 .6839868\n _cons | -8.412887 3.133074 -2.69 0.008 -14.59816 -2.227618\n------------------------------------------------------------------------------\n<\/pre>\n
\n------------------------------------------------------------------------------\n | Robust HC2\n satell | Coef. Std. Err. t P>|t| [95% Conf. Interval]\n-------------+----------------------------------------------------------------\n width | .4633951 .1028225 4.51 0.000 .2604045 .6663856\n _cons | -8.412887 2.939015 -2.86 0.005 -14.21505 -2.610726\n------------------------------------------------------------------------------\n<\/pre>\n
\n------------------------------------------------------------------------------\n | Robust\n satell | Coef. Std. Err. t P>|t| [95% Conf. Interval]\n-------------+----------------------------------------------------------------\n width | .4633951 .0466564 9.93 0.002 .3149135 .6118767\n _cons | -8.412887 1.258169 -6.69 0.007 -12.41694 -4.408833\n------------------------------------------------------------------------------\n<\/pre>\n
\nsource("tools_reg.R")\ntbl <- read.csv("data.csv")\n\nmod <- lm(satell ~ width + factor(color), data=tbl)\n\nsummary(mod)\nget.coef.robust(mod, estimator="HC2")\nget.coef.clust(mod, tbl$color)\n<\/pre>\n\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) -8.4129 3.1331 -2.685 0.00798 ** \nwidth 0.4634 0.1117 4.147 5.33e-05 ***\n<\/pre>\n
-8.412887 + .4633951 w<\/code> and R gives us -8.4129 + 0.4634 w<\/code>. These are identical, but with different reporting precision, as are the standard errors.<\/p>\n\n(Intercept) -8.41289 2.93902 -2.8625 0.004739 ** \nwidth 0.46340 0.10282 4.5067 1.229e-05 ***\n<\/pre>\n
\n Estimate Std. Error t value Pr(>|t|) \n(Intercept) -8.412887 1.258169 -6.6866 3.235e-10 ***\nwidth 0.463395 0.046656 9.9321 < 2.2e-16 ***\n<\/pre>\n