clubSandwich provides specialized methods (called sandwich estimators) for assessing uncertainty in estimated statistical models
π Sandwich estimators are statistical methods for computing standard errors, confidence intervals, and hypothesis tests that work under conditions where more common or basic methods fail (such as in regression models with correlated errors)
π― clubSandwich has an object-oriented design, allowing for it to be used with a variety of statistical models from base R and other packages, such as linear regression models, logistic and other generalized linear models, hierarchical linear models, and instrumental variables models
π‘ Want to learn more? Visit the official documentation at https://jepusto.github.io/clubSandwich/.
estimatr is a popular package for fitting linear regression models. Itβs used by many researchers who run randomized field experiments in political science, economics, and education
π§ Key Challenge:
lm_robustobjects are very similar tolmobjects from base R, but because of their key differences, they were not compatible with clubSandwich.
To demonstrate, we can load in both packages and some data:
library(estimatr)
library(clubSandwich)## Registered S3 method overwritten by 'clubSandwich':
## method from
## bread.mlm sandwichdata(Duncan, package = "carData")
Duncan$cluster <- sample(LETTERS[1:8], size = nrow(Duncan), replace = TRUE)
head(Duncan)## type income education prestige cluster
## accountant prof 62 86 82 A
## pilot prof 72 76 83 D
## architect prof 75 92 90 A
## author prof 55 90 76 A
## chemist prof 64 86 90 E
## minister prof 21 84 87 BNow, we can also make two models using lm and
lm_robust, which will be very similar, since they will be
using the same data and formula.
Weβll be using the following columns from the Duncan dataset:
| Variable | Description |
|---|---|
| π prestige | The percentage of respondents in a social survey who rated the occupation as βgoodβ or better in prestige |
| π type | The type of occupation, with the levels professional & managerial, white-collar, and blue-collar |
| π° income | Percentage of occupational incumbents in the 1950 US Census who earned $3,500 or more per year (roughly $47,000 now) |
| π education | Percentage of occupation incumbents in 1950 who were high school graduates |
Duncan_lm <- lm(
formula = prestige ~ 0 + type + income + type:income + type:education,
data=Duncan
)
summary(Duncan_lm)##
## Call:
## lm(formula = prestige ~ 0 + type + income + type:income + type:education,
## data = Duncan)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.2629 -5.5337 -0.2431 5.1065 22.5198
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## typebc -3.9505 6.7940 -0.581 0.5645
## typeprof 28.0573 12.3657 2.269 0.0294 *
## typewc -10.9937 19.4880 -0.564 0.5762
## income 0.7834 0.1307 5.992 7.12e-07 ***
## typeprof:income -0.3691 0.2039 -1.811 0.0786 .
## typewc:income -0.3603 0.2596 -1.388 0.1736
## typebc:education 0.3196 0.2798 1.142 0.2608
## typeprof:education 0.3382 0.1519 2.226 0.0323 *
## typewc:education 0.4264 0.2300 1.854 0.0719 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.647 on 36 degrees of freedom
## Multiple R-squared: 0.9771, Adjusted R-squared: 0.9713
## F-statistic: 170.3 on 9 and 36 DF, p-value: < 2.2e-16Duncan_rob <- lm_robust(
formula = prestige ~ 0 + type + income + type:income + type:education,
data=Duncan
)
summary(Duncan_rob)##
## Call:
## lm_robust(formula = prestige ~ 0 + type + income + type:income +
## type:education, data = Duncan)
##
## Standard error type: HC2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## typebc -3.9505 6.14631 -0.6428 5.245e-01 -16.4158 8.51475 36
## typeprof 28.0573 19.73527 1.4217 1.637e-01 -11.9677 68.08225 36
## typewc -10.9937 14.61098 -0.7524 4.567e-01 -40.6262 18.63870 36
## income 0.7834 0.09391 8.3425 6.212e-10 0.5930 0.97386 36
## typeprof:income -0.3691 0.35275 -1.0465 3.023e-01 -1.0846 0.34627 36
## typewc:income -0.3603 0.20507 -1.7570 8.742e-02 -0.7762 0.05559 36
## typebc:education 0.3196 0.24116 1.3253 1.934e-01 -0.1695 0.80873 36
## typeprof:education 0.3382 0.24353 1.3888 1.734e-01 -0.1557 0.83211 36
## typewc:education 0.4264 0.14410 2.9590 5.428e-03 0.1341 0.71864 36
##
## Multiple R-squared: 0.9771 , Adjusted R-squared: 0.9713
## F-statistic: 491.1 on 9 and 36 DF, p-value: < 2.2e-16β¨ Key Observation: As you can see from the above summaries, our
lm_robustmodel displays slightly more information, such as confidence intervals, but the Estimate column of both models is the same, which means the models are effectively identical for the purpose I will be using them for here.
Wald_test is one of the most central functions of
clubSandwich, but when I started working on these
changes, an lm_robust object could not be used with it.
Letβs test whether the effect of income on prestige depends on the job type:
Wald_test(Duncan_lm,
constraints = constrain_pairwise(":education", reg_ex = TRUE),
vcov = "CR2", cluster = Duncan$cluster,
tidy = TRUE)## hypothesis test Fstat df_num df_denom p_val sig
## typeprof:education - typebc:education HTZ 0.00147 1 3.44 0.971
## typewc:education - typebc:education HTZ 0.07123 1 3.01 0.807
## typewc:education - typeprof:education HTZ 0.04658 1 2.66 0.845Now, if we repeat the same test using our lm_robust
objectβ¦
Wald_test(Duncan_rob,
constraints = constrain_pairwise(":education", reg_ex = TRUE),
vcov = "CR2", cluster = Duncan$cluster,
tidy = TRUE)## hypothesis test Fstat df_num df_denom p_val sig
## typeprof:education - typebc:education HTZ 0.00147 1 3.44 0.971
## typewc:education - typebc:education HTZ 0.07123 1 3.01 0.807
## typewc:education - typeprof:education HTZ 0.04658 1 2.66 0.845π Success! The results are the same!
π§ Behind the Scenes: Since
Wald_testis one of the main methods of clubSandwich, it depends on many other functions that also did not work withlm_robustobjects, so this required a lot of work under the hood.
You may be wondering, why would you want this as opposed to just
using a regular lm object?
πͺ Enhanced Features: Among other differences,
lm_robusthas arguments clusters and fixed_effects, whichlmdoes not have. While this caused additional incompatibilities, they were well worth ironing out.
Iβve created models using the fixed_effects argument to help demonstrate the appeal:
data("ChickWeight", package = "datasets")
N <- nrow(ChickWeight)
ChickWeight$wt <- 1 + rpois(N, 3)
ChickWeight$Chick_ordered <- ChickWeight$Chick # James' suggestion 4/16
ChickWeight$Chick <- factor(ChickWeight$Chick, ordered = FALSE)
ChickWeight$Chick_int <- as.integer(ChickWeight$Chick)
lm_fit_fe <- lm(
weight ~ 0 + Time:Diet + Chick,
data = ChickWeight
)
summary(lm_fit_fe)##
## Call:
## lm(formula = weight ~ 0 + Time:Diet + Chick, data = ChickWeight)
##
## Residuals:
## Min 1Q Median 3Q Max
## -79.06 -14.78 -1.71 14.71 87.54
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Chick18 30.3094 17.9425 1.689 0.091766 .
## Chick16 9.5705 9.7155 0.985 0.325041
## Chick15 13.2906 9.1528 1.452 0.147077
## Chick13 -5.2060 7.8542 -0.663 0.507733
## Chick9 8.1274 7.8542 1.035 0.301250
## Chick20 5.3774 7.8542 0.685 0.493869
## Chick10 10.0440 7.8542 1.279 0.201530
## Chick8 25.0938 8.0790 3.106 0.001999 **
## Chick17 19.4607 7.8542 2.478 0.013536 *
## Chick19 13.7107 7.8542 1.746 0.081458 .
## Chick4 26.2940 7.8542 3.348 0.000873 ***
## Chick6 40.7107 7.8542 5.183 3.12e-07 ***
## Chick11 56.8774 7.8542 7.242 1.59e-12 ***
## Chick3 42.7940 7.8542 5.449 7.83e-08 ***
## Chick1 38.6274 7.8542 4.918 1.17e-06 ***
## Chick12 41.0440 7.8542 5.226 2.51e-07 ***
## Chick2 46.8774 7.8542 5.968 4.42e-09 ***
## Chick5 53.6274 7.8542 6.828 2.39e-11 ***
## Chick14 78.2940 7.8542 9.968 < 2e-16 ***
## Chick7 76.9607 7.8542 9.799 < 2e-16 ***
## Chick24 -27.7331 8.2201 -3.374 0.000796 ***
## Chick30 9.5169 8.2201 1.158 0.247485
## Chick22 10.2669 8.2201 1.249 0.212219
## Chick23 17.4336 8.2201 2.121 0.034402 *
## Chick27 16.4336 8.2201 1.999 0.046102 *
## Chick28 35.9336 8.2201 4.371 1.49e-05 ***
## Chick26 37.0169 8.2201 4.503 8.26e-06 ***
## Chick25 49.1003 8.2201 5.973 4.30e-09 ***
## Chick29 47.8503 8.2201 5.821 1.02e-08 ***
## Chick21 90.5169 8.2201 11.012 < 2e-16 ***
## Chick33 -14.9497 8.2201 -1.819 0.069530 .
## Chick37 -22.1997 8.2201 -2.701 0.007144 **
## Chick36 10.2170 8.2201 1.243 0.214448
## Chick31 3.8837 8.2201 0.472 0.636794
## Chick39 9.5503 8.2201 1.162 0.245832
## Chick38 17.6337 8.2201 2.145 0.032396 *
## Chick32 32.8837 8.2201 4.000 7.23e-05 ***
## Chick40 32.8837 8.2201 4.000 7.23e-05 ***
## Chick34 44.1337 8.2201 5.369 1.19e-07 ***
## Chick35 68.4670 8.2201 8.329 7.15e-16 ***
## Chick44 15.1966 8.6159 1.764 0.078350 .
## Chick45 14.1727 8.2554 1.717 0.086609 .
## Chick43 37.5894 8.2554 4.553 6.57e-06 ***
## Chick41 23.0060 8.2554 2.787 0.005516 **
## Chick47 22.5060 8.2554 2.726 0.006621 **
## Chick49 32.3394 8.2554 3.917 0.000101 ***
## Chick46 28.6727 8.2554 3.473 0.000557 ***
## Chick50 42.0894 8.2554 5.098 4.79e-07 ***
## Chick42 43.6727 8.2554 5.290 1.80e-07 ***
## Chick48 52.2560 8.2554 6.330 5.27e-10 ***
## Time:Diet1 6.6906 0.2598 25.752 < 2e-16 ***
## Time:Diet2 8.6091 0.3418 25.185 < 2e-16 ***
## Time:Diet3 11.4229 0.3418 33.417 < 2e-16 ***
## Time:Diet4 9.6559 0.3489 27.676 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25.37 on 524 degrees of freedom
## Multiple R-squared: 0.9706, Adjusted R-squared: 0.9676
## F-statistic: 320.9 on 54 and 524 DF, p-value: < 2.2e-16lm_rob_fe <- lm_robust(
weight ~ Time:Diet,
data = ChickWeight,
fixed_effects = ~ Chick
)
summary(lm_rob_fe)##
## Call:
## lm_robust(formula = weight ~ Time:Diet, data = ChickWeight, fixed_effects = ~Chick)
##
## Standard error type: HC2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## Time:Diet1 6.691 0.2964 22.57 2.649e-79 6.108 7.273 524
## Time:Diet2 8.609 0.4951 17.39 8.321e-54 7.637 9.582 524
## Time:Diet3 11.423 0.4837 23.61 1.669e-84 10.473 12.373 524
## Time:Diet4 9.656 0.3054 31.61 1.581e-123 9.056 10.256 524
##
## Multiple R-squared: 0.8843 , Adjusted R-squared: 0.8726
## Multiple R-squared (proj. model): 0.8585 , Adjusted R-squared (proj. model): 0.8442
## F-statistic (proj. model): 592.2 on 4 and 524 DF, p-value: < 2.2e-16π― Key Advantage: As you can see, while the bottom 4 rows of the Estimate column of our first model are equal to their respective places in the second, the second does not show all the dummy variables and is much better suited to this use case, which carries over to using clubSandwich more effectively since many of clubSandwichβs functions take models as input.
π Interesting Discovery: estimatr has a function called
lm_linwhich can be used to createlm_robustobjects as well, though it functions a bit differently fromlm_robust.
Once again, I have loaded in some data and created our models, but
this time one using lm_robust and one using
lm_lin, and these models are both weighted
and clustered:
data("CO2")
N <- nrow(CO2)
CO2$wt <- 1 + rpois(N, 3)
CO2$Plant_ordered <- CO2$Plant # remove rows with NA in uptake
CO2$Plant <- factor(CO2$Plant, ordered = FALSE)
CO2$Plant_int <- as.integer(CO2$Plant)
lin <- lm_lin(
uptake ~ Treatment,
covariates = ~ conc,
data = CO2,
weights = wt,
clusters = Plant
)
CO2$conc_c <- CO2$conc - lin$scaled_center[["conc"]]
rob <- lm_robust(
uptake ~ Treatment +
conc_c + Treatment:conc_c,
data = CO2,
weights = wt,
clusters = Plant
)Now we can use vcovCR on both, which gives a
cluster-robust variance-covariance matrix for each.
vcov_lin <- vcovCR(lin, CO2$Plant, type = "CR2")
vcov_rob <- vcovCR(rob, CO2$Plant, type = "CR2")
vcov_lin## (Intercept) Treatmentchilled conc
## (Intercept) 5.972414231 -5.972414231 -1.054855e-03
## Treatmentchilled -5.972414231 21.663561528 1.054855e-03
## conc -0.001054855 0.001054855 1.031405e-05
## Treatmentchilled:conc 0.001054855 0.020779160 -1.031405e-05
## Treatmentchilled:conc
## (Intercept) 1.054855e-03
## Treatmentchilled 2.077916e-02
## conc -1.031405e-05
## Treatmentchilled:conc 5.063725e-05vcov_rob## (Intercept) Treatmentchilled conc_c
## (Intercept) 5.972414231 -5.972414231 -1.054855e-03
## Treatmentchilled -5.972414231 21.663561528 1.054855e-03
## conc_c -0.001054855 0.001054855 1.031405e-05
## Treatmentchilled:conc_c 0.001054855 0.020779160 -1.031405e-05
## Treatmentchilled:conc_c
## (Intercept) 1.054855e-03
## Treatmentchilled 2.077916e-02
## conc_c -1.031405e-05
## Treatmentchilled:conc_c 5.063725e-05π Perfect Match! Once again, they are identical, but previously, this would not have worked (correctly) with either of these objects.