estimatr 2.0: what changes, what does not, and why

The short version

estimatr 2.0 is a ground-up rewrite, with the main goals of speed and improved code robustness. Almost everything from 1.x runs unchanged and returns identical estimates, identical standard errors, identical confidence intervals, and identical objects. That claim is underpinned by 5,635 test assertions, 695 of which compare against answers recorded from an installed estimatr 1.0.6, coefficient by coefficient and standard error by standard error, and a further 808 of which compare against implementations that share no lineage with estimatr at all. See How this was checked. No call is slower, with gains concentrated where the algorithms changed. See Speed.

estimatr 2.0 introduces a handful of breaking changes:

horvitz_thompson() takes one probability argument instead of five.
Two stargazer helpers are gone.
the default standard error in a fixed effect model with clusters like lm_robust(Y ~ Z, fixed_effects = ~blocks, clusters = clusters, data = dat) changes from CR2 to CR0

estimatr 2.0.0 was rewritten with serious AI assistance (Claude Opus 5.x and Fable 5.x). Agents reworked the algorithms, explored areas of the estimation surface we had neglected, found many speedups, wrote thousands of tests, compared to many other regression implementations. The trade is some residual robotic prose in some of the documentation (though we tried our best to rewrite and enliven).

vignette("mathematical-notes") sets out the mathematical basis for the package: every estimator is defined in math and then validated to machine precision against that definition.

The release is checked twice over besides, against a recording of an installed estimatr 1.0.6, which shows the rewrite moved no answer, and against implementations that share no lineage with it (sandwich, clubSandwich, ivreg, Stata, fixest, plm, blkvar), which shows the answers are right. How this was checked sets out both layers.

What does not change

Five of the six estimators keep their 1.x signatures exactly; horvitz_thompson() is the subject of the next section:

function	signature
`lm_robust()`	`formula, data, weights, subset, clusters, fixed_effects, se_type, ci, alpha, return_vcov, try_cholesky`
`lm_lin()`	`formula, covariates, data, weights, subset, clusters, se_type, ci, alpha, return_vcov, try_cholesky`
`iv_robust()`	`formula, data, weights, subset, clusters, fixed_effects, se_type, ci, alpha, diagnostics, return_vcov, try_cholesky`
`difference_in_means()`	`formula, data, blocks, clusters, weights, subset, se_type, condition1, condition2, ci, alpha`
`lh_robust()`	`..., data, linear_hypothesis`

