estimatr is built for workflows that fit the same model thousands of times: simulation studies, randomization inference, and the design diagnosis DeclareDesign performs. What follows is what each estimator costs, measured against estimatr 1.0.6 so that the gains have a scale to be read against.
The page lives on the website rather than in the installed package, because every number below is specific to one machine and will drift as hardware and compilers move. Keeping it here means a stale figure can be corrected the day someone notices rather than waiting for a release.
How these were measured
Median times on an Apple M4 Pro under R 4.6.0, comparing estimatr
1.0.6 to 2.0.0. Both versions ran on the one machine, so the platform
does not limit the comparison. Everything below comes from
data-raw/benchmark_vignette.R in the package repository,
which is the source of these tables and reproduces them.
Every cell was measured in its own fresh R process, at 200 replications where a call is cheap and as few as 5 where one call takes 41 seconds. That is not fussiness. The 1.0.6 cells that build a large matrix are bound by garbage collection rather than by arithmetic, so their timings depend on how much the heap had already grown when they ran: measured in sequence after the fixed-effects rows, Horvitz-Thompson at N = 1,000 reports 9.6 ms, and measured in a fresh process it reports 53.1 ms. The second number is the one a user meets. A benchmark whose value depends on which benchmark ran before it is not a measurement, so each row gets a process of its own.
The two cells with the largest allocations, 2,000 blocks and Horvitz-Thompson at N = 3,000, move by a few percent from run to run even in fresh processes, so their ratios are given to two significant figures.
The two versions return the same answers. Fifteen of the 18 rows are bit-identical in both coefficients and standard errors, including every row of the first table, and the largest disagreement anywhere is 1.4e-17, in one Horvitz-Thompson standard error. No row is slower in 2.0, and the number that bounds the gain is 1.3x, at weighted CR2, where the two versions do nearly the same work.
Ordinary fits
| task | 1.0.6 | 2.0.0 | speedup |
|---|---|---|---|
lm_robust, HC2, n = 1,000 |
0.90 ms | 0.36 ms | 2.5x |
lm_robust, classical, n = 1,000 |
0.84 ms | 0.34 ms | 2.5x |
lm_robust, CR2, 100 clusters |
1.21 ms | 0.66 ms | 1.8x |
lm_robust, CR2, weighted |
1.00 ms | 0.75 ms | 1.3x |
lm_robust, CR2, n = 5,000, 500 clusters |
4.96 ms | 1.78 ms | 2.8x |
lm_robust, HC2, n = 100 |
0.50 ms | 0.33 ms | 1.5x |
lm_lin |
1.09 ms | 0.49 ms | 2.2x |
iv_robust |
1.20 ms | 0.72 ms | 1.7x |
difference_in_means |
0.70 ms | 0.42 ms | 1.7x |
A simulation that calls lm_robust() 10,000 times saves
about five seconds at n = 1,000. That is the boring, everywhere-at-once
gain. The interesting gains show up in the two places where the
algorithms changed.
Absorbed fixed effects
2.0 absorbs fixed effects by alternating projections in C++, the
algorithm fixest::feols uses, rather than by building the
full design matrix of dummies. Both versions default to HC2 here and
return bit-identical numbers, so the default column is a like-for-like
comparison.
| design | 1.0.6, default (HC2) | 2.0.0, default (HC2) | speedup | 1.0.6, HC1 | 2.0.0, HC1 | speedup |
|---|---|---|---|---|---|---|
| 50 blocks, n = 1,000 | 2.27 ms | 0.54 ms | 4.2x | 1.68 ms | 0.54 ms | 3.1x |
| 500 blocks, n = 10,000 | 563 ms | 1.19 ms | 470x | 14.2 ms | 1.16 ms | 12.2x |
| 2,000 blocks, n = 40,000 | 41,002 ms | 4.17 ms | 9,800x | 124 ms | 3.82 ms | 32.4x |
The default column is where the algorithms differ most, because HC2 requires the hat matrix of the full design including every dummy, and 1.x builds it. That is why 40,000 observations with 2,000 blocks takes 41 seconds in 1.0.6. The HC1 column shows the same computation with a standard error that never needed the full hat matrix, where 2.0 is still 3x to 32x faster because the demeaning itself is faster.
Memory is the more telling half. Peak resident set size for the whole
R process on a two-way design of 50,000 observations across 1,000 by 30
groups is 1,564 MB in 1.0.6 against 292 MB here, and 265 MB of that 292
is an empty R session with the package loaded, so the fit itself costs
about 27 MB where 1.x needed 1.3 GB. The dummy matrix is never
allocated. Because it is built in C++ rather than in R,
gc() never saw it and only the process’s resident size
shows the difference.
Horvitz-Thompson
1.x builds the full N-by-N joint inclusion probability matrix. 2.0 derives the same variance from scalar sufficient statistics, so the complete-randomization case is O(1) in N rather than O(N²).
| design | 1.0.6 | 2.0.0 | speedup |
|---|---|---|---|
| complete, N = 200 | 1.01 ms | 0.11 ms | 8.9x |
| complete, N = 1,000 | 53.1 ms | 0.17 ms | 308x |
| complete, N = 3,000 | 362 ms | 0.32 ms | 1,100x |
Standard errors are identical to floating-point precision for simple, complete, blocked, clustered, and blocked-and-clustered designs.
Reproducing these numbers
Because the two versions cannot share a session, each is timed in its
own. Install 1.0.6 into a separate library, then run the timing script
twice: once as it stands, which loads the installed 2.0, and once with
the first line uncommented, which puts the 1.0.6 library ahead of the
default one so that library(estimatr) loads 1.0.6
instead.
dir.create("~/R/estimatr-1.0.6", recursive = TRUE)
remotes::install_version("estimatr", version = "1.0.6", lib = "~/R/estimatr-1.0.6")
# .libPaths(c("~/R/estimatr-1.0.6", .libPaths())) # uncomment for the 1.0.6 pass
library(estimatr)
library(microbenchmark)
packageVersion("estimatr")
set.seed(42)
n <- 1000
dat <- data.frame(
y = rnorm(n), x = rnorm(n), z = rbinom(n, 1, 0.5),
cl = rep(1:100, each = 10), bl = rep(1:50, each = 20)
)
microbenchmark(
clustered = lm_robust(y ~ z + x, data = dat, clusters = cl),
absorbed = lm_robust(y ~ z + x, data = dat, fixed_effects = ~ bl),
times = 200
)