knitr::opts_chunk$set(echo = TRUE)

Whether analyzing a block-randomized experiment or adding fixed effects for a panel model, absorbing group means can speed up estimation time. The `fixed_effects`

argument in both `lm_robust`

and `iv_robust`

allows you to do just that, although the speed gains are greatest with "HC1" standard errors. Specifying fixed effects is really simple.

library(estimatr) lmr_out <- lm_robust(mpg ~ hp, data = mtcars, fixed_effects = ~ cyl) lmr_out lmr_out$fixed_effects

Before proceeding, three quick notes:

- Most of the speed gains occur when estimating "HC1" robust standard errors, or "stata" standard errors when there is clustering. This is because most of the speed gains come from avoiding inverting a large matrix of group dummies, but this step is still necessary for "HC2", "HC3", and "CR2" standard errors.
- While you can specify multiple sets of fixed effects, such as
`fixed_effects = ~ year + country`

, please ensure that your model is well-specified if you do so. If there are dependencies or overlapping groups across multiple sets of fixed effects, we cannot guarantee the correct degrees of freedom. - For now, weighted "CR2" estimation is not possible with fixed_effects.

In general, our speed gains will be greatest as the number of groups/fixed effects is large relative to the number of observations. Imagine we have 300 matched-pairs in an experiment.

# Load packages for comparison library(microbenchmark) library(sandwich) library(lmtest) # Create matched-pairs dataset using fabricatr set.seed(40) library(fabricatr) dat <- fabricate( blocks = add_level(N = 300), indiv = add_level(N = 2, z = sample(0:1), y = rnorm(N) + z) ) head(dat) # With HC2 microbenchmark( `base + sandwich` = { lo <- lm(y ~ z + factor(blocks), dat) coeftest(lo, vcov = vcovHC(lo, type = "HC2")) }, `lm_robust` = lm_robust(y ~ z + factor(blocks), dat), `lm_robust + fes` = lm_robust(y ~ z, data = dat, fixed_effects = ~ blocks), times = 50 )

Speed gains are *considerably* greater with HC1 standard errors. This is because we need to get the hat matrix for HC2, HC3, and CR2 standard errors, which requires inverting that large matrix of dummies we previously avoided doing. HC0, HC1, CR0, and CRstata standard errors do not require this inversion.

# With HC1 microbenchmark( `base + sandwich` = { lo <- lm(y ~ z + factor(blocks), dat) coeftest(lo, vcov = vcovHC(lo, type = "HC1")) }, `lm_robust` = lm_robust( y ~ z + factor(blocks), dat, se_type = "HC1" ), `lm_robust + fes` = lm_robust( y ~ z, data = dat, fixed_effects = ~ blocks, se_type = "HC1" ), times = 50 )

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.