In many scenarios, data are structured in groups or clusters, e.g. pupils within classes (within schools), survey respondents within countries or, for longitudinal surveys, survey answers per subject. Simply ignoring this structure will likely lead to spuriously low standard errors, i.e. a misleadingly precise estimate of our coefficients. This in turn leads to overlynarrow confidence intervals, overlylow pvalues and possibly wrong conclusions.
Clustered standard errors are a common way to deal with this problem. Unlike Stata, R doesn’t have builtin functionality to estimate clustered standard errors. There are several packages though that add this functionality and this article will introduce three of them, explaining how they can be used and what their advantages and disadvantages are. Before that, I will outline the theory behind (clustered) standard errors for linear regression. The last section is used for a performance comparison between the three presented packages. If you’re already familiar with the concept of clustered standard errors, you may skip to the handson part right away.
Data
We’ll work with the dataset nlswork that’s included in Stata, so we can easily compare the results with Stata. The data comes from the US National Longitudinal Survey (NLS) and contains information about more than 4,000 young working women. As for this example, we’re interested in the relationship between wage (here as logscaled GNPadjusted wage) as dependent variable (DV) ln_wage
and survey participant’s current age
, job tenure
in years and union
membership as independent variables. It’s a longitudinal survey, so subjects were asked repeatedly between 1968 and 1988 and each subject is identified by an unique idcode
.
The example data is used for illustrative purposes only and we skip many things that we’d normally do, such as investigating descriptive statistics and exploratory plots. To keep the data size limited, we’ll only work with a subset of the data (only subjects with IDs 1 to 100) and we also simply dismiss any observations that contain missing values.
library(webuse)
library(dplyr)
nlswork_orig < webuse('nlswork')
nlswork < filter(nlswork_orig, idcode <= 100) %>%
select(idcode, year, ln_wage, age, tenure, union) %>%
filter(complete.cases(.)) %>%
mutate(union = as.integer(union),
idcode = as.factor(idcode))
str(nlswork)
tibble [386 Ă— 6] (S3: tbl_df/tbl/data.frame) $ idcode : Factor w/ 82 levels "1","2","3","4",..: 1 1 1 1 1 1 1 2 2 2 â€¦ $ year : num [1:386] 72 77 80 83 85 87 88 71 77 78 â€¦ $ ln_wage: num [1:386] 1.59 1.78 2.55 2.42 2.61 â€¦ $ age : num [1:386] 20 25 28 31 33 35 37 19 25 26 â€¦ $ tenure : num [1:386] 0.917 1.5 1.833 0.667 1.917 â€¦ $ union : int [1:386] 1 0 1 1 1 1 1 0 1 1 â€¦
Let’s have a look at the first few observations. They contain data from subject #1, who was surveyed several times between 1972 and 1988, and a few observations from subject #2.
idcode year ln_wage age tenure union 1 72 1.59 20 0.917 1 1 77 1.78 25 1.5 0 1 80 2.55 28 1.83 1 1 83 2.42 31 0.667 1 1 85 2.61 33 1.92 1 1 87 2.54 35 3.92 1 1 88 2.46 37 5.33 1 2 71 1.36 19 0.25 0 2 77 1.73 25 2.67 1 2 78 1.69 26 3.67 1
A summary of all but the idcode variable gives the following:
year ln_wage age tenure union Min. :70.00 Min. :0.4733 Min. :18.0 Min. : 0.000 Min. :0.0000 1st Qu.:73.00 1st Qu.:1.6131 1st Qu.:25.0 1st Qu.: 1.167 1st Qu.:0.0000 Median :80.00 Median :1.9559 Median :31.0 Median : 2.417 Median :0.0000 Mean :79.61 Mean :1.9453 Mean :30.8 Mean : 3.636 Mean :0.2591 3rd Qu.:85.00 3rd Qu.:2.2349 3rd Qu.:36.0 3rd Qu.: 4.958 3rd Qu.:1.0000 Max. :88.00 Max. :3.5791 Max. :45.0 Max. :19.000 Max. :1.0000
We have 82 subjects in our subset and the number of times each subject was surveyed ranges from only once to twelve times:
summary(as.integer(table(nlswork$idcode)))
Min. 1st Qu. Median Mean 3rd Qu. Max. 1.000 2.000 4.000 4.707 7.000 12.000
In more than one quarter of the observations, the subject answered to be currently member of a trade union.
The following shows the distribution of the DV in our data. It is roughly normally distributed with mean 1.95 and standard deviation (SD) 0.45.
We can calculate the mean and SD of the DV separately for each subject. A histogram of these subjectspecific means reveals more variability:
y_mean_sd_cl < sapply(levels(nlswork$idcode), function(idcode) {
y_cl < nlswork$ln_wage[nlswork$idcode == idcode]
c(mean(y_cl), sd(y_cl))
})
hist(y_mean_sd_cl[1,], breaks = 20,
main = 'Histogram of DV means per subject', xlab = NA)
We can compare the SD of the subjectspecific means with the mean of the SDs calculated from each subjects’ repeated measures.
c(sd(y_mean_sd_cl[1,]), mean(y_mean_sd_cl[2,], na.rm = TRUE))
[1] 0.4038449 0.2221142
The SD between the subjectspecific means is almost twice as large as the mean of the SD from each subjects’ values. This shows that there’s much more variability between each subject than within each subject’s repeated measures regarding the DV.
Fixedeffects model, not adjusting for clustered observations
Our data contains repeated measures for each subject, so we have panel data in which each subject forms a group or cluster. We can use a fixedeffects (FE) model to account for unobserved subjectspecific characteristics. We do so by including the subject’s idcode
in our model formula. It’s important to note that idcode
is of type factor (we applied idcode = as.factor(idcode)
when we prepared the data) so that for each factor level (i.e. each subject) an FE coefficient will be estimated that represents the subjectspecific mean of our DV.^{ 1}
Let’s specify and fit such a model using lm
. We include job tenure
, union
membership and an interaction between both (the latter mainly for illustrative purposes later when we estimate marginal effects). We also control for age
and add idcode
as FE variable.
m1 < lm(ln_wage ~ age + tenure + union + tenure:union + idcode,
data = nlswork)
summary(m1)
Call: lm(formula = ln_wage ~ age + tenure + union + tenure:union + idcode, data = nlswork) Residuals: Min 1Q Median 3Q Max 0.96463 0.09405 0.00000 0.11460 1.23525 Coefficients: Estimate Std. Error t value Pr(>t) (Intercept) 1.882e+00 1.314e01 14.325 < 2e16 *** age 5.631e03 3.110e03 1.811 0.071193 . tenure 2.076e02 6.964e03 2.980 0.003115 ** union 1.746e01 6.065e02 2.879 0.004272 ** idcode2 6.174e01 1.285e01 4.803 2.47e06 *** idcode3 5.625e01 1.329e01 4.234 3.05e05 *** [...]
We’re not really interested in the subjectspecific means (the FE coefficients), so let’s filter them out and only show our coefficients of interest:
m1coeffs_std < data.frame(summary(m1)$coefficients)
coi_indices < which(!startsWith(row.names(m1coeffs_std), 'idcode'))
m1coeffs_std[coi_indices,]
Estimate Std..Error t.value Pr...t.. (Intercept) 1.882478232 0.131411504 14.325064 8.022367e36 age 0.005630809 0.003109803 1.810664 7.119315e02 tenure 0.020756426 0.006964417 2.980353 3.114742e03 union 0.174619394 0.060646038 2.879321 4.272027e03 tenure:union 0.014974113 0.009548509 1.568215 1.178851e01
Unsurprisingly, job tenure and especially union membership are positively associated with wage. The coefficient of the interaction term shows that with union membership the job tenure effect is even a bit higher, though not significantly.
In the next two sections we’ll see how standard errors for our estimates are usually computed and how this fits into a framework called “sandwich estimators.” Using this framework, we’ll see how the standard error calculations can be adjusted for clustered data.
Standard errors
In ordinary least squares (OLS) regression, we assume that the regression model errors are independent. This is not the case here: Each subject may be surveyed several times so within each subject’s repeated measures, the errors will be correlated. Although that is not a problem for our regression estimates (they are still unbiased [Roberts 2013]), it is a problem for for the precision of our estimates — the precision will typically be overestimated, i.e. the standard errors (SEs) will be lower than they should be [Cameron and Miller 2013]. The intuition behind this regarding our example is that within our clusters we usually have lower variance since the answers come from the same subject and are correlated. This lowers our estimates’ SEs.
We can deal with this using clustered standard errors with subjects representing our clusters. But before we do this, let’s first have a closer look on how “classic” OLS estimates’ SEs are actually computed.
In matrix notation, a linear model has the form
This model has parameters (including the intercept parameter ) expressed as parameter vector . The parameters will be estimated from observations in our data. The DV is (an vector), the independent variables form an matrix . Finally, the error term is an vector that captures everything that influences but cannot be explained by . In classic OLS, we assume that has a mean of zero and a variance of .
By minimizing an estimation for our parameters, , can be found as . The estimated variance for these parameters is then
(1)
where is the estimated variance of the error or residual variance. This is calculated as
with being the residuals. The numerator is also called the residual sum of squares and the denominator is the degrees of freedom.
You can see for example in Sheather 2009 for how the formulas for and are derived or watch this video for a clear stepbystep derivation.
We now have all the pieces together to replicate the standard errors from model m1
with our own calculations. To translate these formulae to R, we use model.matrix
to get the design matrix , residuals
for the residual vector , nobs
for the number of observations , ncol(X)
for the number or parameters, solve
to calculate the inverse of and diag
to extract the diagonal of a square matrix.
X < model.matrix(m1)
u < residuals(m1)
n < nobs(m1)
p < ncol(X)
sigma2 < sum(u^2) / (n  p)
# solve (X^T X) A = I, where I is identity matrix > A is (X^T X)^1
crossXinv < solve(t(X) %*% X, diag(p))
m1se < sqrt(diag(sigma2 * crossXinv))
m1se
[1] 0.131411504 0.003109803 0.006964417 0.060646038 0.128532897 ...
Let’s check if this is equal to the standard errors calculated by lm
(using near
because of minor deviations due to floating point number imprecision):
all(near(m1se, m1coeffs_std$Std..Error))
[1] TRUE
We extracted our parameter estimates’ variance from the diagonal of the (variance)covariance or vcov matrix and R has the vcov
function to calculate it from a fitted model. It’s exactly what we computed before using eq. 1:
all(near(sigma2 * crossXinv, vcov(m1)))
[1] TRUE
Clustered standard errors
Classic OLS SEs can be generalized so that some assumptions, namely that the regression model errors are independent, can be relaxed. The foundation for this is the sandwich estimator^{ 2}
(2)
Let’s first understand how the above equation relates to eq. 1, the classic OLS parameter variance: One assumption of classic OLS is constant variance (or homoscedasticity) in the errors across the full spectrum of our DV. This implicates that is a diagonal matrix with identical elements, i.e. . Plugging this into eq. 2 gives eq. 1 which shows that the classic OLS parameter variance is a specialized form of the sandwich estimator.^{ 3}
When we want to obtain clustered SEs, we need to consider that in the “meat” part of eq. 2 is not a diagonal matrix with identical elements anymore, hence this can’t be simplified to eq. 1. Instead, we can assume that is blockdiagonal with the clusters forming the blocks. This means we assume that the variance in the errors is constant only within clusters and so we first calculate per cluster and then sum the . Cameron and Miller 2013 (p. 11) show how is calculated in detail and also which finitesample correction factor is applied. From this article we get the equation
(3)
where is the number of clusters. It’s interesting to see how the residuals are added up per cluster and then averaged. As Cameron and Miller 2013 (p. 13) notes, this implicates an important limitation: With a low number of clusters, this averaging is imprecise.
Let’s translate this formula to R. We already have as u
(the residuals) and the design matrix X
. We can generate a list of matrices , sum them and multiply the correction factor:
omegaj < lapply(levels(nlswork$idcode), function(idcode) {
j < nlswork$idcode == idcode
# drop = FALSE: don't drop dimensions when we have only one obs.
X_j < X[j, , drop = FALSE]
# tcrossprod is outer product x * x^T
t(X_j) %*% tcrossprod(u[j]) %*% X_j
})
n_cl < length(levels(nlswork$idcode)) # num. clusters
# correction factor * sum of omega_j
omega < (n1) / (np) * (n_cl / (n_cl1)) * Reduce('+', omegaj)
# sandwich formula; extract diagonal and take square root to get SEs
m1clse < sqrt(diag(crossXinv %*% omega %*% crossXinv))
m1clse
[1] 0.157611390 0.006339777 0.011149190 0.101970509 0.020561516 ...
We will later check that this matches the estimates calculated with R packages that implement clustered SE estimation. For now, let’s compare the classic OLS SEs with the clustered SEs:
m1coeffs_with_clse < cbind(m1coeffs_std, ClustSE = m1clse)
m1coeffs_with_clse[coi_indices, c(1, 2, 5)]
Estimate Std..Error ClustSE (Intercept) 1.882478232 0.131411504 0.157611390 age 0.005630809 0.003109803 0.006339777 tenure 0.020756426 0.006964417 0.011149190 union 0.174619394 0.060646038 0.101970509 tenure:union 0.014974113 0.009548509 0.009646023
We can see that, as expected, the clustered SEs are all a bit higher than the classic OLS SEs.
The above calculations were used to show what’s happening “under the hood” and also how the formulas used for these calculations are motivated. However, doing these calculations “by hand” is errorprone and slow. It’s better to use well trusted packages for daily work and so next we’ll have a look at some of these packages and how they can be used. Still, it’s helpful to understand some background and the limitations for this approach. See Cameron and Miller 2013 for a much more thorough guide (though only with examples in Stata) that also considers topics like which variable(s) to use for clustering, what to do when dealing with a low number of clusters or how to implement multiway clustering.
Option 1: sandwich and lmtest
The sandwich package implements several methods for robust covariance estimators, including clustered SEs. Details are explained in Zeileis et al. 2020. The accompanying lmtest package provides functions for coefficient tests that take into account the calculated robust covariance estimates.
As explained initially, the parameter estimates from our model are fine despite the clustered structure of our data. But the SEs are likely biased downward and need to be corrected. This is why we can resume to work with our initially estimated model m1
from lm
. There’s no need to refit it and sandwich works with lm model objects (and also some other types of models such as some glm models). We only have to adjust how we test our coefficient estimates in the following way:
 We need to use
