Standard error clustering in R (either manually or in plm)

I am trying to understand standard error "clustering" and how to execute in R (it is trivial in Stata). In R I have been unsuccessful using either plm or writing my own function. I'll use the diamonds data from the ggplot2 package.

I can do fixed effects with either dummy variables

> library(plyr)
> library(ggplot2)
> library(lmtest)
> library(sandwich)
> # with dummies to create fixed effects
> fe.lsdv <- lm(price ~ carat + factor(cut) + 0, data = diamonds)
> ct.lsdv <- coeftest(fe.lsdv, vcov. = vcovHC)
> ct.lsdv

t test of coefficients:

Estimate Std. Error  t value  Pr(>|t|)
carat                 7871.082     24.892  316.207 < 2.2e-16 ***
factor(cut)Fair      -3875.470     51.190  -75.707 < 2.2e-16 ***
factor(cut)Good      -2755.138     26.570 -103.692 < 2.2e-16 ***
factor(cut)Very Good -2365.334     20.548 -115.111 < 2.2e-16 ***
factor(cut)Premium   -2436.393     21.172 -115.075 < 2.2e-16 ***
factor(cut)Ideal     -2074.546     16.092 -128.920 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

or by de-meaning both left- and right-hand sides (no time invariant regressors here) and correcting degrees of freedom.

