Absorbing Fixed Effects with estimatr

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)

## Warning: Assigning non-quosure objects to quosure lists is deprecated as of rlang 0.3.0.
## Please coerce to a bare list beforehand with `as.list()`
## This warning is displayed once per session.

lmr_out

##       Estimate Std. Error   t value  Pr(>|t|)    CI Lower    CI Upper DF
## hp -0.02403883 0.01503818 -1.598521 0.1211523 -0.05484314 0.006765475 28

lmr_out$fixed_effects

##     cyl4     cyl6     cyl8 
## 28.65012 22.68246 20.12927

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.

Speed gains

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)
)

## Warning: `is_lang()` is deprecated as of rlang 0.2.0.
## Please use `is_call()` instead.
## This warning is displayed once per session.

## Warning: `lang_name()` is deprecated as of rlang 0.2.0.
## Please use `call_name()` instead.
## This warning is displayed once per session.

## Warning: `lang_modify()` is deprecated as of rlang 0.2.0.
## Please use `call_modify()` instead.
## This warning is displayed once per session.

head(dat)

##   blocks indiv z          y
## 1    001   001 1  1.4961828
## 2    001   002 0 -0.8595843
## 3    002   003 1  0.1709400
## 4    002   004 0 -0.3215731
## 5    003   005 1 -0.3037704
## 6    003   006 0 -1.4214866

# 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
)

## Unit: milliseconds
##             expr       min        lq      mean    median        uq
##  base + sandwich 223.29708 228.49406 230.72521 229.36891 232.98812
##        lm_robust  84.48679  85.38703  90.32777  86.87279  89.65503
##  lm_robust + fes  49.38458  49.74730  52.95657  50.20422  52.57565
##       max neval
##  246.2106    50
##  162.7008    50
##  142.9698    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
)

## Unit: milliseconds
##             expr        min         lq      mean     median         uq
##  base + sandwich 223.360276 228.206416 231.23346 229.162694 232.744454
##        lm_robust  70.638566  71.621418  74.93372  72.218917  76.005272
##  lm_robust + fes   8.516474   9.030435   9.18524   9.207212   9.339308
##        max neval
##  297.74939    50
##  146.60404    50
##   10.06247    50