coeftest
from the lmtest package;  we need to pass it our model and either a function to calculate the covariance matrix or an already estimated covariance matrix to the
vcov
parameter;  we need to specify a cluster variable in the
cluster
parameter.
The sandwich package provides several functions for estimating robust covariance matrices. We need vcovCL
for clustered covariance estimation and will pass this function as vcov
parameter. Furthermore, we cluster by subject ID, so the cluster variable is idcode
.
library(sandwich)
library(lmtest)
m1coeffs_cl < coeftest(m1, vcov = vcovCL, cluster = ~idcode)
m1coeffs_cl[coi_indices,]
Estimate Std. Error t value Pr(>t) (Intercept) 1.882478232 0.157611390 11.9437956 3.667970e27 age 0.005630809 0.006339777 0.8881715 3.751601e01 tenure 0.020756426 0.011149190 1.8616981 6.362342e02 union 0.174619394 0.101970509 1.7124500 8.784708e02 tenure:union 0.014974113 0.009646023 1.5523613 1.216301e01
The calculated SE values seem familiar and they are indeed equal to what we calculated before as m1clse
“by hand”:
all(near(m1clse, m1coeffs_cl[,2]))
[1] TRUE
The lmtest package provides several functions for common postestimation tasks, for example coefci
to calculate confidence intervals (CIs). If we use these, we need to make sure to specify the same type of covariance estimation, again by passing the appropriate vcov
and cluster
parameters:
(m1cis < coefci(m1, parm = coi_indices, vcov = vcovCL,
cluster = ~idcode))
2.5 % 97.5 % (Intercept) 1.572314302 2.19264216 age 0.006845258 0.01810688 tenure 0.001184099 0.04269695 union 0.026048678 0.37528746 tenure:union 0.004008324 0.03395655
This is really important, as otherwise the classic (nonclustered) covariance estimation is applied by default. This, due to lower SEs, leads to narrower CIs:
coefci(m1, parm = coi_indices)
2.5 % 97.5 % (Intercept) 1.6238731375 2.14108333 age 0.0004889822 0.01175060 tenure 0.0070511278 0.03446172 union 0.0552738733 0.29396491 tenure:union 0.0038164264 0.03376465
Here, the tenure
and union
CIs suddenly don’t include zero any more!
Instead of passing vcovCL
as function to the vcov
parameter, it’s more convenient and computationally more efficient to calculate the covariance matrix only once using vcovCL
and then passing this matrix to functions like coeftest
and coefci
instead:
cl_vcov_mat < vcovCL(m1, cluster = ~idcode)
Now we pass this matrix for the vcov
parameter. We don’t need to specify the cluster
parameter anymore, since this information was only needed in the previous step.
m1coeffs_cl2 < coeftest(m1, vcov = cl_vcov_mat)
all(near(m1coeffs_cl[,2], m1coeffs_cl2[,2])) # same SEs?
[1] TRUE
m1cis2 < coefci(m1, parm = coi_indices, vcov = cl_vcov_mat)
all(near(m1cis, m1cis2)) # same CIs?
[1] TRUE
Another example would be to calculate marginal effects, for example with the margins package. Again, to arrive at clustered SEs we will need to pass the proper covariance matrix via the vcov
parameter. We do this for the marginal effect of tenure
at the two levels of union
:
library(margins)
margins(m1, vcov = cl_vcov_mat, variables = 'tenure',
at = list(union = 0:1)) %>%
summary()
factor union AME SE z p lower upper tenure 0.0000 0.0208 0.0111 1.8617 0.0626 0.0011 0.0426 tenure 1.0000 0.0357 0.0083 4.3089 0.0000 0.0195 0.0520
Otherwise classic SEs are estimated, which are smaller:
margins(m1, variables = 'tenure', at = list(union = 0:1)) %>% summary()
factor union AME SE z p lower upper tenure 0.0000 0.0208 0.0070 2.9804 0.0029 0.0071 0.0344 tenure 1.0000 0.0357 0.0081 4.3846 0.0000 0.0198 0.0517
As you can see, the combination of lm
and the packages sandwich and lmtest are all you need for estimating clustered SEs and inference. However, you really need to be careful to include the covariance matrix at all steps of your calculations.
Option 2: lm.cluster
from miceadds
There’s also lm.cluster
from the package miceadds^{ 4}, which may be a bit more convenient to use. Internally, it basically does the same that we’ve done before by employing sandwich’s vcovCL
(see source code parts here and there for example).
Instead of using lm
, we fit the model with lm.cluster
and specify a cluster variable (this time as string, not as formula). The model summary then contains the clustered SEs:
library(miceadds)
m2 < lm.cluster(ln_wage ~ age + tenure + union + tenure:union + idcode,
cluster = 'idcode',
data = nlswork)
m2coeffs < data.frame(summary(m2))
m2coeffs[!startsWith(row.names(m2coeffs), 'idcode'),]
Estimate Std..Error t.value Pr...t.. (Intercept) 1.882478232 0.157611390 11.9437956 6.995639e33 age 0.005630809 0.006339777 0.8881715 3.744485e01 tenure 0.020756426 0.011149190 1.8616981 6.264565e02 union 0.174619394 0.101970509 1.7124500 8.681378e02 tenure:union 0.014974113 0.009646023 1.5523613 1.205758e01
An object m
that is returned from lm.cluster
is a list that contains the lm
object as m$lmres
and the covariance matrix as m$vcov
. Again, these objects need to be “dragged along” if we want to do further computations. For margins
, we also need to pass the data again via data = nlswork
:
margins(m2$lm_res, vcov = m2$vcov, variables = 'tenure',
at = list(union = 0:1), data = nlswork) %>%
summary()
factor union AME SE z p lower upper tenure 0.0000 0.0208 0.0111 1.8617 0.0626 0.0011 0.0426 tenure 1.0000 0.0357 0.0083 4.3089 0.0000 0.0195 0.0520
The result is consistent with our former computations. The advantage over the “lm + sandwich + lmtest” approach is that you can do clustered SE estimation and inference in one go. For further calculations you still need to be careful to supply the covariance matrix from m$vcov
.
Option 3: lm_robust
from estimatr
Another option is to use lm_robust
from the estimatr package which is part of the DeclareDesign framework [Blair et al. 2019]. Like lm.cluster
, it’s more convenient to use, but it doesn’t rely on sandwich and lmtest in the background and instead comes with an own implementation for model fitting and covariance estimation. This implementation is supposed to be faster than the other approaches and we’ll check that for our example later.
But first, let’s fit a model with clustered SEs using lm_robust
. We use the same formula as with lm
or lm.cluster
, but also specify the clusters
parameter:^{ 5}
library(estimatr)
m3 < lm_robust(ln_wage ~ age + tenure + union + tenure:union + idcode,
clusters = idcode,
data = nlswork)
summary(m3)
Call: lm_robust(formula = ln_wage ~ age + tenure + union + tenure:union + idcode, data = nlswork, clusters = idcode) Standard error type: CR2 Coefficients: Estimate Std. Error t value Pr(>t) CI Lower CI Upper DF (Intercept) 1.882e+00 0.141424 13.310872 8.596e14 1.5930761 2.171880 28.642 age 5.631e03 0.005726 0.983442 3.342e01 0.0061215 0.017383 26.790 tenure 2.076e02 0.010089 2.057345 5.771e02 0.0007714 0.042284 14.812 union 1.746e01 0.093790 1.861817 7.735e02 0.0209918 0.370231 20.049 idcode2 6.174e01 0.019094 32.334722 2.810e07 0.6656871 0.569086 5.283 idcode3 5.625e01 0.078607 7.156034 9.687e07 0.7273116 0.397718 18.550
Unlike lm
, lm_robust
allows to specify fixed effects in a separate fixed_effects
formula parameter which, according to the documentation, should speed up computation for many types of SEs. Furthermore, this cleans up the summary output since there are no more FE coefficients:
m3fe < lm_robust(ln_wage ~ age + tenure + union + tenure:union,
clusters = idcode,
fixed_effects = ~idcode,
data = nlswork)
summary(m3fe)
Call: lm_robust(formula = ln_wage ~ age + tenure + union + tenure:union, data = nlswork, clusters = idcode, fixed_effects = ~idcode) Standard error type: CR2 Coefficients: Estimate Std. Error t value Pr(>t) CI Lower CI Upper DF age 0.005631 0.005726 0.9834 0.33419 0.0061215 0.01738 26.790 tenure 0.020756 0.010089 2.0573 0.05771 0.0007714 0.04228 14.812 union 0.174619 0.093790 1.8618 0.07735 0.0209918 0.37023 20.049 tenure:union 0.014974 0.009043 1.6558 0.14062 0.0062982 0.03625 7.187 Multiple Rsquared: 0.7554 , Adjusted Rsquared: 0.6861 [...]
When we compare the results from lm_robust
with lm
, we can see that the point estimates are the same. The lm_robust
SEs are, as expected, higher than the “classic” SEs from lm
. However, the lm_robust
SEs are also a bit smaller than those calculated from sandwich::vcovCL
:
m3fe_df < tidy(m3fe) %>% rename(est.lm_robust = estimate,
se.lm_robust = std.error)
m3fe_df$se.sandwich < m1coeffs_cl[coi_indices,2][2:5]
m3fe_df$est.classic < m1coeffs_std[coi_indices,1][2:5]
m3fe_df$se.classic < m1coeffs_std[coi_indices,2][2:5]
m3fe_df[c('term', 'est.classic', 'est.lm_robust', 'se.classic',
'se.lm_robust', 'se.sandwich')]
term est.classic est.lm_robust se.classic se.lm_robust se.sandwich age 0.005630809 0.005630809 0.003109803 0.005725617 0.006339777 tenure 0.020756426 0.020756426 0.006964417 0.010088940 0.011149190 union 0.174619394 0.174619394 0.060646038 0.093789776 0.101970509 tenure:union 0.014974113 0.014974113 0.009548509 0.009043429 0.009646023
This is because lm_robust
by default uses a different clusterrobust variance estimator “to correct hypotheses tests for small samples and work with commonly specified fixed effects and weights” as explained in the Getting started vignette. Details can be found in the Mathematical notes for estimatr.
As with the lm
and lm.cluster
results, we can also estimate marginal effects with a lm_robust
result object. However, this doesn’t seem to work when you specify FEs via fixed_effects
parameter as done for m3fe
:
margins(m3fe, variables = 'tenure', at = list(union = 0:1)) %>%
summary()
Error in predict.lm_robust(model, newdata = data, type = type, se.fit = TRUE, :
Can't set `se.fit` == TRUE with `fixed_effects`
With m3
(where FEs were directly specified in the model formula), marginal effects estimation works and we don’t even need to pass a separate vcov
matrix, since this information already comes with the lm_robust
result object m3
.^{ 6}
margins(m3, variables = 'tenure', at = list(union = 0:1)) %>% summary()
factor union AME SE z p lower upper tenure 0.0000 0.0208 0.0101 2.0573 0.0397 0.0010 0.0405 tenure 1.0000 0.0357 0.0077 4.6174 0.0000 0.0206 0.0509
As already said, lm_robust
uses a different variance estimator than sandwich’s vcovCL
and Stata. However, by setting se_type
to 'stata'
we can replicate these “Stata Clustered SEs”:
m3stata < lm_robust(ln_wage ~ age + tenure + union + tenure:union + idcode,
clusters = idcode,
se_type = 'stata',
data = nlswork)
m3stata_se < tidy(m3stata) %>% pull(std.error)
all(near(m3stata_se, m1clse)) # same SEs?
[1] TRUE
In summary, lm_robust
is as convenient to use as lm
or lm.cluster
, but offers similar flexibility as sandwich
for estimating clustered SEs. A big advantage is that you don’t need to care about supplying the right covariance matrix to further postestimation functions like margins
. The proper covariance matrix is directy attached to the fitted lm_robust
object (and can be accessed via model$vcov
or vcov(model)
if you need to). Is parameter estimation also faster with lm_robust
?
Performance comparison
We’ll make a rather superficial performance comparison only using the nlswork
dataset and microbenchmark. We will compare the following implementations for estimating model coefficients and clustered SEs:
lm
andvcovCL
from sandwichlm.cluster
lm_robust
with default SEs (se_type = 'CR2'
)lm_robust
with Stata SEs (se_type = 'stata'
)lm_robust
withfixed_effects
parameter and Stata SEs (fixed_effects = idcode, se_type = 'stata'
)
For a fair comparison, we don’t calculate CIs (which lm_robust
by default does). These are the results for 100 test runs:
As expected, lm/sandwich
and lm.cluster
have similar run times. lm_robust
is faster for all three configurations (3. to 5.) and is especially fast when estimating Stata SEs (4. and 5.). With our example data, specifying fixed_effects
(5.) doesn’t seem to speed up the calculations.
Conclusion
We’ve seen that it’s important to account for clusters in data when estimating model parameters, since ignoring this fact will likely result in overestimated precision which in turn can lead to wrong inference. R provides many ways to estimate clustered SEs. The packages sandwich
and lmtest
come with a rich set of tools for this task (and also for other types of robust SEs) and work with lm
and other kinds of models. lm.cluster
from the miceadds package provides a more convenient wrapper around sandwich
and lmtest
. However, users should be careful to not forget passing along the separate cluster robust covariance matrix for postestimation tasks. This is something users don’t need to care for when using lm_robust
from the estimatr package, since the covariance matrix is not separate from the fitted model object. Another advantage is that lm_robust
seems to be faster than the other options.
The source code for this article can be found on GitHub.
References
Blair, Graeme, Jasper Cooper, Alexander Coppock, and Macartan Humphreys. 2019. â€śDeclaring and Diagnosing Research Designs.â€ť American Political Science Review 113 (3): 838â€“59.
Cameron, A Colin, and Douglas L Miller. 2013. â€śA Practitionerâ€™s Guide to ClusterRobust Inference.â€ť
Roberts, Molly. 2013. â€śRobust and Clustered Standard Errors.â€ť
Sheather, Simon. 2009. A Modern Approach to Regression with R. Springer Texts in Statistics. New York: SpringerVerlag.
Zeileis, Achim, Susanne KĂ¶ll, and Nathaniel Graham. 2020. â€śVarious Versatile Variances: An ObjectOriented Implementation of Clustered Covariances in R.â€ť Journal of Statistical Software 95 (1).

It is not always necessary to use an FE model and you can very well estimate robust SEs from clustered data also without FEs. ↩

This is called sandwich estimator because of the structure of the formula: Between two slices of bread there is the meat . ↩

Plugging in in eq. 2 gives which reduces to eq. 1. because ↩

Once again a strange package name. Unlike sandwich this package name is not derived from a formula, but simply stands for additional functionality for imputation with the mice package — and just happens to include clustered SE estimation. ↩

Note that this time we have to use the “bare” unquoted variable name — not a formula, not a string. Every package has a different policy! ↩

The standard
vcov
function will return the correct (cluster robust) covariance matrix for a fittedlm_robust
model, whereas it will return the “classic” covariance matrix for a fittedlm
model. ↩
Pingback: Clustered standard errors with R  File4Me.com
https://rdrr.io/cran/rms/man/robcov.html