> # by demeaning with degrees of freedom correction
> diamonds <- ddply(diamonds, .(cut), transform, price.dm = price - mean(price), carat.dm = carat  .... [TRUNCATED]
> fe.dm <- lm(price.dm ~ carat.dm + 0, data = diamonds)
> ct.dm <- coeftest(fe.dm, vcov. = vcovHC, df = nrow(diamonds) - 1 - 5)
> ct.dm

t test of coefficients:

Estimate Std. Error t value  Pr(>|t|)
carat.dm 7871.082     24.888  316.26 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

I can't replicate these results with plm, because I don't have a "time" index (i.e., this isn't really a panel, just clusters that could have a common bias in their error terms).

> plm.temp <- plm(price ~ carat, data = diamonds, index = "cut")
duplicate couples (time-id)
Error in pdim.default(index[], index[]) :

I also tried to code my own covariance matrix with clustered standard error using Stata's explanation of their cluster option (explained here), which is to solve $$\hat V_{cluster} = (X'X)^{-1} \left( \sum_{j=1}^{n_c} u_j'u_j \right) (X'X)^{-1}$$ where $u_j = \sum_{cluster~j} e_i * x_i$, $n_c$ si the number of clusters, $e_i$ is the residual for the $i^{th}$ observation and $x_i$ is the row vector of predictors, including the constant (this also appears as equation (7.22) in Wooldridge's Cross Section and Panel Data). But the following code gives very large covariance matrices. Are these very large values given the small number of clusters I have? Given that I can't get plm to do clusters on one factor, I'm not sure how to benchmark my code.

> # with cluster robust se
> lm.temp <- lm(price ~ carat + factor(cut) + 0, data = diamonds)
>
> # using the model that Stata uses
> stata.clustering <- function(x, clu, res) {
+     x <- as.matrix(x)
+     clu <- as.vector(clu)
+     res <- as.vector(res)
+     fac <- unique(clu)
+     num.fac <- length(fac)
+     num.reg <- ncol(x)
+     u <- matrix(NA, nrow = num.fac, ncol = num.reg)
+     meat <- matrix(NA, nrow = num.reg, ncol = num.reg)
+
+     # outer terms (X'X)^-1
+     outer <- solve(t(x) %*% x)
+
+     # inner term sum_j u_j'u_j where u_j = sum_i e_i * x_i
+     for (i in seq(num.fac)) {
+         index.loop <- clu == fac[i]
+         res.loop <- res[index.loop]
+         x.loop <- x[clu == fac[i], ]
+         u[i, ] <- as.vector(colSums(res.loop * x.loop))
+     }
+     inner <- t(u) %*% u
+
+     #
+     V <- outer %*% inner %*% outer
+     return(V)
+ }
> x.temp <- data.frame(const = 1, diamonds[, "carat"])
> summary(lm.temp)

Call:
lm(formula = price ~ carat + factor(cut) + 0, data = diamonds)

Residuals:
Min       1Q   Median       3Q      Max
-17540.7   -791.6    -37.6    522.1  12721.4

Coefficients:
Estimate Std. Error t value Pr(>|t|)
carat                 7871.08      13.98   563.0   <2e-16 ***
factor(cut)Fair      -3875.47      40.41   -95.9   <2e-16 ***
factor(cut)Good      -2755.14      24.63  -111.9   <2e-16 ***
factor(cut)Very Good -2365.33      17.78  -133.0   <2e-16 ***
factor(cut)Premium   -2436.39      17.92  -136.0   <2e-16 ***
factor(cut)Ideal     -2074.55      14.23  -145.8   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1511 on 53934 degrees of freedom
Multiple R-squared: 0.9272, Adjusted R-squared: 0.9272
F-statistic: 1.145e+05 on 6 and 53934 DF,  p-value: < 2.2e-16

> stata.clustering(x = x.temp, clu = diamonds$cut, res = lm.temp$residuals)
const diamonds....carat..
const                11352.64           -14227.44
diamonds....carat.. -14227.44            17830.22

Can this be done in R? It is a fairly common technique in econometrics (there's a brief tutorial in this lecture), but I can't figure it out in R. Thanks!

• @ricardh, curse upon all economists for not checking whether the term they want to use is already used in statistics. Cluster in this context means group and is totally unrelated to cluster analysis, this is why rseek gave you unrelated results. Hence I removed the clustering tag. For panel data analysis check out package plm. It has a nice vignette, so you may find what you want. As for your question it is not clear what you want. Within group standard errors? – mpiktas Apr 27 '11 at 7:03
• @ricardh, it would help a lot if you could link to some manual of Stata where this cluster option is explained. I am sure it would be possible to replicate in R. – mpiktas Apr 27 '11 at 7:04
• +1 for that comment. economists colonize terminology like crazy. Though sometimes it's hard to pick the villain. Ii took a while e.g. until I realized factor had nothing to do with factanal but refers to categorized variables. However cluster in R refers to cluster analysis, k-means is just THE partitioning method: statmethods.net/advstats/cluster.html . I don't get your question, but I also guess cluster has nothing to do with it. – hans0l0 Apr 27 '11 at 9:17
• @mpiktas, @ran2 -- Thanks! I hope I clarified the question. In short, why does "standard error clustering" exist if it is just fixed effects, which already existed? – Richard Herron Apr 27 '11 at 9:33
• cluster.vcov function in the "multiwayvcov" package does what you are looking for. – user84051 Aug 4 '15 at 2:59

For White standard errors clustered by group with the plm framework try

coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0",cluster="group"))

where model.plm is a plm model.

http://www.inside-r.org/packages/cran/plm/docs/vcovHC or the plm package documentation

EDIT:

For two-way clustering (e.g. group and time) see the following link:

http://people.su.se/~ma/clustering.pdf

Here is another helpful guide for the plm package specifically that explains different options for clustered standard errors:

http://www.princeton.edu/~otorres/Panel101R.pdf

Clustering and other information, especially for Stata, can be found here:

http://www.kellogg.northwestern.edu/faculty/petersen/htm/papers/se/se_programming.htm

EDIT 2:

Here are examples that compare R and stata: http://www.richard-bluhm.com/clustered-ses-in-r-and-stata-2/

Also, the multiwayvcov may be helpful. This post provides a helpful overview: http://rforpublichealth.blogspot.dk/2014/10/easy-clustered-standard-errors-in-r.html

From the documentation:

library(multiwayvcov)
library(lmtest)
data(petersen)
m1 <- lm(y ~ x, data = petersen)

# Cluster by firm
vcov_firm <- cluster.vcov(m1, petersen$firmid) coeftest(m1, vcov_firm) # Cluster by year vcov_year <- cluster.vcov(m1, petersen$year)
coeftest(m1, vcov_year)
# Cluster by year using a formula
vcov_year_formula <- cluster.vcov(m1, ~ year)
coeftest(m1, vcov_year_formula)

# Double cluster by firm and year
vcov_both <- cluster.vcov(m1, cbind(petersen$firmid, petersen$year))
coeftest(m1, vcov_both)
# Double cluster by firm and year using a formula
vcov_both_formula <- cluster.vcov(m1, ~ firmid + year)
coeftest(m1, vcov_both_formula)
• for me coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0")) as well as coeftest(model.plm, vcov=vcovHC(model.plm,type="HC0",cluster="group")) produce identical results. Do you know why this is the case? – Peter Pan Dec 17 '17 at 15:30
• The link people.su.se/~ma/clustering.pdf is no longer working. Do you remember the title of the page? – MERose Jan 16 '18 at 20:50

After a lot of reading, I found the solution for doing clustering within the lm framework.

There's an excellent white paper by Mahmood Arai that provides a tutorial on clustering in the lm framework, which he does with degrees-of-freedom corrections instead of my messy attempts above. He provides his functions for both one- and two-way clustering covariance matrices here.

Finally, although the content isn't available free, Angrist and Pischke's Mostly Harmless Econometrics has a section on clustering that was very helpful.

Update on 4/27/2015 to add code from blog post.

api=read.csv("api.csv") #create the variable api from the corresponding csv
attach(api) # attach of data.frame objects
api1=api[c(1:6,8:310),] # one missing entry in row nr. 7
modell.api=lm(API00 ~ GROWTH + EMER + YR_RND, data=api1) # creation of a simple linear model for API00 using the regressors Growth, Emer and Yr_rnd.

##creation of the function according to Arai:
clx <- function(fm, dfcw, cluster) {
library(sandwich)
library(lmtest)
library(zoo)
M <- length(unique(cluster))
N <- length(cluster)
dfc <- (M/(M-1))*((N-1)/(N-fm$rank)) # anpassung der freiheitsgrade u <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum)) vcovCL <-dfc * sandwich (fm, meat = crossprod(u)/N) * dfcw coeftest(fm, vcovCL) } clx(modell.api, 1, api1$DNUM) #creation of results.
• Arai's paper is no longer online. Can you provide the actual link? – MERose Apr 27 '15 at 11:48
• @MERose -- Sorry! Unfortunately we can't attach pdfs! I found this blog post that benchmarks the code. I will edit this answer to include the code. – Richard Herron Apr 27 '15 at 16:09
• This might be an updated version of the white paper: Mahmood Arai, Cluster-robust standard errors using R. – gung Dec 12 '16 at 0:31

The easiest way to compute clustered standard errors in R is to use the modified summary function.

lm.object <- lm(y ~ x, data = data)
summary(lm.object, cluster=c("c"))

There's an excellent post on clustering within the lm framework. The site also provides the modified summary function for both one- and two-way clustering. You can find the function and the tutorial here.

If you don't have a time index, you don't need one: plm will add a fictitious one by itself, and it won't be used unless you ask for it. So this call should work:

> x <- plm(price ~ carat, data = diamonds, index = "cut")
Error in pdim.default(index[], index[]) :
duplicate couples (time-id)

Except that it doesn't, which suggests you've hit a bug in plm. (This bug has now been fixed in SVN. You can install the development version from here.)

But since this would be a fictitious time index anyway, we can create it by ourselves:

diamonds\$ftime <- 1:NROW(diamonds)  ##fake time

Now this works:

x <- plm(price ~ carat, data = diamonds, index = c("cut", "ftime"))
coeftest(x, vcov.=vcovHC)
##
## t test of coefficients:
##
##       Estimate Std. Error t value  Pr(>|t|)
## carat  7871.08     138.44  56.856 < 2.2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Important note: vcovHC.plm() in plm will by default estimate Arellano clustered by group SEs. Which is different from what vcovHC.lm() in sandwich will estimate (e.g. the vcovHC SEs in the original question), that is heteroscedasticity-consistent SEs with no clustering.

A separate approach to this is sticking to lm dummy variable regressions and the multiwayvcov package.

library("multiwayvcov")
fe.lsdv <- lm(price ~ carat + factor(cut) + 0, data = diamonds)
coeftest(fe.lsdv, vcov.= function(y) cluster.vcov(y, ~ cut, df_correction = FALSE))
##
## t test of coefficients:
##
##                      Estimate Std. Error t value  Pr(>|t|)
## carat                 7871.08     138.44  56.856 < 2.2e-16 ***
## factor(cut)Fair      -3875.47     144.83 -26.759 < 2.2e-16 ***
## factor(cut)Good      -2755.14     117.56 -23.436 < 2.2e-16 ***
## factor(cut)Very Good -2365.33     111.63 -21.188 < 2.2e-16 ***
## factor(cut)Premium   -2436.39     123.48 -19.731 < 2.2e-16 ***
## factor(cut)Ideal     -2074.55      97.30 -21.321 < 2.2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

In both cases you will get the Arellano (1987) SEs with clustering by group. The multiwayvcov package is a direct and significant evolution of Arai's original clustering functions.

You can also look at the resulting variance-covariance matrix from both approaches, yielding the same variance estimate for carat:

vcov.plm <- vcovHC(x)
vcov.lsdv <- cluster.vcov(fe.lsdv, ~ cut, df_correction = FALSE)
vcov.plm
##          carat
## carat 19165.28
diag(vcov.lsdv)
##                carat      factor(cut)Fair      factor(cut)Good factor(cut)Very Good   factor(cut)Premium     factor(cut)Ideal
##            19165.283            20974.522            13820.365            12462.243            15247.584             9467.263