Every S3 method 1.x provided is present and behaves the same way: tidy(), glance(), summary(), print(), predict(), coef(), confint(), vcov(), nobs(), update(), plus the texreg and emmeans hooks. One return class changes: tidy(), glance(), and the new augment() return tibbles, as broom’s methods do, where 1.x returned plain data frames. $, [[, and row indexing are unaffected; tidy(fit)[, "estimate"] now gives a one-column tibble rather than a vector, so write tidy(fit)$estimate.

What changes, and why

1. `horvitz_thompson()` takes one probability argument

1.x asks you to say what you know about the randomization in five different ways, and the right combination depends on the design:

# estimatr 1.x. Not run: decl is a randomizr declaration, pr a vector of
# per-unit probabilities, and pr_mat a matrix of joint probabilities.
horvitz_thompson(y ~ z, data = dat, ra_declaration = decl)
horvitz_thompson(y ~ z, data = dat, condition_prs = 0.4, simple = TRUE)
horvitz_thompson(y ~ z, data = dat, blocks = bl, condition_prs = pr)
horvitz_thompson(y ~ z, data = dat, clusters = cl, condition_prs = pr)
horvitz_thompson(y ~ z, data = dat, condition_pr_mat = pr_mat)

2.0 has one argument, condition_prs, which takes whichever object you actually have:

# estimatr 2.0. Not run; the same objects as above.
horvitz_thompson(y ~ z, data = dat, condition_prs = decl)                     # ra_declaration
horvitz_thompson(y ~ z, data = dat, condition_prs = c("0" = 0.6, "1" = 0.4))  # named vector
horvitz_thompson(y ~ z, data = dat, condition_prs = pr_mat)                   # per-unit matrix

The blocks, clusters, simple, ra_declaration, condition_pr_mat, subset and return_condition_pr_mat arguments are gone, and se_type = "constant" is gone.

The rationale. An ra_declaration already carries the block structure, the cluster structure, the per-unit marginal probabilities, and whether the randomization was simple or complete. Passing blocks and clusters separately restates information the declaration already holds, and gives the estimator two sources of truth that can disagree. Any parametric design you can describe in words is a declare_ra() call:

declare_ra(blocks = bl, clusters = cl, prob = pi, simple = FALSE)

and any design you cannot describe parametrically is a permutation matrix:

declare_ra(permutation_matrix = perm)

which replaces condition_pr_mat and, in the process, replaces 538 lines of matrix-construction helper code.

The upside, part one: multi-arm designs. 1.x refuses an ra_declaration with more than two arms. 2.0 contrasts any two arms of one, with condition1 and condition2 picking the contrast. The estimand stays the average treatment effect over all N units of the design, so data must carry one row per unit, including the arms outside the contrast; a declaration whose size does not match nrow(data) is an error rather than a silent misalignment.

library(randomizr)
set.seed(3)
decl3 <- declare_ra(N = 300, conditions = c("control", "T1", "T2"))
Z3 <- conduct_ra(decl3)
dat3 <- data.frame(Y = rnorm(300) + 0.4 * (Z3 == "T1") + 0.8 * (Z3 == "T2"), Z = Z3)

horvitz_thompson(Y ~ Z, data = dat3, condition_prs = decl3,
                 condition1 = "control", condition2 = "T2")
#> Horvitz-Thompson estimator
#>     Estimate Std. Error  t value     Pr(>|t|)  CI Lower CI Upper DF
#> T2 0.8289969   0.145345 5.703648 1.172701e-08 0.5441258 1.113868 NA

The upside, part two: speed. Passing the declaration is what makes the O(1) variance possible. When 2.0 knows the design is complete randomization within blocks, it computes the Young’s inequality bound from six scalar sums per block rather than from an N-by-N matrix of joint probabilities. Passing a bare probability vector instead still works, and gives you the conservative simple-randomization bound, valid for any design but exact only for Bernoulli assignment. The two paths are now visibly different at the call site: an ra_declaration buys you the design-aware variance, and a plain vector buys you the conservative one. In 1.x that distinction was buried in which combination of arguments you happened to supply.

set.seed(2)
decl <- declare_ra(blocks = rep(c("a", "b", "c", "d"), each = 50), prob = 0.4)
Z <- conduct_ra(decl)
dat_ht <- data.frame(Y = rnorm(200) + 0.5 * Z, Z = Z)

horvitz_thompson(Y ~ Z, data = dat_ht, condition_prs = decl)
#> Horvitz-Thompson estimator
#>    Estimate Std. Error  t value    Pr(>|t|)  CI Lower  CI Upper DF
#> 1 0.4071871  0.1414074 2.879532 0.003982658 0.1300337 0.6843405 NA

2. Clustered fixed effects get a different default standard error

One default moves, and it is the only one. fixed_effects used together with clusters now defaults to se_type = "CR0", where 1.x defaulted to "CR2". Every other fixed_effects call gives you exactly what 1.x gave you, to the last bit, only much faster. Point estimates are unchanged everywhere, including in the clustered case.

set.seed(343)
dat <- data.frame(
  y = rnorm(1000), x = rnorm(1000), z = rbinom(1000, 1, 0.5),
  cl = rep(1:100, each = 10), bl = rep(1:50, each = 20), bl2 = rep(1:20, times = 50)
)

# unclustered, one factor: HC2, as in 1.x
lm_robust(y ~ z + x, data = dat, fixed_effects = ~ bl)$se_type
#> [1] "HC2"

# unclustered, two factors: still HC2, as in 1.x
lm_robust(y ~ z + x, data = dat, fixed_effects = ~ bl + bl2)$se_type
#> [1] "HC2"

The change warns once per session rather than once per call, since absorbed fixed effects are usually fitted in a loop. Writing se_type = "CR0" accepts the new default and removes the warning; writing se_type = "CR2" gets the 1.x number back exactly, also without warning.

Why everything else got faster. HC2 and HC3 are built from the leverage values of the full design matrix, the one with every fixed-effect dummy in it. Absorbing fixed effects by demeaning is precisely the decision not to build that matrix, which is why 1.x had to expand the dummies to compute them, and why it took 41 seconds to fit 40,000 observations across 2,000 blocks. The way out is an identity. The projection onto the full design splits exactly into the projection onto the dummies and the projection onto the demeaned covariates:

P_[X | D] = P_D + P_{M_D X}

so each leverage value is the demeaned-X hat value plus diag(P_D), and that second piece is cheap: with one factor it is just each unit’s weight share within its own group, and with several it costs a factorisation the size of the design’s narrowest dimension. No dummy matrix is ever built, at any number of factors. Same numbers, without the matrix that made them slow.

Why CR2 is the exception. Its adjustment comes from cluster-level blocks of the hat matrix rather than from the leverage diagonal, and blocks do not decompose the way the diagonal does. So CR2 still has to write the dummies out, which is roughly cubic in the number of levels and gives back the whole reason to absorb fixed effects in the first place. That is why it is the one default that moves. Ask for it by name and you still get it, exact and equal to 1.0.6:

lm_robust(y ~ z + x, data = dat, fixed_effects = ~ bl,
          clusters = cl, se_type = "CR2")$std.error
#>          z          x 
#> 0.06167608 0.03165573

What the speed buys. Absorbing a large number of fixed effects and still getting HC2 was impractical before. 40,000 observations across 2,000 blocks takes 41 seconds in 1.x and about 4 milliseconds in 2.0, for a bit-identical number. Two-way tells most: 50,000 observations across 1,000 x 30 groups takes 1.x 12.5 seconds and 2.0 7 milliseconds. The memory is the more telling half. Peak resident set size for the whole R process is 1,564 MB in 1.x 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 that made wide fixed effects impractical is never allocated, and because it is built in C++ rather than in R, gc() never saw it and only the process’s resident size does.

A rank-deficiency bug goes with it. If one FE factor is spanned by the others (a nested factor, or a disconnected design), the FE design is rank deficient. 1.x expanded the dummies, let a pivoted QR drop the redundant columns, and read the hat values off the padded design, so its absorbed answer disagreed with its own explicit-dummy fit. Taking diag(P_D) through a pseudo-inverse is exact at any rank, so 2.0 returns the dummy-regression answer, and uses the design’s true rank for the residual degrees of freedom rather than the nominal level count.

The two-way example just above is one of these designs, which is worth knowing because it was not constructed to be. Its 1,029 fixed-effect columns have rank 1,020, so nine levels are spanned by the others. 2.0’s absorbed standard error matches its own explicit-dummy fit to 4.3e-17; 1.0.6’s disagrees with its own by 1.9e-6, in the fourth significant digit. That is the only row in this vignette where the two versions do not agree to floating point, and the disagreement is 1.x being wrong rather than a change of convention.

3. Functions that are gone

starprep() and commarobust() were stargazer conveniences. Both are removed, and both remain as names that error and name the replacement, so a 1.x script says what happened rather than failing with “could not find function”.

starprep(lm(y ~ z, data = dat))
#> Error:
#> ! `starprep()` was removed in estimatr 2.0.
#> It prepared fits for stargazer, which is no longer maintained.
#> Use modelsummary, which reads `tidy()` and `glance()` and so works on every estimator in this package:
#>   modelsummary::modelsummary(list(fit1, fit2))

declaration_to_condition_pr_mat(), gen_pr_matrix_cluster() and permutations_to_condition_pr_mat() are not exported. They built the condition_pr_mat that horvitz_thompson() no longer accepts, and the new variance does not need one. These three are simply absent rather than deprecated, because unlike the stargazer helpers they have no user-facing replacement to point at.

Porting an existing script

Four greps cover every breaking change. Whatever they do not find needs no attention.

Grep for	Replace with	Where it bites
`ra_declaration`, `condition_pr_mat`, `simple =`, or `blocks`/`clusters` inside `horvitz_thompson()`	one `condition_prs =`, holding a declaration, a named vector, or a matrix	Horvitz-Thompson designs
`fixed_effects` with no `se_type` named, and `clusters`	nothing, or name the `se_type` you want	absorbed fixed effects with clusters: the default moved, see below
`starprep`, `commarobust`, `declaration_to_condition_pr_mat`, `gen_pr_matrix_cluster`, `permutations_to_condition_pr_mat`	the first two have no replacement here; the last three are unnecessary	table output, HT internals
`fixed_effects =` followed by anything that is not a `~`	wrap the grouping variable in a formula: `fixed_effects = ~ block`	code that passed a bare column name or a vector

The first row stops with R’s own unused argument error, which names the argument you passed but not what to use instead, so it is worth grepping for rather than waiting for. The third row is the other way round: starprep() and commarobust() still exist as functions that error and name their replacement, while the three matrix builders are simply gone and fail with could not find function. The fourth still runs, and still gives the 1.x answer, but warns: enforcing the formula is the resolution of issue #304, and a warning enforces it without breaking a working script.

The second row is the one to read carefully, because it moves a number, and it is the only place in the release where a default does. With no se_type named, lm_robust(y ~ z, fixed_effects = ~ block, clusters = cl) gets CR2 from 1.x and CR0 from 2.0. The point estimate is identical; only the standard error moves. It warns once per session, so nothing changes silently and nothing spams a simulation loop. Writing se_type = "CR0" accepts the new default and removes the warning; writing se_type = "CR2" gets you the 1.x number exactly, at the cost of expanding the dummies.

Fixed effects without clusters are unaffected, at any number of factors. The default is HC2 in both versions and the numbers are bit-identical, whether you absorb one factor or five. Earlier drafts of 2.0 did move this default to HC1 for two or more factors, on the belief that the leverage identity held only for one; it holds for any number, so there is nothing to trade away and the default stays where 1.x had it. If a published fixed-effects standard error turns on which adjustment you used, name the se_type rather than taking the default.

What is new

Five things 1.x cannot do, all of which fall out of the rewrite rather than being bolted on.

Multi-arm Horvitz-Thompson. 1.x refuses an ra_declaration with more than two arms. 2.0 contrasts any two arms of one, with the estimand still defined over all N units of the design. See What changes, and why for the call.

A joint hypothesis test. lh_robust() returns a joint_hypothesis element carrying a Wald F statistic on cluster-adjusted degrees of freedom, which 1.x declines to compute.

Residuals. residuals(fit) returns them, on the scale of the data, in the original row order for clustered fits, and structural rather than first-stage for iv_robust(). In 1.x the slot was NULL for every estimator.

A warning when a regressor is dropped. Collinear terms come back as NA coefficients in both versions; only one of them tells you which terms it dropped.

Blocked designs whose blocks are not all the same shape. 1.x either errors on such a design or applies the matched-pairs estimator to every block, big ones included. 2.0 uses the estimators of Pashley and Miratrix (2021). See below, which is worth a section of its own.

Blocked designs whose blocks are not all the same shape

1.x handles two kinds of blocked design and refuses everything between them. If every block has at least two treated and two control units, each block carries its own Neyman variance. If every block is a matched pair, the variance comes from the variation across pairs. A design with both, or a block holding one treated unit and three control units, either errors or silently applies the matched-pairs estimator to every block after a warning.

Such designs are not exotic. Coarsened exact matching, full matching, and multisite trials with one or two sites per stratum all produce them.

The substantive correction is to classify blocks by how many units each arm holds, not by how large the block is. A block of eight units with one of them treated has no more estimable within-block variance than a matched pair does; the singleton arm gives you one number, and one number has no variance. 2.0 implements the three estimators of Pashley and Miratrix (2021) on that classification:

Blocks with at least two units in each arm contribute their own Neyman variance, their equation 4.
Blocks with a singleton arm contribute through the variation across such blocks, their equation 8, the “unified” estimator, which requires no two blocks to share a size. With equal sizes it reduces to the usual matched-pairs estimator, their equation 5, which is used directly because equation 8 is undefined at two equal-sized blocks.
A design with both kinds combines the parts by squared share of the sample, their section 3.3.

Here is a design with both: twenty blocks of ten units with five treated, and twelve blocks of four units with one treated.

set.seed(7)
big <- data.frame(
  bl = rep(paste0("big", 1:20), each = 10),
  z = rep(rep(0:1, each = 5), times = 20)
)
small <- data.frame(
  bl = rep(paste0("sm", 1:12), each = 4),
  z = rep(c(1, 0, 0, 0), times = 12)
)
dat_bl <- rbind(big, small)
dat_bl$y <- rnorm(nrow(dat_bl)) + 0.3 * dat_bl$z

difference_in_means(y ~ z, data = dat_bl, blocks = bl)
#> Design:  Hybrid blocked 
#>   Estimate Std. Error  t value   Pr(>|t|)    CI Lower  CI Upper       DF
#> z 0.222323  0.1175722 1.890949 0.06073167 -0.01015715 0.4548032 137.7128

1.x will not fit that design at all. It stops with

Error: If design is not pair-matched, every block must have at least two
treated and control units.

Replace the twelve four-unit blocks with twelve matched pairs and 1.x does fit it, by warning and then treating all thirty-two blocks as pairs. The point estimates agree, because the point estimator was never at issue. The standard error comes out 23% larger and the degrees of freedom are 31 rather than 168.9, because twenty blocks that could have supplied their own within-block variance are thrown into the across-block calculation instead.

pairs <- data.frame(bl = rep(paste0("pr", 1:12), each = 2), z = rep(c(1, 0), times = 12))
dat_pr <- rbind(big, pairs)
dat_pr$y <- rnorm(nrow(dat_pr)) + 0.3 * dat_pr$z

difference_in_means(y ~ z, data = dat_pr, blocks = bl)
#> Design:  Hybrid blocked 
#>    Estimate Std. Error  t value   Pr(>|t|)   CI Lower CI Upper       DF
#> z 0.3615922  0.1414606 2.556134 0.01146546 0.08233344 0.640851 168.8879

Under estimatr 1.0.6 the same call warns that some blocks hold two units while others hold more, reports design = "Matched-pair", and returns the same estimate of 0.3616 with a standard error of 0.1741 on 31 degrees of freedom, against 0.1415 on 168.9 here.

The design element reports which case applied, so you never have to infer it from the block sizes: "Blocked" when every block has two units per arm, "Matched-pair" when every block is a pair, "Small blocks" when every block has a singleton arm without all being pairs, and "Hybrid blocked" when the design mixes the two.

difference_in_means(y ~ z, data = dat_bl, blocks = bl)$design
#> [1] "Hybrid blocked"

Degrees of freedom are not treated in the paper, which stops at the variance. 2.0 combines the two components by Welch-Satterthwaite, which reduces to n - 2K for an all-big design and to K - 1 for an all-small one, matching what each literature uses on its own.

Two designs are refused, because the variance genuinely cannot be estimated rather than because the software is unwilling. The first is a design with exactly one block holding a singleton arm, which leaves nothing to compare that block against and would contribute a variance of zero:

dat_one <- rbind(big, data.frame(bl = "sm1", z = c(1, 0, 0, 0)))
dat_one$y <- rnorm(nrow(dat_one))
difference_in_means(y ~ z, data = dat_one, blocks = bl)
#> Error in `blocked_variance_pm()`:
#> ! Only one block has a single treated or control unit (block sm1).
#> The variance contributed by such blocks is estimated from the variation across them, so at least two are needed.
#> Merge that block with another, drop it, or use `lm_robust()` with block fixed effects.

The second is a set of differently-sized singleton-arm blocks in which one holds half or more of their units, which is the condition equation 8 needs to stay defined and conservative. Both messages name the offending block and suggest merging blocks or using lm_robust() with block fixed effects.

Standard errors agree to 1e-10 with blkvar::block_estimator(method = "hybrid_p"), the authors’ own implementation, across all-big, all-small, matched-pair, and hybrid designs. All-big designs and matched pairs return exactly what they always did.

Blocks of clusters are a different problem, and are not covered. Pashley and Miratrix treat treatment assigned within blocks, not blocks of clusters; clusters appear once in their paper, to be set aside. Blocked designs with clusters therefore keep the earlier estimators, and a block with a single treated or single control cluster is refused outright. See Current status for what turned up when that boundary was checked.

Speed

No call is slower in 2.0, and the gains concentrate where the algorithms changed. Ordinary lm_robust() is about 2.5x faster. Absorbed fixed effects are about 470x faster at 500 blocks and about 9,800x faster at 2,000, returning bit-identical numbers, and the two-way fit that cost 1.3 GB in 1.x costs about 27 MB here. Horvitz-Thompson under complete randomization is about 1,100x faster at N = 3,000, because the N-by-N joint inclusion probability matrix is never built.

Full tables, the measurement method, and a script that reproduces them are on the Performance page.

How this was checked

The first layer asks whether the rewrite changed an answer. 695 assertions compare 2.0 against answers recorded from an installed estimatr 1.0.6, coefficient by coefficient and standard error by standard error, across every supported standard error type, weighted and unweighted, clustered and unclustered, single and multivariate outcomes. A separate file pins the entire returned surface of sixteen fit types, names as well as values, because a package’s compatibility surface is what it returns rather than what it exports: during development six fields went missing from fitted objects and a seventh returned a wrong value, and neither a NAMESPACE diff nor a search of every reverse dependency’s source could see it, since no export changed and no call site changed either.

The second layer asks whether the answer is right, which the first cannot. Anything estimatr inherited from 1.0.6, error included, passes a comparison against 1.0.6 in silence. So 808 further assertions compare against implementations built independently of this one:

checked against	what it covers	how
`sandwich`	HC0 through HC3 and both cluster corrections, weighted and not	live, same session, 1e-10
`clubSandwich`	CR2 and its Satterthwaite degrees of freedom, including under absorption	live, same session, 1e-10
`ivreg`	2SLS HC2 and HC3, via `sandwich` on an `ivreg` fit	live, same session, 1e-10
Stata `regress`, `areg`, `ivregress`	`se_type = "stata"`, with and without absorbed fixed effects; the first-stage, endogeneity, and robust over-identification tests, via `estat`	frozen output, tolerance per value
`fixest`, `plm`	absorption, by two independent routes	recorded fixture, versions recorded
`blkvar`	the blocked-design variance, from the authors of the estimator	live
a hand-built Lin specification	`lm_lin`, including its weighted paths	live, all 36 cells

estimatr matches sandwich, clubSandwich and ivreg to machine precision everywhere they overlap, weighted included, with the CR2 Satterthwaite degrees of freedom exact.

sandwich and clubSandwich are compared live, in one session, because both sides then run on one BLAS and can be held far tighter than any recording. fixest and plm are recorded instead, with their versions, because both change small-sample defaults between releases and a live test would fail on somebody else’s release note rather than on anything here. Stata tolerances are derived per value from the digits Stata actually printed, since the printed precision spans four orders of magnitude across those tables and one constant would be either far too loose or far too tight.

Where the answers genuinely differ, the difference is asserted rather than dropped. A comparison quietly excluded because it disagrees looks exactly like one that was never written, on any green run. Weighted HC2 against Stata is therefore pinned twice, as equal to the R reference to machine precision and as different from Stata by a bounded amount, and weighted 2SLS root MSE the same way.

estimatr agrees with the maintained reference implementation, and the agreement is exact. 2SLS admits two candidate leverage values, and which one you get from sandwich depends on which 2SLS object you hand it, since sandwich has no leverage convention of its own: it calls hatvalues() on the fit. estimatr uses the leverage of the second-stage regression, h = xhat'(Xhat'Xhat)^-1 xhat, the diagonal of an orthogonal projection. So does the ivreg package, and sandwich::vcovHC() applied to an ivreg::ivreg() fit returns estimatr’s standard errors to machine precision, HC2 and HC3 alike.

The other candidate is diag(H*), where H* is the matrix taking y to the fitted values. Belsley et al. (1980) considered it for 2SLS diagnostics, observed that H* is idempotent but asymmetric, and recommended the second-stage hatvalues instead, on the ground that the diagonal of an asymmetric matrix is not a leverage. Fox et al. (2026) adopt that recommendation as the default in ivreg, whose vignette puts it plainly: the diagonal elements of H* “can’t be treated as summary measures of leverage, that is, as hatvalues.” Achim Zeileis is also the author of sandwich.

Where you will see a difference is against AER::ivreg(), whose hatvalues() method predates the ivreg package and returns diag(H*). On mtcars the two answers differ by 8.6% at HC2 and 18.5% at HC3, which is worth knowing if you are reconciling estimatr against an older script. The practical argument for the projection is the one Belsley, Kuh and Welsch gave: it lies in [0, 1] by construction, so 1 - h is never negative and HC2 is always defined. diag(H*) carries no such bound. It is already negative for one observation of mtcars, and with a weak first stage it routinely exceeds 1, at which point the correction has no square root and the standard errors are NaN.

Stata declines the question entirely: ivregress accepts unadjusted, robust, cluster, bootstrap, jackknife and hac, and refuses vce(hc2) and vce(hc3) outright. We confirmed that on Stata 17, which returns error 198 for every such request.

What none of this covers. Every number quoted here is from one machine, and a local check cannot see the cross-platform floor: the first CI run after the comparison fixture was frozen failed eleven assertions on Ubuntu and Windows and none on macOS, where the fixture had been recorded, and every difference was too small to print. Recorded comparisons therefore run at 1e-9, set from the worst case in the fixture rather than by taste. The multi-arm Horvitz-Thompson contrast has no reference implementation anywhere, since 1.x refuses the case outright, so it is ours to defend rather than something checked against a second opinion; see Current status.

The estimatr issue list

All 71 issues open on DeclareDesign/estimatr at the time of the rewrite were run against both versions, one reproducer at a time. Where they landed:

status	n
Fixed in 2.0	26
Feature requests and discussions, unchanged	23
Out of scope by design, or not actionable	7
Not reproducible without the reporter’s data	6
Superseded by the rewrite	5
Still open, real work	4

The fixes worth naming. Blocked designs whose blocks are not all the same shape (#336) now use the Pashley and Miratrix (2021) estimators, where 1.x either errors or applies the matched-pairs estimator to every block; this is the largest single addition. residuals() returns something (#345). Rank detection matches lm(), so a constant regressor comes back as NA rather than as a coefficient of 1e11 (#351, #395). predict() works with fixed effects, with factors, and with no newdata (#403, #404). lh_robust() uses cluster-adjusted degrees of freedom and returns a joint test (#405, #320, #390). Dropped collinear terms are named rather than silently returned as NA (#411). augment() exists, which opens the broom-aware packages downstream (#377). Multi-arm Horvitz-Thompson works at all (#183).

The four still open are honest holdouts. Two of them, #412 and #337, asked whether se_type = "stata" reproduces Stata’s clustered standard errors under fixed effects; the one-way case is settled and pinned in test_vs_stata.R against a restored 2019 areg fixture, which agrees to Stata’s print precision. #289 needs the reporter to say what remains after the clustered case was fixed. #123 is a running catalogue of S3 methods lm has: residuals and variable.names are now present, anova and simulate are not.

One defect turned up that has no issue behind it, and it is worth knowing if you run block-clustered designs. A block holding a single treated or single control cluster has no estimable within-block variance, and both versions returned one anyway. Exhaustive enumeration puts the estimate at 0.12 to 0.25 of the truth, depending on the block’s size, for 88% coverage of a nominal 95% interval. 2.0 refuses those designs and names the offending blocks. Matched-pair clustered designs are unaffected, since their variance comes from across blocks rather than within them.

Current status

The test suite is 5,906 assertions with none failing, R CMD check --as-cran gives 0 errors, 0 warnings and 1 NOTE for the maintainer change, and the numerical comparisons against 1.0.6 run across 50 seeds per design type for every Horvitz-Thompson design family.

What follows is the list of places where 2.0 deliberately departs from 1.x: a different number, or a refusal where 1.x answered. In each one the 1.x behaviour is wrong, and every item was reproduced by running an installed 1.0.6 rather than read off this implementation’s own code. If you have results from 1.x, these are the calls worth re-running. Nothing else in this release changes an answer.

lh_robust() used residual degrees of freedom with clustered models. With n = 100 and 10 clusters, lh_robust() tested hypotheses against 97 degrees of freedom rather than the cluster-adjusted 9, producing confidence intervals far too narrow and inconsistent with the lm_robust() fit it was built from. 2.0 looks up the per-coefficient degrees of freedom, and under CR2 gives each hypothesis its own Satterthwaite degrees of freedom, as clubSandwich::linear_contrast() does. (Issue #405.)
lh_robust() reported no joint test. 2.0 returns a joint_hypothesis element with a Wald F statistic on the same cluster-adjusted degrees of freedom. (Issue #320.)
iv_robust(diagnostics = TRUE) returned a first-stage F test with no p-value, and computed the Wu-Hausman numerator degrees of freedom from the column count of the first-stage residual matrix, which overcounts by one because the intercept residuals are collinear and get dropped. 2.0 uses the actual rank increase.
Ordered factors could not be used as clusters. class(x) %in% c("factor", "integer") returns a length-2 logical for an ordered factor, which if() rejects. (Issue #421.)
A formula stored in a variable could not be passed to fixed_effects. The quosure captured the variable name rather than the formula. (Issue #348.)
lm_robust(y ~ 1, fixed_effects = ~ block) crashed after accumulating 50 convergence warnings on a zero-column design matrix. 2.0 returns a well-formed intercept-only result with the correct residual degrees of freedom. (Issue #303.)
Residuals were not returned at all. residuals(fit) was NULL for every estimator, though fitted.values was present. 2.0 returns residuals on the scale of the data, in the original row order for clustered fits, and structural rather than first-stage residuals for iv_robust(). (Issue #345.)
```
fit <- lm_robust(y ~ z + x, data = dat)
head(round(residuals(fit), 4))
#> [1] -1.3919  0.5352  0.1360 -0.8539 -0.4197 -1.4462
all.equal(residuals(fit) + fitted(fit), dat$y, check.attributes = FALSE)
#> [1] TRUE
```
A collinear regressor was dropped silently. The dropped term comes back as an NA coefficient with no warning, which is what made an lm() user open the issue after finding coefficients that disagreed between the two functions. 2.0 warns and names the terms it dropped. (Issue #411.)
```
dat$x_copy <- dat$x
fit <- lm_robust(y ~ x + x_copy, data = dat)
#> Warning in lm_return(return_list, model_data = model_data, formula = formula):
#> Some coefficients are collinear with other regressors and were dropped, and are
#> returned as NA: x_copy.
```
A joint hypothesis test errored. lh_robust(Y ~ X1 + X2 + X3, linear_hypothesis = c("X1", "X2")) should return one joint test; 1.x replies that it “implements tests for hypotheses involving linear combinations of variables but not joint hypotheses.” The joint_hypothesis element above covers it. (Issue #390.)
An all-NA outcome errored when weights were given, though the same model without weights returned an NA coefficient. The asymmetry bites when one model is fitted across many subgroups and one subgroup has no observed outcome. (Issue #370.)
HC3 returned a silently inflated standard error, and HC2 returned NaN, on a near-saturated design. A hat value is a projection diagonal and cannot exceed 1, but rounding puts it marginally above on designs close to saturation. On Y ~ Z * factor(x) with 25 levels of x in 40 rows, 1.x gives se(Z) = 39.5 against a classical standard error of 1.42, with no warning, and roughly 39.5 on nearly every other coefficient. HC2 gives NaN throughout, because sqrt() of the resulting negative term poisons the whole variance matrix however small it is. The direction is upward, so what it costs is a finding rather than a false one, but saturated specifications and treatment-by-stratum interactions reach it easily. 2.0 drops those observations from the variance and warns with a count.
horvitz_thompson() left the unidentifiable pairs out of a custom design’s variance. With a permutation_matrix declaration, two units that can never appear together contribute a term no design of that shape identifies. 1.x dropped those terms. Enumerating all ten assignments of a 5-cluster, m = 2 design, 1.x’s mean estimated variance is 0.838 of the true sampling variance, and it returns NA on one of the ten. The direction here is downward: intervals too narrow. 2.0 bounds those terms by Young’s inequality instead (Aronow and Samii 2013), so the standard error is larger on purpose.
lm_robust() counted zero-weight rows as observations. nobs and df.residual included rows whose weight is 0, so classical standard errors came back too small by an amount that grows with the share zeroed: 2.6% at 5% of rows, 13.5% at 25%, 29.6% at 50%. HC0, HC2 and HC3 are untouched to within 0.1%, so the default was mostly safe; classical, HC1 and stata were not. Trimming an inverse-probability weight to zero is the usual way in. 2.0 counts them the way lm() does, and the rows still appear in residuals and fitted.values.
A clusters variable with one level gave CR2 standard errors of about 1.5e-17, silently. CR2’s degrees of freedom are Satterthwaite, so the J - 1 = 0 guard that stops the other cluster-robust estimators never fired. One level usually means a cluster variable collapsed upstream, which turns a data-handling mistake into an apparently infinitely precise estimate. 2.0 errors.
offset() in a formula was silently ignored, so the coefficient was the one from the model without the term: 0.904333 where the offset model gives 0.953607. 2.0 errors and names the rewrite, y - <offset> ~ x.
glance() reported the first coefficient’s degrees of freedom in a column named df.residual. Under CR2 that is Satterthwaite: 4.83 on a fit whose residual degrees of freedom are 98. No estimate moves; a reported table is wrong.
A clustered iv_robust() fit reported an overidentification test that ignored the clustering. With diagnostics = TRUE and clusters, 1.x returned the heteroskedasticity-robust score statistic, the number the same model gives unclustered. Where the instruments and the errors both vary by cluster, that statistic is too large: in 1,000 draws of 50 clusters of 10 with valid instruments, it rejected at the 5% level 36.1% of the time. The direction is toward rejection, so valid instruments were declared invalid. 2.0 sums the score’s variance within clusters, which rejects 4.0% of the time on the same draws.
A weighted over-identified iv_robust() fit returned NA for the overidentification test without saying why. 2.0 computes it: under se_type = "classical" as Sargan’s statistic on the weighted model, valid when the classical weighted standard errors are, and under a robust se_type as Wooldridge’s score test on the weighted score’s HC0 or CR0 sandwich, the form Stata computes too, valid when the robust ones are.
emmeans::emmeans() on an lm_robust fit failed unless emmeans was attached. An error rather than a wrong number.

Still open. HAC (Newey-West) and CR3 standard errors are not implemented, and neither is se_type = "constant" for Horvitz-Thompson. The multi-arm Horvitz-Thompson contrast is ours to defend rather than a port: 1.x refuses the case outright, so there is no reference implementation to check against. It is conservative in every design simulated so far, with the estimated standard error running 1.08 to 1.32 times the true sampling standard deviation, but nobody has published the bound for a two-of-K contrast.

References

Aronow, Peter M., and Cyrus Samii. 2013. “Conservative Variance Estimation for Sampling Designs with Zero Pairwise Inclusion Probabilities.” Survey Methodology 39 (1): 231–41.

Belsley, David A., Edwin Kuh, and Roy E. Welsch. 1980. Regression Diagnostics: Identifying Influential Data and Sources of Collinearity. John Wiley & Sons. https://doi.org/10.1002/0471725153.

Fox, John, Christian Kleiber, and Achim Zeileis. 2026. Ivreg: Instrumental-Variables Regression by 2SLS, 2SM, or 2SMM, with Diagnostics. https://CRAN.R-project.org/package=ivreg.

Pashley, Nicole E., and Luke W. Miratrix. 2021. “Insights on Variance Estimation for Blocked and Matched Pairs Designs.” Journal of Educational and Behavioral Statistics 46 (3): 271–96. https://doi.org/10.3102/1076998620946272.