Estimation of one-way and multi-way clustered covariance matrices using an object-oriented approach.

vcovCL(x, cluster = NULL, type = NULL, sandwich = TRUE, fix = FALSE, ...)
meatCL(x, cluster = NULL, type = NULL, cadjust = TRUE, multi0 = FALSE, ...)



a fitted model object.


a variable indicating the clustering of observations, a list (or data.frame) thereof, or a formula specifying which variables from the fitted model should be used (see examples). By default (cluster = NULL), either attr(x, "cluster") is used (if any) or otherwise every observation is assumed to be its own cluster.


a character string specifying the estimation type (HC0--HC3). The default is to use "HC1" for lm objects and "HC0" otherwise.


logical. Should the sandwich estimator be computed? If set to FALSE only the meat matrix is returned.


logical. Should the covariance matrix be fixed to be positive semi-definite in case it is not?


logical. Should a cluster adjustment be applied?


logical. Should the HC0 estimate be used for the final adjustment in multi-way clustered covariances?


arguments passed to meatCL.


Clustered sandwich estimators are used to adjust inference when errors are correlated within (but not between) clusters. vcovCL allows for clustering in arbitrary many cluster dimensions (e.g., firm, time, industry), given all dimensions have enough clusters (for more details, see Cameron et al. 2011). If each observation is its own cluster, the clustered sandwich collapses to the basic sandwich covariance.

The function meatCL is the work horse for estimating the meat of clustered sandwich estimators. vcovCL is a wrapper calling sandwich and bread (Zeileis 2006). vcovCL is applicable beyond lm or glm class objects.

bread and meat matrices are multiplied to construct clustered sandwich estimators. The meat of a clustered sandwich estimator is the cross product of the clusterwise summed estimating functions. Instead of summing over all individuals, first sum over cluster.

A two-way clustered sandwich estimator \(M\) (e.g., for cluster dimensions "firm" and "industry" or "id" and "time") is a linear combination of one-way clustered sandwich estimators for both dimensions (\(M_{id}, M_{time}\)) minus the clustered sandwich estimator, with clusters formed out of the intersection of both dimensions (\(M_{id \cap time}\)): $$M = M_{id} + M_{time} - M_{id \cap time}$$ Additionally, each of the three terms can be weighted by the corresponding cluster bias adjustment factor (see below and Equation 20 in Zeileis et al. 2020). Instead of subtracting \(M_{id \cap time}\) as the last subtracted matrix, Ma (2014) suggests to subtract the basic HC0 covariance matrix when only a single observation is in each intersection of \(id\) and \(time\). Set multi0 = TRUE to subtract the basic HC0 covariance matrix as the last subtracted matrix in multi-way clustering. For details, see also Petersen (2009) and Thompson (2011).

With the type argument, HC0 to HC3 types of bias adjustment can be employed, following the terminology used by MacKinnon and White (1985) for heteroscedasticity corrections. HC0 applies no small sample bias adjustment. HC1 applies a degrees of freedom-based correction, \((n-1)/(n-k)\) where \(n\) is the number of observations and \(k\) is the number of explanatory or predictor variables in the model. HC1 is the most commonly used approach, and is the default, though it is less effective than HC2 and HC3 when the number of clusters is relatively small (Cameron et al. 2008). HC2 and HC3 types of bias adjustment are geared towards the linear model, but they are also applicable for GLMs (see Bell and McCaffrey 2002, and Kauermann and Carroll 2001, for details). A precondition for HC2 and HC3 types of bias adjustment is the availability of a hat matrix (or a weighted version therof for GLMs) and hence these two types are currently only implemented for lm and glm objects.

The cadjust argument allows to switch the cluster bias adjustment factor \(G/(G-1)\) on and off (where \(G\) is the number of clusters in a cluster dimension \(g\)) See Cameron et al. (2008) and Cameron et al. (2011) for more details about small-sample modifications.

The cluster specification can be made in a number of ways: The cluster can be a single variable or a list/data.frame of multiple clustering variables. If expand.model.frame works for the model object x, the cluster can also be a formula. By default (cluster = NULL), attr(x, "cluster") is checked and used if available. If not, every observation is assumed to be its own cluster. If the number of observations in the model x is smaller than in the original data due to NA processing, then the same NA processing can be applied to cluster if necessary (and x$na.action being available).

Cameron et al. (2011) observe that sometimes the covariance matrix is not positive-semidefinite and recommend to employ the eigendecomposition of the estimated covariance matrix, setting any negative eigenvalue(s) to zero. This fix is applied, if necessary, when fix = TRUE is specified.


A matrix containing the covariance matrix estimate.


Bell RM, McCaffrey DF (2002). “Bias Reduction in Standard Errors for Linear Regression with Multi-Stage Samples”, Survey Methodology, 28(2), 169--181.

Cameron AC, Gelbach JB, Miller DL (2008). “Bootstrap-Based Improvements for Inference with Clustered Errors”, The Review of Economics and Statistics, 90(3), 414--427. doi: 10.3386/t0344

Cameron AC, Gelbach JB, Miller DL (2011). “Robust Inference with Multiway Clustering”, Journal of Business & Ecomomic Statistics, 29(2), 238--249. doi: 10.1198/jbes.2010.07136

Kauermann G, Carroll RJ (2001). “A Note on the Efficiency of Sandwich Covariance Matrix Estimation”, Journal of the American Statistical Association, 96(456), 1387--1396. doi: 10.1198/016214501753382309

Ma MS (2014). “Are We Really Doing What We Think We Are Doing? A Note on Finite-Sample Estimates of Two-Way Cluster-Robust Standard Errors”, Mimeo, Availlable at SSRN: URL

MacKinnon, JG, White, H (1985). “Some Heteroskedasticity-Consistent Covariance Matrix Estimators with Improved Finite Sample Properties” Journal of Econometrics, 29(3), 305--325. doi: 10.1016/0304-4076(85)90158-7

Petersen MA (2009). “Estimating Standard Errors in Finance Panel Data Sets: Comparing Approaches”, The Review of Financial Studies, 22(1), 435--480. doi: 10.1093/rfs/hhn053

Thompson SB (2011). “Simple Formulas for Standard Errors That Cluster by Both Firm and Time”, Journal of Financial Economics, 99(1), 1--10. doi: 10.1016/j.jfineco.2010.08.016

Zeileis A (2004). “Econometric Computing with HC and HAC Covariance Matrix Estimator”, Journal of Statistical Software, 11(10), 1--17. doi: 10.18637/jss.v011.i10

Zeileis A (2006). “Object-Oriented Computation of Sandwich Estimators”, Journal of Statistical Software, 16(9), 1--16. doi: 10.18637/jss.v016.i09

Zeileis A, Köll S, Graham N (2020). “Various Versatile Variances: An Object-Oriented Implementation of Clustered Covariances in R.” Journal of Statistical Software, 95(1), 1--36. doi: 10.18637/jss.v095.i01

See also


## Petersen's data data("PetersenCL", package = "sandwich") m <- lm(y ~ x, data = PetersenCL) ## clustered covariances ## one-way vcovCL(m, cluster = ~ firm)
#> (Intercept) x #> (Intercept) 4.490702e-03 -6.473517e-05 #> x -6.473517e-05 2.559927e-03
vcovCL(m, cluster = PetersenCL$firm) ## same
#> (Intercept) x #> (Intercept) 4.490702e-03 -6.473517e-05 #> x -6.473517e-05 2.559927e-03
## one-way with HC2 vcovCL(m, cluster = ~ firm, type = "HC2")
#> (Intercept) x #> (Intercept) 4.494487e-03 -6.592912e-05 #> x -6.592912e-05 2.568236e-03
## two-way vcovCL(m, cluster = ~ firm + year)
#> (Intercept) x #> (Intercept) 4.233313e-03 -2.845344e-05 #> x -2.845344e-05 2.868462e-03
vcovCL(m, cluster = PetersenCL[, c("firm", "year")]) ## same
#> (Intercept) x #> (Intercept) 4.233313e-03 -2.845344e-05 #> x -2.845344e-05 2.868462e-03
## comparison with cross-section sandwiches ## HC0 all.equal(sandwich(m), vcovCL(m, type = "HC0", cadjust = FALSE))
#> [1] TRUE
## HC2 all.equal(vcovHC(m, type = "HC2"), vcovCL(m, type = "HC2"))
#> [1] TRUE
## HC3 all.equal(vcovHC(m, type = "HC3"), vcovCL(m, type = "HC3"))
#> [1] TRUE
## Innovation data data("InstInnovation", package = "sandwich") ## replication of one-way clustered standard errors for model 3, Table I ## and model 1, Table II in Berger et al. (2017), see ?InstInnovation ## count regression formula f1 <- cites ~ institutions + log(capital/employment) + log(sales) + industry + year ## model 3, Table I: Poisson model ## one-way clustered standard errors tab_I_3_pois <- glm(f1, data = InstInnovation, family = poisson) vcov_pois <- vcovCL(tab_I_3_pois, InstInnovation$company) sqrt(diag(vcov_pois))[2:4]
#> institutions log(capital/employment) log(sales) #> 0.002406388 0.135953255 0.041523405
## coefficient tables if(require("lmtest")) { coeftest(tab_I_3_pois, vcov = vcov_pois)[2:4, ] }
#> Estimate Std. Error z value Pr(>|z|) #> institutions 0.009687237 0.002406388 4.025634 5.682195e-05 #> log(capital/employment) 0.482883549 0.135953255 3.551835 3.825545e-04 #> log(sales) 0.820317600 0.041523405 19.755548 7.187199e-87
if (FALSE) { ## model 1, Table II: negative binomial hurdle model ## (requires "pscl" or alternatively "countreg" from R-Forge) library("pscl") library("lmtest") tab_II_3_hurdle <- hurdle(f1, data = InstInnovation, dist = "negbin") # dist = "negbin", zero.dist = "negbin", separate = FALSE) vcov_hurdle <- vcovCL(tab_II_3_hurdle, InstInnovation$company) sqrt(diag(vcov_hurdle))[c(2:4, 149:151)] coeftest(tab_II_3_hurdle, vcov = vcov_hurdle)[c(2:4, 149:151), ] }