The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

estimatr 2.0.1

A rank-deficient fit is reported once

The fit messages that a column was dropped, naming it; tidy() and summary() then messaged again with the count, so a rank-deficient fit inside a tidy() pipeline announced itself twice for one event. stats::lm() says it once, in the printed summary. The tidy() and summary() message is gone; the named message at the fit and the count in the printed summary stay.

The message also says what the dropped column was collinear with, which the name alone does not. In four-person cells of Y ~ Z + X_woman + X_income a covariate was sometimes constant, aliased with the intercept and leaving the treatment coefficient meaning what it did, and sometimes equal to the treatment indicator, which makes the two inseparable; both printed the same sentence. The spanning set is a regression of the dropped column on the kept ones, reported as constant, identical to a named column, or a linear combination of named columns, for up to three dropped columns.

One rule for a full-leverage observation under HC2 and HC3

Before this release, an observation at or near leverage 1 was handled by three sets that had each been right when written and had come to disagree. The variance discarded an observation at 1 - h <= 0; the count that decided whether the fit warned, and the rule that returned NA for a coefficient such observations alone identify, both used sandwich’s tolerant h > 1 - sqrt(eps). An observation at h = 1 - 5e-9 was therefore counted and made its coefficient NA, under a message saying it had been dropped, while its term stayed in the meat inflated by about 2e8. The discard now uses the same tolerance as the count and the NA rule, so the three are one set and the message is true of it.

The warning also fired on every absorbed fixed-effects fit with a singleton group, which is the most common fixed-effects structure there is, and told the user to change estimator on a fit that was correct as it stood: a singleton’s demeaned row is zero, so it carries no share of any reported coefficient’s variance and every standard error is the one the design without it gives, to machine precision. The fit now warns if and only if some standard error is NA because the discarded observations alone identified its coefficient. Where none is, the count is reported as a message, as fixest does, and it is stored on the fit as n_leverage_near_one either way. A second warning behind the first, raised for a NaN standard error with no observation near leverage 1, diagnosed it as a leverage problem; every NaN source in the variance code other than a negative variance diagonal is now trapped, and that case already warns with the right reason, so the second warning is gone.

A standard error of numerically zero is reported as NA

A variance computed as 6e-17 is rounding error, and reported as a number it carries a t statistic of 1e16 and a p-value of 0. Three shapes reach it and the entry above catches none of them, because none of them needs an observation at leverage 1. An outcome the regressors reproduce exactly, with residual degrees of freedom left over, returns that number for every coefficient; a constant outcome returns exactly 0. A factor level confined to one arm of a treatment returns it for that arm’s contrast alone, with three stars on an estimate of 2e-17, while the rest of the fit is ordinary: on a covariate balance test over 88 studies this printed a significant imbalance on levels holding one or two observations. A cluster-robust fit reaches the same number by a third route, where the within-cluster sums of scores cancel although no residual is near zero: a paired conjoint whose outcome sums to 1 within respondent returns it for every respondent-level covariate, with residuals running to plus and minus 0.5.

A coefficient’s standard error is now NA when its variance is below .Machine$double.eps times the classical variance it would carry if the regressors explained nothing, tss/df times XtX_inv(j,j). The ratio is dimensionless, so it is comparable across columns of unlike scale, and it is about 1 - R^2 for a classical fit, so a fit has to reproduce its outcome to the last bit to be caught: a design with an R-squared of 1 - 1e-9 is untouched. The rule applies under every se_type. Where the outcome is reproduced exactly, a classical fit returns the same 6e-17, so the rule is not about robust variances; where the fit is clustered, the residuals can be ordinary and it is their within-cluster sums that vanish, which se_type = "HC2" on the same data does not see. The fit messages once, saying whether the outcome is reproduced exactly, and otherwise naming the coefficients and which of the two mechanisms holds.

lm() reports these designs silently, as 1.0.6 did.

The F statistic is NA where the variance of the coefficients it tests cannot be inverted

The model F is a quadratic form in the inverse of the tested coefficients’ variance block, guarded until now by chol() raising an error. A block at a reciprocal condition number of 1e-32 does not raise one: the factorization succeeds and the inverse is assembled out of cancellation. On the balance test above, glance() returned F = 389 on p = 1e-146 for a contrast with no variance at all, and F = 1275 for another. The statistic is now NA when rcond() of that block is below the double precision floor, and when any coefficient under test has the NA standard error of the entry above. The classical branch reaches the same place through an R-squared of exactly 1 and is covered by the same rule.

Reporting NA for a coefficient’s own standard error while returning a number for the joint test over the same coefficients was the inconsistency this removes.

A factor clusters no longer loses the Horvitz-Thompson standard error

horvitz_thompson() with a clustered or blocked-and-clustered declaration returned NA for the standard error, on every draw, whenever the declaration’s clusters was a factor. The estimate was right; only the variance was lost. The cluster-level aggregation is tapply(), which groups by a factor’s levels rather than by the values present, so every cluster absent from an arm summed to NA. Under two-arm randomization every control cluster is absent from the treated arm, so the failure needed nothing unusual: an ordinary factor cluster column was enough.

It was also reported in the wrong direction. The NA reached the guard that reports an unestimable variance bound, whose message tells the user their declaration is the problem and suggests replacing it with one carrying the design’s block and cluster structure, which is what they had already supplied.

1.0.6 answers these designs and never had this bug, so it is a 2.0.0 regression. The fix restores 1.0.6’s answer to fifteen digits, and a factor, character, and integer coding of the same clusters now give the same standard error.

lm_lin() says what a dropped treatment interaction costs

A rank-deficient fit already names the columns it returned as NA. In lm_lin() that is not enough, because the treatment coefficient is the effect at the covariate means only while every covariate is interacted with the treatment, and the name of a dropped column does not say that the number the reader takes for the effect has changed meaning.

The two ways a Lin design goes rank deficient are not the same event. A duplicated or constant covariate takes its own main effect down with its interaction, and the fit that remains is lm_lin() on the covariates that are left. A covariate that is constant within one treatment arm keeps its main effect and loses only its interaction: centering happens before interacting, so x_c - Z * x_c equals -mean(x) * (1 - Z), and the intercept, the treatment indicator, the centered covariate, and its interaction are an exact four-column dependency. Which of the four is dropped is then an answer rather than a naming question. On the subgroup cells of one published replication the fitted values and r.squared agree to twelve printed digits under either drop while the treatment coefficient moves from -1.56 to -2.00, so nothing in the fit’s own goodness of fit reveals it. The effect at the covariate means is not identified in that cell under any version of this package or any other.

The message now adds a sentence in the second case only, naming the treatment coefficient the drop affects and saying that it is no longer the effect at the covariate means. It stays a message rather than a warning for the reasons in the entry below: stats::lm() signals nothing here, and a warning raised inside nested grouped dplyr verbs, which is how a subgroup analysis reaches this code, crashes dplyr 1.2.1.

A dropped treatment indicator is different again, and is now a warning. The other two cases leave a treatment coefficient that means something narrower than the caller asked for; this one leaves no treatment coefficient at all, printed as an NA under a name that is still there. A treatment indicator can be collinear with the rest of a Lin design only where its arm is empty or where the covariates reproduce it exactly, so there is no reading of such a fit that recovers an effect. Reaching it takes a constant treatment column under the order-preserving drop rule, since the only column ahead of the treatment is the intercept; 2.0’s norm-ranked pivot reached it on 6 of 30 rank-deficient fits of one published replication, which is the pairing that makes the warning worth having.

Which collinear column is dropped now follows stats::lm()

A rank-deficient design has no unique fit, so an estimator must choose which columns to drop. stats::lm() makes that choice with LINPACK’s dqrdc2, which is order preserving: it walks the columns left to right and compares each column’s residual norm, after the columns already kept are projected out, against that column’s own original norm, moving a failing column to the end. The later column of a collinear set is therefore the one dropped, and the ordering the modeller wrote is respected. 2.0 detected rank with Eigen’s ColPivHouseholderQR, which pivots by largest remaining norm and reads nothing but the numbers, so it could drop a column written first and keep one written last. 2.0 had reproduced dqrdc2’s threshold criterion without its ordering.

The consequence is that 2.0 silently changed which coefficient comes back NA on every rank-deficient fit, relative to both 1.0.6 and lm(). On a published replication whose subgroup models fit ten covariates in cells of 29 to 370 observations, 30 of 56 fits are rank deficient; lm() drops the treatment column in 0 of those 30, and 2.0 drops it in 6. Over 300 random rank-deficient designs, 2.0’s drop set matched lm()’s on 57 of 300, and 2.0.1’s matches on 300 of 300 with the surviving coefficients equal to lm()’s.

lm_solver() and the HC2, HC3, and CR2 meat share the rule, so the variance is read off the same columns the coefficients were fitted at.

None of this makes an unidentified coefficient identified. Where a regressor lies in the span of the others no unique coefficient exists, the reported number still moves when the columns are reordered, and lm() reports such numbers with no indication that they are arbitrary. What changes is that the choice is now predictable and follows the order the model was written in.

try_cholesky = TRUE could return a full set of coefficients at a rank-deficient design

try_cholesky = TRUE asks for the faster arithmetic, and it is meant to select the arithmetic rather than the answer. That path tests rank by reading the diagonal of the Cholesky factor of the normalised Gram matrix, where each entry stands for one column’s residual after the earlier columns are projected out, and it falls back to the QR when an entry is small. A Gram matrix carries about half the digits of the residual its factor describes, so an exactly dependent column, whose residual share is 1e-16, comes back with a diagonal entry of order 1e-8 to 1e-6 depending on how the cancellation falls, and the test was against dqrdc2’s own 1e-7. Whether the fallback fired was therefore decided by rounding. On y ~ a1 + a2 + a3 + b1 + b2 + b3 at 200 observations, with a2 equal to a1 plus 5% noise, a3 = a1 - a2, and b1, b2, and b3 built the same way, the smallest entry lands at 2.4e-7, no fallback fired, and the fit returned seven finite coefficients at a design of rank five, where lm() and the QR path both return two NA.

The threshold is now the Cholesky’s own resolution rather than dqrdc2’s, so a design whose factor cannot separate the columns goes to the QR, which makes the rank decision at 1e-7 in full precision. No design above the threshold changes, since the fallback selects the path and not the answer. The near-dependencies that reach this code elsewhere are built as a + eps * noise, which does not cancel and so lands below 1e-7, which is why the gap survived a test that walks a design from full rank to singular in twelve steps and compares both paths against lm() at each one.

HC2 and HC3 return NA for a coefficient that a full-leverage observation alone identifies

An observation the fit reproduces exactly has leverage 1 and residual 0, so HC2’s contribution for it is a 0/0. 2.0 resolved that by setting the denominator to 0 wherever 1 - h <= 0, dropping those observations from the variance. Dropping them is free for every coefficient the observation carries no information about, and by Frisch-Waugh-Lovell a singleton dummy’s neighbours are exactly that: they keep the standard error of the design with the singleton’s row and column removed. It is not free for the coefficient those observations alone identify. On a 50-row design with a one-member factor level, 2.0 returns 0.2303 for the singleton’s coefficient against a classical 0.8501, which is 0.271 of it, and that number is assembled entirely from rows that say nothing about the coefficient it describes. Its standard error is now NA, and the other four on that design are unchanged to machine precision.

Which coefficients those are is computed rather than assumed. Observation i contributes a_ij^2 * sigma_i^2 to Var(beta_j), where a_i = (X'X)^-1 x_i, so discarding it is right exactly when a_ij is 0 and wrong otherwise. The criterion coincides with the partial leverage of Kranz (2024), h_ki = xtilde_ki^2 / sum_j xtilde_kj^2 with xtilde_k the residual from regressing x_k on the other regressors, which is positive on the same set of coefficients. The set depends on the parametrisation, as it should, since it is a claim about coefficients as written: under treatment coding one singleton level makes one coefficient non-estimable, and under sum coding the same data makes five, the intercept among them.

The result is narrower than sandwich::vcovHC(), not wider. A single NaN in the meat propagates through bread %*% meat %*% bread, so sandwich returns NaN for every coefficient on the design above, where estimatr returns four finite standard errors and one NA. Only HC2 and HC3 are affected. HC0, HC1, stata, CR0 and CR2 never form 1 - h, return a finite number here, and are left alone: what they report at full leverage is a property of those estimators rather than a computational artefact.

A dropped collinear regressor is a message, not a warning

2.0 added a warning naming the coefficients a rank-deficient fit returns as NA (#411), because a user comparing lm() and lm_robust() output had no way to see that a term had been dropped. stats::lm() signals nothing at all in that situation, at the fit, at summary(), at predict(), or through broom::tidy(): it prints Coefficients: (1 not defined because of singularities) in the summary and raises no condition. The notice is now a message, which keeps the information and matches lm() on the one axis a caller can intercept.

A caller who wrapped a rank-deficient fit in suppressWarnings() to silence the 2.0 notice will now see it printed; suppressMessages() silences it instead.

The warning also reached code no one was looking at. A warning raised inside a grouped dplyr verb that is itself nested inside another grouped verb crashes dplyr 1.2.1: cur_group_label() formats the group key with map_chr(keys, pillar::format_glimpse), one key formats to length zero, and vapply rejects it. The shape is ordinary, group_by() |> nest() |> mutate(map(data, f)) where f fits a model inside its own group_by() |> reframe(), and it is how the 2026-09-23 corpus sweep found a published archive script that ran clean under 1.0.6 and died under 2.0. A message never enters that machinery.

An underidentified iv_robust() model returns NA rather than an error

A first stage that does not reproduce every regressor stopped iv_robust(). The two neighbours it should match do not stop: lm_robust() on a rank-deficient design and AER::ivreg() on the same model each drop a column, return its coefficient as NA, and report the rest. iv_robust() now does the same. Of the two parts of the change, which column is dropped is the substantive one.

The second stage regresses the outcome on the first-stage fitted values, and lm_solver()’s pivoted QR then ranked columns by norm rather than by position. On y ~ x | z with a flat first stage it dropped the intercept and returned the intercept’s own value as the coefficient on x, which is the number a reader takes for the LATE. The drop set is now chosen before the fit: the intercept and the exogenous regressors are kept, since the first stage reproduces them exactly, and the endogenous regressors the instruments cannot reproduce come back as NA, in reverse formula order, as stats::lm() drops a collinear column. The surviving coefficients, standard errors, and degrees of freedom are the reduced model’s to 2e-15.

Formula order alone would not choose the drop set. In y ~ x + w | z + w the endogenous regressor precedes the exogenous one, and in y ~ en + x1 | x1 + dup with dup equal to x1 it is en that has no instrument. AER::ivreg() drops by position and so keeps en there, reporting a coefficient of -11.9 whose fitted values lie in the span of the intercept and x1 to 6e-16, and forming its residuals from en rather than from what was fitted. estimatr returns NA for en and the ordinary least squares fit on the intercept and x1, which is the only fit the instruments support. On mpg ~ hp + cyl | am and on y ~ x + w | z + w the two packages agree.

The message is the one lm_robust() already gives for a dropped collinear regressor, naming the coefficients returned as NA.

A rank-deficient fit reported its F statistic as NA unless the dropped column was last

Under every robust se_type, summary() on a fit with a column dropped for collinearity reported the F statistic as NA whenever that column was not the last one. On y ~ x1 + mid + en with mid collinear, lm() and estimatr’s own classical branch both returned the statistic and every robust type returned NA. The coefficient vector carries an NA where the column went, while the statistic’s indices and its variance matrix both count positions in the kept set, so the two agreed only when the gap fell at the end. The statistic is now the reduced model’s under every se_type, which is what it already was when the dropped column happened to be last. lm_robust(), lm_lin(), and iv_robust() are all affected, and iv_robust() reaches the case more often now that an underidentified regressor returns NA in place.

Absorbed fixed-effect estimates were all NA on a rank-deficient fit

lm_robust(..., fixed_effects = ~ g) returned every entry of fixed_effects as NA whenever any regressor was dropped for collinearity, although the coefficients, standard errors, and degrees of freedom were all correct. Each group effect is formed as the fitted value at a representative row minus that row’s regressors times the coefficient vector, and the coefficient vector carries an NA wherever a column went, so the product was NA for every group. A dropped column contributes nothing to the fitted values, so both sides now drop it. Unlike the F statistic above, this did not depend on where in the formula the dropped column sat.

The case is reachable whenever an absorbed factor also appears among the regressors, for example a group-level covariate alongside fixed_effects = ~ g.

A fit no longer draws from the random number stream

clean_model_data() in 1.0.6 named the hidden variables it writes into the model frame with sprintf(".__%s%%%d__", da, sample.int(.Machine$integer.max, 1)): one draw for every argument evaluated there, weights, clusters and condition_prs among them, and one more for fixed_effects. Every such call advanced the caller’s random number stream, while a fit carrying only a formula and data drew nothing. 2.0.0 numbers those names with a counter, which is as unlikely to collide with a user’s own variables and draws nothing at all.

No estimate changes. What changes is every random number a script draws after its first such fit, so a seeded simulation whose loop contains one gets a different sequence of draws under 2.0 and will not reproduce a result recorded under 1.x. Drawing one sample.int(.Machine$integer.max, 1) per such argument in a wrapper around the fit restores the old stream exactly, which is worth doing once to identify what moved a number rather than keeping in a maintained script.

The 2026-09-23 maintained-corpus sweep gives the size of it. Offer-Westort, Coppock and Green (2021) builds its randomization-inference p-values from simulations that fit a weighted or clustered model on every draw, and its footnote 15 from twenty clustered fits taken before ten million draws. Six of those p-values moved: 0.429 to 0.417, 0.199 to 0.189, 0.027 to 0.021, 0.011 to 0.015, 0.008 to 0.012, and, from a simulated null at ten million draws, 0.017851 to 0.017795. An appendix figure of simulated comparisons moved by as much as a factor of 3.6 on a coefficient. Emulating the 1.0.6 draws under 2.0.0 returns all of them to the values the reproduction repository records, which are the published ones. The repository’s other 27 output files are identical or differ in their last two digits, and across five seeds under 2.0 the first of those p-values runs from 0.417 to 0.437, so what moved is which draws the null distribution contains rather than any estimate.

lh_robust() honours ci = FALSE, and refuses se_type = "none" in its own terms

lh_robust(ci = FALSE) returned a p-value and a confidence interval anyway, and raised “no non-missing arguments to min; returning Inf” on the way. The fit with no intervals carries no degrees of freedom, and the combination took the smallest of them, so the minimum was over an empty set and the hypothesis got an infinite df. The flag is now honoured as every other estimator honours it: the estimate, its standard error, and the test statistic are returned, and the p-value, the interval, and the df are NA. 1.0.6 errored on the same call, so neither version did what was asked.

lh_robust(se_type = "none") stopped with “Object must have vcov matrix. Try setting return_vcov = TRUE in the estimator function”, which names an argument no estimator here takes and not the one the caller set. A fit with no variance gives a linear combination no standard error, and the refusal now says that.

difference_in_means() refuses a multivariate outcome, and update() refits one

difference_in_means(cbind(y1, y2) ~ z) was refused with “Must have units with both treatment conditions within each block”, which describes a different problem and sends the reader to look at a blocking variable they may not have. It now says that the estimator takes one outcome at a time.

update() on a difference_in_means fit stopped with “need an object with call component”. The fit stored no call, where horvitz_thompson() and the regression fits all store one. It does now, so update() refits as it does elsewhere.

emmeans() on a fit no longer depends on emmeans’ handling of a missing model frame

The emmeans method borrowed emmeans’ own method for lm, which reads the fit’s stored model frame. An lm_robust fit stores none, and the development emmeans 2.0.5 builds a formula from that frame’s names, so every emmeans() call on a fit errored with “Perhaps a ‘data’ or ‘params’ argument is needed”. The method now recovers the data from the call and the terms directly, as emmeans’ lm method does once the frame is set aside, and gives the same answers under emmeans 2.0.4 and both versions of 2.0.5.

estimatr 2.0.0

estimatr 2.0.0 is a ground-up rewrite aimed at the DeclareDesign use case: OLS, Lin-adjusted OLS, 2SLS IV, difference-in-means, Horvitz-Thompson, and linear hypothesis tests with heteroskedasticity- and cluster-robust standard errors. It fixes several long-standing correctness bugs, improves performance on the critical path, adds feols-style fixed effects absorption, and replaces the O(N²) Horvitz-Thompson variance with a design-aware O(1) computation. The six estimators keep their signatures, with one exception: horvitz_thompson(), whose five probability arguments consolidate into condition_prs. Removals are listed under What is dropped. Run side by side on one machine, the numbers agree to 1e-12 wherever both versions answer.

See vignette("estimatr2.0") for a user-facing tour of what changes and what does not.

Maintainer and authorship

The maintainer changes with this release, from Graeme Blair to Alexander Coppock. Graeme Blair remains an author. Macartan Humphreys, an author of the 1.x releases, has asked to be listed as a contributor instead, and his role in Authors@R moves from aut to ctb. Nothing else about the authorship changes.

How this version was written, and how to check it

estimatr 2.0.0 was written with AI. vignette("mathematical-notes") says so in full and, in the same document, gives every estimator its mathematical definition and then measures estimatr against that definition when the vignette is built, to machine precision, immediately underneath. What that guarantee is, and what it deliberately is not, is set out there rather than here.

The evidence behind the release, the counts and the two checking layers, is in that vignette and in vignette("estimatr2.0") under “How this was checked”. Wherever one of those counts appears it comes from data-raw/count_assertions.R, which measures it rather than remembering it. That script exists because the counts were four remembered numbers and drifted twice.

Where a divergence is real, it is asserted from both sides

The weighted HC2 comparison against Stata is pinned twice, as equal to the R reference to machine precision and as different from Stata by a bounded amount. The same is done for weighted 2SLS root MSE. The iv_robust HC2 leverage convention is pinned the other way round: estimatr matches sandwich::vcovHC() on an ivreg::ivreg() fit exactly, and differs from AER::ivreg(), whose hatvalues() method predates the ivreg package and returns a quantity ivreg’s own authors decline to treat as leverage. See vignette("estimatr2.0"). A comparison that is quietly excluded because it disagrees is indistinguishable, on a green run, from one that was never written.

Found by revdepcheck, and fixed here

Running the 35 CRAN reverse dependencies turned up six fields that had gone missing from fitted objects, and one that was returned with a wrong value. None of it was visible to a NAMESPACE diff, since no export changed, nor to the estimator comparisons, which name the fields they check, nor to a grep of the reverse dependencies’ sources, since the calls are unchanged. It surfaced because eventstudyr reads felevels on a fixed-effects fit and got NULL.

Restored: felevels and the unprojected tss on every fixed-effects fit; fixed_effects (the absorbed group effects) and proj_fstatistic on iv_robust() fixed-effects fits; treatment_levels on lm_lin() fits, which holds every level including the baseline and is not the same thing as the treatment_vals this version added; and nclusters, reported as 0, on an unclustered blocked difference_in_means().

proj_fstatistic was wrong, not merely missing. With absorbed fixed effects the intercept is already out of the design, and the numerator degrees of freedom subtracted one anyway, so the statistic tested a hypothesis one restriction short of the intended one. A two-regressor fixed-effects model reported F on 1 numerator degree of freedom rather than 2, with a correspondingly wrong value. A single-regressor fixed-effects model took the count to zero and dropped the statistic altogether, which is how iv_robust() came to have no proj_fstatistic at all.

felevels also came back in the wrong order. The fixed-effects matrix is coerced to character before demeaning, which re-derives the levels by string sort, so a factor with levels 1 to 30 reported them as 1, 10, 11. The names are now taken before that coercion.

felevels also had to keep listing the levels that are in the fit rather than the levels that are in the data. The first attempt at the ordering fix above read the level names off the model frame before rows were dropped for missingness, which counted levels that no observation survives into. That is not cosmetic: downstream code sizes a degrees-of-freedom correction off length(fit$felevels[[term]]). On eventstudyr’s example data, 1,700 complete rows of 2,000 leave 34 of 40 time periods in the estimation sample, and counting all 40 moved its standard errors in the fourth significant figure. felevels now carries the estimation-sample levels, in the order the data gave them.

One difference from 1.0.6 does remain, and it moves no estimate. 1.0.6 named each element after its term, ~ get(idvar) giving felevels[["get(idvar)"]], except with a single fixed-effect factor on a model fitted with missing data, where it fell back to V1. 2.0 names it consistently. Code that reads felevels$V1 will need the term name instead.

A note on the comparison tolerance

The suite compares against a fixture of answers recorded from an installed 1.0.6 rather than calling that version live, which is what lets the package carry its own name. One consequence was not obvious until CI ran: the recording is a fixed set of numbers made on one machine, so comparing it to values computed on another machine has a floor set by BLAS and LAPACK rather than by this package. Before the freeze both sides were computed in the same process, so several comparisons could reasonably ask for exact equality; they cannot any more. The first CI run after the rename failed 11 assertions on Ubuntu and Windows and none on macOS, where the fixture was recorded, and every difference was too small for waldo to render.

Recorded comparisons therefore run at REF_TOL = 1e-9, set from the worst case rather than by taste: the tightest recorded quantity is the weighted adj.r.squared at -7.1e-4, whose one ulp is a relative difference of 3.1e-13, leaving about 3,200 ulps of headroom, and nothing else in the fixture comes within 100 ulps of the tolerance. Comparisons between two things computed in the same session keep their tight tolerances, since those have no platform floor.

tests/testthat/test_return_surface.R now pins the whole returned surface of sixteen fit types against answers recorded from 1.0.6, so a field cannot go missing again without a test failing. A missing field is the failure mode worth guarding: fit$felevels returns NULL rather than erroring, so the reader gets a wrong answer instead of a stop.

Two capabilities that had been dropped, and are now restored

Both were removed on the understanding that 1.x had restricted them too. It had not: 1.0.6 answered every one of these calls, and answered them correctly. The reverse-dependency run is what surfaced the mistake.

HC2, HC3 and CR2 work with fixed_effects again. 1.0.6 built the full dummy matrix and took the hat values off the [X | FE dummies] design whenever the requested se_type needed them. HC2 and HC3 now reach the same numbers without building it, at any number of factors, through the leverage decomposition below. CR2 is the exception and still expands, because its adjustment is built from cluster-level blocks of the hat matrix rather than from h_ii. Every case is exact and equal to the explicit-dummy fit, and equal to 1.0.6, which the fixture now records. The one combination 1.0.6 also refused, weighted CR2 with fixed_effects, is still refused.

A bare grouping vector passed to fixed_effects warns instead of erroring. 1.x accepted a bare column name or an already-evaluated vector. Issue #304 asked for the formula to be enforced; a warning that names the argument and shows the expected form does that without breaking working code, so the vector is still accepted and still names the term the way 1.0.6 named it. Anything that is neither a formula nor a grouping vector is an error, as before.

One default that moves, and says so

fixed_effects with clusters defaults to "CR0"; 1.0.6 defaulted to "CR2". CR2 is the one estimator here still built from cluster-level blocks of the hat matrix rather than from h_ii, so it is the one case that still has to expand the dummies, and a default that quietly costs O(g^3) in the number of levels would undo the reason for absorbing fixed effects at all. The warning fires once per session, since absorbed fixed effects are usually fitted in a loop or a simulation and a per-call warning would be noise. Writing se_type = "CR0" accepts the default and removes the warning; writing se_type = "CR2" returns the 1.0.6 number exactly, also without warning.

Unclustered fixed_effects keeps the "HC2" default at any number of factors, which is what 1.0.6 returned. An earlier draft of 2.0 fell back to "HC1" for two or more factors and warned, on the understanding that HC2 there meant expanding the dummies. It does not; see below.

What the reverse-dependency run leaves broken

Of estimatr’s 35 CRAN reverse dependencies, 34 check clean against 2.0.0. (hbal at first would not install on the test machine at all, for a local gfortran path problem unrelated to this package; corrected, it checks clean too.)

clubSandwich, RCT and statuser were broken by the two removals above and now check clean.

One remains:

• eventstudyr has three assertions that read felevels$V1 on a one-way fixed-effects fit and want the 1.0.6 fallback name, in test-EventStudyFHS.R and test-EventStudyOLS.R. Nothing else in its suite changes and no estimate moves: its one-way branches size their degrees-of-freedom correction as 1 + rank without consulting felevels, and its comparisons against recorded Stata output all hold.

What is kept

All public functions from 1.x that are relevant to the DeclareDesign workflow are present with identical interfaces:

• lm_robust(): OLS with HC0/HC1/HC2/HC3/CR0/CR2/stata SEs
• lm_lin(): Lin (2013) covariate-adjusted estimator
• iv_robust(): two-stage least squares with the same SE menu
• difference_in_means(): Neyman variance, paired, blocked, and clustered designs
• horvitz_thompson(): inverse probability weighting with Young’s inequality variance
• lh_robust(): linear hypothesis tests via car::linearHypothesis with robust variance
• S3 methods: tidy, glance, summary, print, predict, coef, confint, vcov, nobs, update

Return objects are structurally compatible with estimatr 1.0.6: field names, classes, and method dispatch are unchanged. Cross-version identity tests confirm coefficient and standard error agreement across all supported SE types, to 1e-12 when both versions are run on the same machine and to 1e-9 against the recorded fixture the suite ships.

horvitz_thompson() fits now carry the same post-estimation methods as every other estimator here. vcov(), confint(), glance(), summary() and nobs() were missing, so the first three errored and the last two fell through to the base defaults and returned something that resembled output without being the estimate. All five now agree with estimatr 1.0.6 exactly, including confint(level = ) at a non-default level, where the interval is normal rather than t because the estimator has no degrees of freedom to spend. The glance() gap was the one worth finding: without it modelsummary() did not fail, it silently dropped the goodness-of-fit rows and printed a coefficient table that looked complete.

HC2 and HC3 with fixed effects

Absorbed fixed effects used to refuse se_type = "HC2" and "HC3", on the grounds that they need hat values from the full [dummies | X] design while absorption leaves only the demeaned ones. The hat values in fact decompose exactly, for any number of FE factors, and the decomposition is written out below.

Verified to machine precision against writing the dummies out, weighted and unweighted, balanced and unbalanced, from one factor to five, and against 1.x. At n = 40,000 with 2,000 groups, one-way HC2 takes about 4 ms here against 41.0 s in 1.0.6. At n = 50,000 across 1,000 x 30 groups, two-way HC2 takes 7 ms against 12.5 s before, and the R process peaks at 292 MB against 1,564 MB, of which 265 MB is an empty session with the package loaded.

The default with fixed effects therefore returns to "HC2", at any number of factors, which is the package default everywhere else and is what 1.x returns for the same call. An earlier draft of 2.0 fell back to "HC1" for two or more factors, on the understanding that HC2 there meant expanding the dummies. It does not.

The identity generalises to any number of factors. What fails beyond one factor is the sum of each factor’s own within-group share, which ignores the cross terms and is wrong in the third decimal place; that is what an earlier draft tested and rejected. The projection itself decomposes exactly:

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

so h_ii is the demeaned-X hat value plus diag(P_D), whatever the number of factors. Writing the FE design with its widest factor in full dummies and the rest contrast-coded makes D'WD block-diagonal in its leading block, so diag(P_D) costs a factorisation of order sum_{k>1}(g_k - 1), as small as the design allows, rather than an n-by-g dummy hat matrix. Inverting that block through its eigendecomposition rather than a solve means a disconnected or nested design, where D is rank deficient, gets the pseudo-inverse and the right projection anyway.

iv_robust() gets the same treatment: its second stage runs on fitted regressors, but those are demeaned by the same fixed effects, so the decomposition applies unchanged. It never used the shortcut before, even for one factor.

"CR2" is the exception. Its adjustment is built from cluster-level blocks of the hat matrix rather than from h_ii, and blocks do not decompose the way the diagonal does, so CR2 with fixed_effects still expands the dummies exactly as 1.0.6 did. It is available and exact; it simply costs what 1.0.6 charged for it.

What this is worth. On 50,000 observations across 1,000 x 30 groups, se_type = "HC2" with two-way fixed_effects took 12.5 seconds in an R process peaking at 1,564 MB; it now takes 7 milliseconds in one peaking at 292 MB, and 265 MB of that is an empty session with the package loaded. The dummy matrix is never built, so the memory that used to make wide fixed effects impractical is simply not allocated. It is built in C++ rather than in R, which is why gc() never reported it and only the process’s resident size does.

A rank-deficiency bug, fixed. When one FE factor is spanned by the others (a nested factor, or a disconnected design), the FE design is rank deficient, and the nominal count sum(levels) - K + 1 overstates its rank. 1.0.6 used that nominal count for the residual degrees of freedom, and for HC2 and HC3 it expanded the dummies, let a pivoted QR drop the redundant columns, and then read the hat values off the padded design. Its absorbed answer therefore disagreed with its own explicit-dummy fit, and disagreed differently depending on the se_type.

2.0 takes the exact rank from the same eigendecomposition that produces the leverage, and uses it for every se_type, so an absorbed fit and the same fit with the dummies written out now agree exactly on standard errors, degrees of freedom, p-values and intervals however degenerate the design is. Absorbing fixed effects is a computational choice and should not move the answer.

This is a deliberate departure from 1.0.6, and it only bites where the FE design is genuinely rank deficient; on a design of full rank the nominal and exact counts agree and nothing moves. Canvassed against the ecosystem on a nested-factor design and a disconnected two-way design: base lm() with explicit dummies and plm (twoways within) both report the exact rank, and estimatr now matches them. fixest::feols reports the nominal one and issues no note, so it differs from its own dummy regression on both designs; that difference is fixest’s to make, and is recorded here only so the discrepancy is not mistaken for a bug in either package.

Fixing this exposed a latent crash: r_fe was left at r + fe_rank when the meat falls back to the plain XtX_inv, so X.leftCols(r_fe) overran and aborted inside Eigen. Unreachable while HC2 was refused with fixed effects, reachable now.

What is dropped

Horvitz-Thompson auxiliary arguments. horvitz_thompson() is included, but the probability specification is consolidated into the single condition_prs argument. The blocks, clusters, simple, condition_pr_mat, subset, and return_condition_pr_mat arguments and se_type = "constant" are dropped. An ra_declaration already carries the block structure, cluster structure, per-unit probabilities, and simple-versus-complete flag, so the separate arguments restated information the declaration holds, and gave the estimator two sources of truth that could disagree. Fully custom designs go through declare_ra(permutation_matrix = perm), which replaces condition_pr_mat.

tidy(), glance(), and augment() return tibbles, as broom’s methods do; 1.x returned plain data frames. $, [[, and row indexing work as before. The constructions that change are the ones that index by position and relied on a data frame dropping to a vector or a scalar: tidy(fit)[, "estimate"] now returns a one-column tibble rather than a vector, and tidy(fit)[2, 2] a one-by-one tibble rather than a number. Write tidy(fit)$estimate and tidy(fit)$estimate[2]. The two-index form is the one that travels furthest before it fails, since abs() and / both accept a one-by-one tibble in silence and pnorm() is where it stops, with “Non-numeric argument to mathematical function” two calls away from the indexing that caused it. tibble moves from Suggests to Imports.

starprep() and commarobust(). Both are removed, but they remain as names that error with a pointer to the replacement rather than being deleted, so a script written against 1.x reports what happened instead of failing with could not find function. commarobust() recomputed robust standard errors on a fitted lm, which is what lm_robust() does. starprep() prepared fits for stargazer, which has not been maintained in years; table building now goes through modelsummary, which reads tidy() and glance() and therefore works on every estimator here with no adapter.

The condition-probability-matrix helpers declaration_to_condition_pr_mat(), gen_pr_matrix_cluster(), and permutations_to_condition_pr_mat() go with the Horvitz-Thompson arguments above: the new variance computation does not build that matrix.

extract.lm_robust() and extract.iv_robust() are kept, for , which is their only consumer. They are exported as plain functions rather than registered as S3 methods because texreg looks them up by name.

The alo_star_men dataset. data/alo_star_men.rda and its documentation are removed, so data(alo_star_men) errors. It was an example dataset for the vignettes, and those go too.

The nine 1.0.6 vignettes, including the widely cited mathematical-notes. vignette("estimatr2.0") is the one that ships. The 1.0.6 versions remain on the frozen 1.0.6 pkgdown site.

HAC (Newey-West) standard errors. Not implemented.

CR3 standard errors. Not implemented.

Errors in 1.0.6, and whether one could have reached a published result

Every item below was reproduced against an installed estimatr 1.0.6, not inferred from the rewrite’s own code. Each says which direction it moves a result, because that decides whether it costs a finding or manufactures one. Anything in this section is worth re-checking in existing work; anything not in it, including everything in the sections above about the rewrite’s internals, was never in a released version.

HC3 returns a silently inflated standard error, and HC2 returns NaN, on a near-saturated design. On Y ~ Z * factor(x) with 25 levels of x in 40 rows, 1.0.6 gives se(Z) = 39.5 against a classical standard error of 1.42, with no warning of any kind, and roughly 39.5 on nearly every other coefficient in the table. HC2 gives NaN for all of them. The direction is upward, so the cost is a finding lost rather than one invented, but the magnitude is large enough that a real effect disappears. Saturated specifications and treatment-by-stratum interactions produce this readily. See the #395 entry below for the mechanism.

Horvitz-Thompson with a permutation_matrix design is anti-conservative. Enumerating all ten assignments of a 5-cluster, m = 2 design, 1.0.6’s mean estimated variance is 0.838 of the true sampling variance, and it returns NA on one of the ten. Pairs of units that can never appear together contribute a term no design of that shape identifies, and 1.0.6 left them out rather than bounding them. The direction is downward: intervals too narrow, over-rejection. 2.0 applies the Aronow and Samii (2013) bound instead, so the standard error is larger on purpose.

Zero weights are counted as observations. nobs and df.residual counted rows whose weight is 0, so classical standard errors come back too small by an amount that grows with the share of rows zeroed:

The direction is downward. Two things limit the blast: HC0, HC2 and HC3 are untouched to within 0.1%, so the package default was mostly safe, and only classical, HC1 and stata move materially. Trimming an inverse-probability weight to zero, or using 0/1 weights to subset, is the way into it.

The same count reached further than nobs. A cluster whose every row has weight 0 was still a cluster: nclusters counted it, the CR0 and stata degrees of freedom were one too many, the stata scale factor used it, and the F statistic’s denominator degrees of freedom did. A zero-weight row left inside a cluster’s CR2 block moved the weighted CR2 standard errors by about 8e-4 relative to the fit without it, where clubSandwich’s CR2 does not move. A fixed-effect group whose every row has weight 0 was still a level, so the absorbed rank was one too many. adj.r.squared used the row count where df.residual used the observation count. difference_in_means() reported nobs as the row count under weights. 2.0 counts a zero-weight row as absent from every one of these: the fit with zero weights equals the fit on the rows with positive weight, field for field, for every estimator and every se_type (test_equivalence.R). Under two or more absorbed factors a zero-weight group’s empty sum reached the eigendecomposition as NaN; 1.0.6 returned NA standard errors there, and 2.0 answers.

A clusters variable with one level gives CR2 standard errors of about 1.5e-17, with no warning. CR2’s degrees of freedom are Satterthwaite, so the J - 1 = 0 guard that catches the other cluster-robust estimators never fires. One level is usually the symptom of a cluster variable that collapsed upstream, so the failure converts a data-handling mistake into an apparently infinitely precise estimate. 2.0 errors.

offset() in a formula is silently ignored. Not a standard error problem: the coefficient itself is the one from the model without the offset. On a probe where the offset model gives 0.953607, 1.0.6 returns 0.904333, which is lm() with the term deleted. 2.0 errors and names the rewrite (y - <offset> ~ x).

glance() reports the first coefficient’s degrees of freedom in a column named df.residual. Under CR2 that is a Satterthwaite quantity: 4.83 on a fit whose residual degrees of freedom are 98. No estimate moves; a reported table is wrong.

emmeans::emmeans() on an lm_robust fit fails unless emmeans happens to be attached, with “Perhaps a ‘data’ or ‘params’ argument is needed”. An error rather than a wrong number.

An underidentified iv_robust() model returns estimates. With fewer excluded instruments than endogenous regressors, 1.0.6 fits the second stage anyway, drops whichever regressor the instruments cannot reproduce, and reports the rest with standard errors. On mpg ~ hp + cyl | am it drops the intercept and returns coefficients on hp and cyl, under a warning that there are more regressors than instruments. When the count is sufficient but an instrument is a copy of an exogenous regressor, as in y ~ en + x1 | x1 + dup with dup equal to x1, it drops en, the endogenous regressor, and warns about nothing. The direction is unknown, which is the problem: what comes back is not a 2SLS estimate of anything, and it looks like one. 2.0 errors, comparing the rank of the first-stage fitted values with the rank of the regressors, which catches both cases and does not fire on regressors that are merely collinear among themselves. 2.0.1 returns NA for the unidentified regressor instead of erroring, choosing which regressor that is rather than leaving it to the pivot; see the 2.0.1 notes above.

A clustered iv_robust() fit reported an overidentification test that ignores the clustering. With diagnostics = TRUE and clusters, 1.0.6 returns Wooldridge’s heteroskedasticity-robust score statistic, the same number the model gives unclustered. Where the instruments and the errors both vary at the cluster level, that variance is too small and the statistic too large: in 1,000 draws of 50 clusters of 10 with valid instruments, the test rejects at the 5% level 36.1% of the time, and the median statistic is 2.7 times the correct one. The direction is toward rejection, so valid instruments are declared invalid. 2.0 sums the score’s variance within clusters, which rejects 4.0% of the time on the same draws, and returns NA with a warning when there are no more clusters than restrictions. Stata’s estat overid computes nothing after vce(cluster). A weighted over-identified fit, where 1.0.6 returned NA for this test without saying so, now gets the robust score test under any robust se_type, built on the weighted fit’s sandwich variance. On mtcars, with 2SLS of mpg on hp and gear instrumented by wt, am, and gear and weighted by drat / 5, that gives 1.702. Stata’s forceweights value of 1.845 is a different statistic rather than a competing answer: it is the frequency-weight score test, equal to the unweighted test on the data with each row repeated as often as its weight and rescaled to N rows, which is right only when a weight counts identical, independent copies of a row. Under se_type = "classical" with weights it gets Sargan’s statistic on the weighted model, which matches Stata’s forceweights value to its printed precision. With weights, every diagnostic is now the test on the model with each row multiplied by the square root of its weight. The first-stage F and Wu-Hausman tests are Wald tests under the fit’s own se_type, so each is valid under exactly the conditions the matching standard errors are: a classical test when the weights are inverse error variances, a robust test under heteroskedasticity. The robust score test uses the score’s HC0 or CR0 sandwich under every robust se_type, as Wooldridge (1995) and Stata define it; the HC1, HC2, HC3, CR2, and "stata" refinements correct a coefficient covariance and have no counterpart in a score test. Stata declines estat endogenous and estat overid after aweights unless forced, since it cannot tell precision weights from sampling weights, and computes the first-stage F regardless.

Bug fixes

Weighted R² formula: investigated, no change (affects lm_robust, lm_lin, iv_robust)

An early version of 2.0 “corrected” the weighted TSS from weights_internal² to weights_internal. That change was wrong and was reverted. estimatr 1.x’s weighted R² formula is correct, and 2.0 reproduces it exactly. The reasoning is recorded here so it is not re-litigated.

The standard weighted TSS is:

TSS = Σᵢ wᵢ (yᵢ - ȳ_w)²

where ȳ_w is the standard weighted mean (weights wᵢ). Inside lm_robust_fit, raw weights are rescaled to weights_internal = sqrt(w / μ_w), so weights_internal² = w / μ_w. Using weights_internal² in the TSS formula recovers the correct standard weighted mean as the centering point. Using weights_internal instead (with sqrt-weights) centers on ȳ_{sqrt(w)} ≠ ȳ_w, which breaks the WLS guarantee that RSS ≤ TSS and allows R² < 0.

2.0 uses weights_internal² throughout, matching estimatr 1.x’s formula exactly. For an intercept-only WLS model, R² = 0 exactly. For any WLS model, R² ∈ [0, 1].

iv_robust diagnostic tests (affects iv_robust with diagnostics = TRUE)

Two bugs in the diagnostic computations were identified and fixed:

First-stage F-test missing p-value. The returned vector from first_stage_ftest() lacked a p.value entry, causing glance() on iv_robust objects with diagnostics to silently drop the p-value or error downstream.

Wu-Hausman endogeneity test wrong numerator df. wu_hausman_reg_ftest() computed numdf as the number of columns in the first-stage residual matrix, which over-counted by one. numdf is now lm_resids$rank - lm_noresids$rank, the actual rank increase from adding the residuals.

The auxiliary regression was also being handed a first-stage residual for every regressor rather than for the endogenous ones only. An exogenous regressor appears in the instrument matrix, so its first-stage residual is zero up to rounding, and the correct count survived only because rank detection dropped that column of noise. It is not safe to rely on: both stats::lm() and a pivoted QR at tol = 1e-7 keep such a column, and the scale-invariant rank detection described below does too, which adds a spurious degree of freedom to the numerator. iv_robust() now builds the auxiliary design from the endogenous residuals alone, so the test no longer depends on how rank is determined.

Which regressors are endogenous is now decided by the column space, not by name. The fix above still took the endogenous regressors to be the ones whose names do not appear in the instrument formula. An exogenous regressor that the instruments span without naming it, as in y ~ en + x1 | v1 + v2 + v3 with v3 = x1 + inst, was therefore counted as endogenous: its first-stage residual is rounding error, the auxiliary regression carried that column, and the test returned 0.144 on 2 degrees of freedom where y ~ en + x1 | inst + inst2 + x1, the same model, returns 0.109 on 1. AER::ivreg() and estimatr 1.0.6 return 0.109 on 1 for both spellings, so this was a regression in the rewrite and never shipped. A regressor now counts as exogenous when the instruments reproduce it, its first-stage residual having norm at most 1e-7 of its own, and the weak-instrument F and the overidentification test are computed on a basis for the excluded instruments rather than on whichever columns happen to be named. All three diagnostics now depend only on the model, which test_invariance.R checks by writing the instruments three ways. Where the instrument formula names every exogenous regressor, no number moves beyond rounding. In the unnamed spelling the first-stage F now reports en alone, where AER::ivreg() also reports a row for x1 with a statistic of 2.7e32.

Multi-way fixed-effects demeaning stopped early on an outcome in small units

With two or more fixed_effects factors the demeaning iterates, and it stopped once a sweep changed nothing by more than 1e-8 * (1 + max|value|). The 1 made that test absolute whenever the values are small, and the outcome is demeaned on its own. An outcome measured in units of 1e-6 stopped after half the sweeps a unit-scale one takes. On a two-way design the coefficients were right, and the residuals came back wrong by 4.4e-4, the fitted values by 3.2e-3, and the standard errors by 7.1e-7, relative to the same fit in unit scale; in units of 1e-9 the residuals were meaningless. The test is now relative to the values’ own magnitude, with a floor for a column the fixed effects absorb entirely, and outcome scales from 1e-9 to 1e6 agree to 3e-14. On the benchmark design the number of sweeps and the time are unchanged. estimatr 1.0.6 has a milder form of the same failure, with standard errors off by 1.4e-5 at 1e-9. Found by test_invariance.R, which rescales the outcome of every estimator.

Ordered factor cluster variable crashes (#421)

Passing an ordered factor as the clusters argument caused the error 'length = 2' in coercion to 'logical(1)'. The class check class(x) %in% c("factor", "integer") returns a length-2 logical vector for ordered factors (whose class is c("ordered", "factor")), which R’s if() rejects. Fixed by replacing the class check with is.factor() and is.integer().

formula() object passed to fixed_effects crashes (#348)

When a pre-built formula object is stored in a variable and passed as fixed_effects = fe_formula, the rlang quosure captures the variable name (a symbol), not the formula value. Passing the raw quosure to stats::model.frame.default() produced the error invalid type (language). Fixed by calling eval_tidy() on the quosure before handing it to model.frame.default().

Intercept-only model with fixed effects crashes (#303)

lm_robust(y ~ 1, fixed_effects = ~block) crashed with length of 'dimnames' [2] not equal to array extent. After FE demeaning the intercept column is dropped, leaving a zero-column design matrix. The alternating-projections loop in demean_fes() ran 50 iterations on the empty matrix (the convergence check max(abs(empty_matrix)) returns -Inf, which never satisfies the criterion), accumulating 50+ warnings before a downstream dimension error.

The fix has two parts. First, demean_fes() skips the demeaning loop entirely when there are no non-intercept columns. Second, lm_robust() detects the zero-column case after demeaning and returns a well-formed result directly: empty coefficient vector, correct df.residual = n - fe_rank, full-model R², and fitted values equal to the within-group means.

lh_robust CIs inconsistent with lm_robust when clusters present (#405)

The individual hypothesis CIs and p-values from lh_robust() did not match those from lm_robust() when clusters were specified. The root cause: lh_robust() used lmr$df.residual (= n − k, e.g., 97 for n = 100, k = 3) as the degrees of freedom for t-tests and confidence intervals. For clustered models, the correct df is the cluster-adjusted value stored per-coefficient in lmr$df (e.g., 9 for a CR2 model with 10 clusters, or 8–9 via Satterthwaite approximation).

Using the residual df of 97 instead of the cluster df of ~9 produces confidence intervals that are far too narrow. The fix parses the hypothesis expression (e.g., "x=0") to extract the coefficient name ("x"), looks up lmr$df["x"], and uses that as the df. For a hypothesis involving several coefficients (e.g., "x - z=0"), CR2 gives the combination its own Satterthwaite degrees of freedom (see the CR2 entry below), and every other se_type gives all coefficients the same degrees of freedom, which the combination takes. An earlier draft of this release used the smallest per-coefficient degrees of freedom instead; that is not a conservative bound, and it was replaced before 2.0.0 shipped.

lh_robust joint hypothesis test absent (#320)

lh_robust() returned only individual hypothesis tests with no joint F-test. The fix adds a joint_hypothesis element to the return object containing a Wald F-statistic:

F = (Lβ̂)ᵀ (L Vcov Lᵀ)⁻¹ (Lβ̂) / m

where m is the number of restrictions, tested against F(m, df) using the same cluster-adjusted df as the individual tests. The vcov matrix L Vcov Lᵀ is obtained directly from car::linearHypothesis() via the vcov attribute of its return value, so the cluster-robust variance is correctly propagated.

The same work closes #390, where lh_robust(Y ~ X1 + X2 + X3, linear_hypothesis = c("X1", "X2")) should return one joint test and 1.x instead errors: “lh_robust currently implements tests for hypotheses involving linear combinations of variables but not joint hypotheses.”

Residuals are returned (#345)

Neither lm_robust() nor iv_robust() nor lm_lin() returned residuals: residuals(fit) was NULL, and fit$residuals did not exist, though fitted.values did. lm_return() dropped the field.

2.0 returns residuals on the scale of the data. Weighted fits return the unweighted residuals, so residuals + fitted.values recovers the outcome. Clustered fits are returned in the original row order, undoing the internal sort by cluster. Fixed-effects models return full-model residuals, which by Frisch-Waugh-Lovell equal the residuals of the demeaned regression. iv_robust() returns structural residuals, y − Xβ̂ on the original endogenous regressors rather than their first-stage projections. Because residuals.default() reads object$residuals, residuals(fit) now works with no new method.

Dropped collinear regressors warn (#411)

A regressor collinear with the others is dropped and returned as an NA coefficient, silently. The reporting user compared lm() and lm_robust() output, found coefficients that differed, and had no way to see that a term had been dropped and the rest were therefore conditional on a different set of regressors.

2.0 warns once, naming the dropped terms. An under-identified iv_robust() call no longer reaches this warning: it is refused, for the reason given under “Errors in 1.0.6”.

augment() for use with broom-aware packages (#377)

augment() returns the model frame with .fitted and .resid appended, which is what downstream packages expect before they will accept a model object. It was not implementable before this release, since no estimator returned residuals (#345). augment(newdata = ) predicts on new data instead. Multivariate outcomes are refused. Methods are provided for lm_robust (covering lm_lin) and iv_robust.

Closed by the rewrite rather than by a fix

• #233, Horvitz-Thompson messages. estimatr 1.x prints `simple` = FALSE, using complete cluster randomization and `condition_prs` not found, estimating probability of treatment to be constant... on every call, which in a simulation means every iteration. 2.0 requires condition_prs, so it never guesses a probability and never prints either message.
• #334, faster fixed effects. The issue asks for fixest’s algorithm. 2.0 uses it: alternating projections in C++, 4x to 10,000x faster on the designs measured above, the range widening with the number of levels.
• #365, replace AER tests with ivreg, and #402, skip tests relying on suggested packages. Both are true by construction here: the suite never used AER, and every test needing randomizr, carData or blkvar is guarded with skip_if_not_installed(). The comparisons against 1.0.6 need no guard at all, because they read recorded answers rather than calling the older version.
• #260, bad defaulting to matched pairs. The reporter has 17 blocks, one of which holds a single treated and single control unit, and estimatr 1.x collapses the whole design to matched pairs on 16 degrees of freedom. With two or more such blocks 2.0 now estimates the design properly; with exactly one, as here, it errors naming the block, because one small block leaves nothing to estimate its contribution against.

R CMD check --as-cran runs with 0 errors, 0 warnings and 1 NOTE, the maintainer change.

Blocked designs with blocks of different sizes (#336)

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 with a block holding one treated and three control units, either errors (“every block must have at least two treated and control units”) or silently applies the matched-pairs estimator to every block, big ones included, 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.

2.0 implements the estimators of Pashley and Miratrix (2021). Blocks are classified by how many units each arm holds, not by how large the block is, which is the substantive correction: a block of eight units with one of them treated has no more estimable within-block variance than a matched pair does.

• Blocks with at least two treated and two control units contribute their own Neyman variance, their equation 4.
• Blocks with a singleton treated or control unit contribute through the variation across such blocks, their equation 8, the “unified” estimator that requires no two blocks to share a size. With equal sizes this 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.

Degrees of freedom are not treated in the paper. 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.

The design element now reports "Blocked", "Matched-pair", "Small blocks", or "Hybrid blocked". 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.

Two designs are refused, because the variance genuinely cannot be estimated rather than because the software is unwilling: exactly one block with a singleton arm, which leaves nothing to compare it against and would contribute a variance of zero, and a set of different-sized such blocks in which one holds half or more of their units, which is the condition equation 8 needs to stay defined and conservative.

Unchanged: all-big designs and matched pairs return exactly what they did before.

Blocks of clusters with a singleton cluster in an arm are now refused

Pashley and Miratrix cover treatment assigned within blocks, not blocks of clusters: clusters appear once in their paper, to be set aside (“we are focusing on treatment assigned within cluster”, section 2). The hybrid above therefore does not extend to cluster-randomized blocks, and blocked designs with clusters keep the previous estimators.

Checking that boundary turned up a separate defect, present in 1.x too. A block with one treated or one control cluster has no estimable within-block variance, because the singleton arm contributes nothing to the between-cluster variability, yet both packages returned a number for it with no warning. Exhaustive enumeration of every assignment, with cluster-level potential outcomes held fixed:

With two or more clusters per arm the estimator is exact. With one it is too small by roughly the block’s cluster count, and it worsens as the block grows. End to end, blocks of 6, 6, 6, and 4 clusters where the last has a single treated cluster gave 88% coverage of a nominal 95% interval.

2.0 now errors, naming the offending blocks and suggesting either merging them so every block has two clusters per arm, or lm_robust() with block fixed effects and clusters. Matched-pair clustered designs are exempt, since every block there has a singleton arm by construction and the variance is estimated across blocks rather than within them.

Rank detection matches lm() (#351, #395)

lm_robust(y ~ x) with x constant returned coefficients of order 1e11 instead of NA. The design matrix is exactly rank deficient, but Eigen’s default ColPivHouseholderQR threshold is roughly epsilon * ncol relative to the largest pivot, tight enough that the collinear column survives as a pivot of order 1e-14. stats::lm() uses LINPACK dqrdc2 with tol = 1e-7, and 2.0 sets the same threshold.

The threshold alone is not the whole of what dqrdc2 does, and setting only the threshold introduced a second failure in the opposite direction, caught before release. Eigen compares every pivot against the largest pivot in the matrix, while dqrdc2 compares each column’s remaining norm against its own original norm. The second test is scale invariant and the first is not, so a single regressor measured in large units pushed the others below 1e-7 times the largest pivot and they were dropped as collinear on a design of full rank. Three independent normals and an intercept, with one column multiplied by 1e9, lost three coefficients to NA; the effect begins at a condition number of about 1e7, and it reached lm_robust() and iv_robust() alike, since try_cholesky defaults to FALSE. A polynomial regression discontinuity in an unscaled running variable is the ordinary case that hits it.

2.0 therefore normalizes the columns of the design to unit norm before the rank-revealing QR and undoes the scaling afterwards, at both the coefficient path and the getMeatXtX() path the hat values are read off. The criterion is then per column, which is what dqrdc2 does and what the 1e-7 threshold was meant to reproduce. Exactly collinear columns are still dropped, and a full-rank design now returns the same coefficients whatever units its columns are measured in. test_invariance.R asserts that property directly, sweeping a column through twelve orders of magnitude and requiring every coefficient and standard error to come back unchanged.

The same normalization fixed a third failure, on the try_cholesky = TRUE path, and this one is a difference from 1.0.6 rather than a pre-release catch. lm_robust(try_cholesky = TRUE) solves through a Cholesky decomposition of X'X instead of the pivoted QR, and it was returning a coefficient for every column of a rank-deficient design, splitting arbitrarily between collinear ones, where the QR path and stats::lm() return NA. It warned about the resulting NaN standard error and never about the collinearity. The guard that was supposed to prevent it tested llt.info(), which cannot catch the case: Eigen’s LLT reports success on a numerically singular Gram matrix. With the columns normalized, each diagonal entry of the Cholesky factor is the norm of that column’s residual after the earlier columns are projected out, measured against its own original norm, which is dqrdc2’s criterion again; the path now falls back to the QR below the same 1e-7. A design that is merely ill conditioned falls back too and pays the QR’s cost rather than returning a bad answer. try_cholesky is off by default, so a user reaches this only by asking for it, and estimatr 1.0.6 has the behaviour still. The two paths now agree to 3e-17, and the argument no longer changes the answer, only the speed.

A near-saturated design also produces observations that the fit reproduces exactly, whose leverage is 1. That case costs nothing in the answer: the residual is exactly 0, HC2’s contribution is a 0/0 that resolves to 0, and the standard error comes back finite, though built from fewer rows than the fit used. The failure is leverage computed marginally above 1, which rounding produces on the same designs and which no projection diagonal can actually take. 1 - h is then a small negative number, and the two estimators failed differently: HC2 divided by it and took the square root of a negative meat, so every standard error in the fit was NaN however small the offending term, while HC3 squared the denominator, which cancels the sign, and returned a finite number carrying a large spurious term without warning. The silent one was the worse failure.

Both failures are 1.0.6’s as well as the rewrite’s, confirmed by running it: on the probe design 1.0.6 returns NaN for HC2 and 39.5 for HC3 against a classical 1.42, silently. 2.0 sets the denominator to 0 wherever 1 - h <= 0, so both estimators drop those observations from the variance instead, and warns with a count of how many sat at or near leverage 1 and which se_type read it. The count is a tolerance, h > 1 - sqrt(.Machine$double.eps), which is sandwich::vcovHC()’s criterion on the same quantity, rather than a test on the sign of 1 - h. The reason is that on an exactly saturated design the standard error is the same whichever side of 1 the rounding puts the hat value, so a sign test would have left the warning to a single unit in the last place: the solver used here returns 1 + 2.2e-16 on such a design where qr() and stats::hatvalues() return exactly 1. A hat value is dimensionless and bounded by 1, so a fixed tolerance on it means the same thing in every design. CR2 never forms 1 - h and needs no guard: it eigendecomposes a per-cluster matrix and keeps only directions above 1e-12, so a cluster that is fitted exactly already contributes 0. One consequence is worth naming: on a design with a singleton fixed-effect group, the absorbed fit and the same model with explicit dummies used to disagree, the first finite and the second NaN. They now agree, in the warning as well as in the number.

predict() with fixed effects (#403, #404)

Three separate failures. predict() with no newdata errored instead of returning the in-sample fit. predict() on a fixed-effects model with a factor regressor failed with “non-conformable arguments”, because the terms object still carries the intercept that fixed-effects absorption removed from the coefficients. And 2.0 did not store the absorbed group effects at all, so no fixed-effects prediction was possible.

lm_robust() now stores the absorbed group effects for a single-outcome model with one set of fixed effects, and predict() adds them back, reproducing the dummy-variable regression exactly. New levels in newdata and multi-way fixed effects both error with an explanation: with several sets of fixed effects the absorbed effects are identified in sum but not separately, so a new observation cannot be assigned its share.

glance() on iv_robust with several endogenous regressors (#389)

The first-stage F test carries one entry per endogenous regressor, named "<var>:value" once there is more than one, so indexing it by "value" produced an NA with an NA name and the data.frame() call failed on the row name. glance() must return one row, so it now reports the weakest of the per-regressor first stages, which is the quantity a weak-instrument diagnostic asks about. The per-regressor tests remain in diagnostic_first_stage_fstatistic.

model.frame() on iv_robust (#397)

model.frame() returned a two-column frame whose second column was named "x + hp | am + hp", because model.frame.default() reads the | of a two-part formula as a logical operator. A model.frame.iv_robust() method now rebuilds a one-part formula naming every variable and defers to the default method.

Smaller fixes

• lm_robust_fit() with an unnamed or integer matrix (#269). The exported fitting function failed on an integer design matrix (“Wrong R type for mapped matrix”) and produced an object tidy() could not handle when y had no column names. Both inputs are now coerced and named.
• variable.names() (#123). Added for lm_robust and iv_robust.
• fixed_effects given a bare column name (#304). It produced “invalid formula” from deep inside model.frame. It now warns, names the argument, shows the expected form, and goes on to fit the model the way 1.x did.
• iv_robust() with multiple outcomes and diagnostics = TRUE died in the score test with “dims do not match the length of object”, in 1.0.6 as well. The diagnostic fields hold one statistic each, so the fit now warns that diagnostics are not available with multiple outcomes and skips them, as it does under fixed_effects.
• horvitz_thompson() with a cbind() outcome answered for the first outcome alone, silently: the matrix response was flattened and the first n values read as the outcome. 1.0.6 errored on the length mismatch. It now refuses and says to fit one model per outcome.
• lh_robust() with multiple outcomes (#297). Coefficients of a multivariate fit are named "<outcome>:<term>", which a hypothesis such as "cyl = 2" cannot refer to, and car::linearHypothesis() reported it as malformed. 2.0 errors with that explanation and tells the user to fit one outcome at a time.
• No observations left to fit. When missing values or subset leave no rows, 1.0.6’s lm_robust(), iv_robust() and lm_lin() returned every coefficient as NA, and an earlier 2.0 added a warning calling them collinear, which was false; difference_in_means() said both treatment conditions were needed “within each block”. All four now say that no observations are left, as lm() does.
• tidy(conf.level =) on a multivariate fit ignored the level. Intervals were matched to rows by name, and a multivariate fit’s are named "<outcome>:<term>", so nothing matched and the fit’s own 95% intervals came back under a request for 90%, silently, in 1.0.6 as well. They are now matched by position, which is how confint() builds them.
• confint() on an lh_robust fit always errored, “numbers of columns of arguments do not match”, in 1.0.6 as well, and nobs() on one returned NULL. Both now answer, and parm can name the hypothesis.
• predict() on an iv_robust() fit with fixed_effects failed with “non-conformable arguments”, where 1.0.6 returned predictions. It now drops the absorbed intercept and adds the group effects back, as predict() on lm_robust() does, and reproduces the fitted values in sample.
• lh_robust() under CR2 gave a combination of coefficients the wrong degrees of freedom. 1.0.6 refused CR2 in lh_robust(). The rewrite allowed it, and gave any hypothesis that was not a single coefficient the smallest of the coefficients’ Satterthwaite degrees of freedom, on the stated grounds that this was conservative. It is not a bound: on the test data it was 16.10 where the combination’s own is 15.75, which narrows the interval. Each hypothesis now gets its own Satterthwaite degrees of freedom, computed from the same CR2 components as the coefficients’, and matches clubSandwich::linear_contrast() unweighted, weighted, and with absorbed fixed effects. The joint test’s denominator takes the smallest of them.
• Horvitz-Thompson probabilities that cannot be probabilities. horvitz_thompson() accepted probabilities outside [0, 1] and returned an ordinary-looking estimate and standard error from them, and a probability of 0 for a condition a unit was observed in returned NaN or Inf without a word. Both are now refused, the second naming the condition and how many units it affects; the check runs after every condition_prs format has been turned into per-unit probabilities, so it covers them all.
• Infinite weights, and weights that are all zero. Both returned every coefficient as NA under a warning calling the regressors collinear, in 1.0.6 as well, where lm() refuses an infinite weight. Both are now refused with a message saying which. A missing weight is still a dropped row.
• Two refusals named the wrong problem. difference_in_means() given a condition1 and condition2 that are not values of the treatment filtered the data down to nothing and then said both conditions were needed “within each block”; it now says the conditions must be values found in the treatment. horvitz_thompson() given a named probability vector whose names miss a condition called it an “Unrecognised condition_prs format”; it now names the vector’s names and the conditions it was looking for.
• A single unit in a treatment arm of an unblocked difference_in_means(). The refusal is right, since that arm’s variance cannot be estimated, but its message said every block needed two treated and two control units, on a design with no blocks, in 1.0.6 as well. It now says what is needed without assuming blocks.

A single fixed-effects factor lost its name when rows were dropped

With one factor in fixed_effects and any observation dropped for missingness, the model frame’s one-column matrix of factor codes degrades to a vector, and rebuilding it as a matrix leaves it without a column name. Three things read that name, and all three failed. se_type = "CR2" errored outright with “‘.’ in formula and no ‘data’ argument”, because CR2 is the one estimator that still expands the fixed effects and it does so through model.matrix(~ 0 + .) on a data frame that had no names to expand. The absorbed group effects in fixed_effects came back labelled 1, 2, 3 instead of g1, g2, g3. And predict() rejected its own data as containing new levels, since the levels it was matching against had lost the term.

estimatr 1.0.6 lost the name in the same place and fell back to V1, so it returned V11, V12, V13 and its predict() failed too, but it did answer the CR2 call. The name is now taken from the level names captured before the model frame is built, which are the same ones felevels reports, so the two cannot disagree and complete and incomplete data give the same labels. The restored CR2 fit agrees with 1.0.6 and with the explicit-dummy fit to twelve digits.

Every fixed-effects test until now used complete data, which is why nothing reported this.

res_var and weights were named by accident on weighted fixed-effects fits

Both fields took whatever dimnames survived data[["y"]] - fitted.values, and R resolves that in favour of the first operand carrying any. The outcome was carrying the model frame’s row names, except on the fixed-effects path, where it has been demeaned into a fresh matrix. So a weighted fixed-effects fit disagreed with a weighted fit without fixed effects, and with 1.0.6, in two fields at once: res_var (and the proj_ R-squareds derived from it) came back named after the outcome column when 1.0.6 leaves it unnamed, and weights came back unnamed when 1.0.6 names it. Now that the row names are carried explicitly rather than inherited, all three cases agree with 1.0.6. Verified against an installed 1.0.6 and pinned in test_return_surface.R.

All-missing outcome with weights (#370)

lm_robust(X ~ 1, weights = w, data = df) where X is entirely NA errors in estimatr 1.x with “‘x’ must be an array of at least two dimensions”, while the same model without weights returns an NA coefficient. The asymmetry matters when one model is fitted across many subgroups and some subgroup has no observed outcome. 2.0 refuses both calls with the same error, “No observations are left to fit once rows with missing values, and any excluded by subset, are dropped.”, which is also what lm() does with no complete rows. An earlier draft of this entry said 2.0 returns the NA coefficient in both cases; 2.0.0 as released errors instead.

Performance improvements

CR2 degrees-of-freedom computation: O(J²) → O(r²·J)

CR2 standard errors require, for each cluster j, a scalar adjustment derived from the hat matrix block H_jj. The original 1.x implementation formed a J×J matrix P_array and computed a trace and Frobenius norm over it, which is O(J²) memory and compute when J is large.

2.0 uses an algebraic identity to reduce this to five r×r Gram matrix products, where r is the number of model parameters (typically ≪ J). For a model with 3 parameters and 500 clusters the allocation drops from a 500×500 matrix to five 3×3 matrices. Measured speedup on a 1,000-observation / 100-cluster CR2 model: approximately 2× end-to-end.

Formula parsing bypass for non-IV calls

iv_robust() requires Formula::as.Formula() to parse the two-sided formula y ~ x | z. lm_robust() and lm_lin() do not, but the original code called as.Formula() unconditionally in clean_model_data(). This path is now skipped for non-IV calls, saving approximately 80 µs per lm_robust() invocation, which is meaningful in simulations that call it thousands of times.

Horvitz-Thompson variance: O(N²) → O(1)

1.x constructs the full 2N×2N joint inclusion probability matrix for every Horvitz-Thompson call. The Young’s inequality variance for any complete-randomization design (including within blocks, and at the cluster level for clustered designs) reduces to six scalar sufficient statistics per block, so the matrix is never needed. The custom-design pathway, where an explicit permutation matrix is the only description of the design, is the one case that still forms a matrix, and does so in one tcrossprod() call in place of 1.x’s 538-line helper_condition_pr_matrix.R.

Measured against 1.x on complete randomization: 10× at N = 200, 73× at N = 1,000, 837× at N = 3,000. Standard errors are identical to floating-point precision for simple, complete, blocked, clustered, and blocked-and-clustered designs.

Multi-arm Horvitz-Thompson

estimatr 1.x refuses an ra_declaration with more than two arms. 2.0 supports the contrast between any two arms of a multi-arm design: condition1 and condition2 pick the contrast, the estimand stays the average treatment effect over all N units, and the variance uses the joint assignment probabilities implied by the arm sizes.

Two generalizations were needed. The denominator of the estimator is N, the size of the design, and not the number of units landing in the two selected conditions; those coincide only when the design has two arms. The Young’s inequality coefficients come from the pairwise joint probabilities m2(m2-1)/(n(n-1)), m1(m1-1)/(n(n-1)), and m2 m1/(n(n-1)), which do not assume m1 = n - m2. The cross coefficient collapses to 1/n for every complete design. Two-arm contrasts continue to run through 1.x’s gen_joint_pr_complete path, including its floor/ceiling mixture for non-integer implied counts, so every two-arm result is numerically unchanged.

Tested by exhaustive enumeration: over all 90 assignments of six units to three arms of two, the average estimate equals the true ATE exactly, and the closed-form variance equals the one computed from the exact joint probability matrix of the same 90 permutations.

Because probabilities are matched to units by row position, data must contain one row per unit of the declaration. Passing a declaration whose size does not match nrow(data) is now an error rather than a silent misalignment.

TSS computation vectorized

The total sum of squares calculation replaced an apply(y, 1, '-', colMeans(y)) column-mean sweep with direct centring and colSums(). 1.x ran an R-level loop over rows; an intermediate version of 2.0 used sweep(), which builds its subtrahend with array() and then transposes it with aperm(), so it materialised two extra n-by-k copies to do what recycling does in one. Speedup on a 1,000-observation multivariate model: approximately 280 µs, and two fewer allocations the size of the outcome on every fit.

Front-end overhead

At n = 100,000 the numerical core of a fit is a small fraction of the call. Most of the rest was the model frame’s row names.

stats::model.frame() turns a data frame’s compact automatic row names into a character vector as long as the data, and from there every object derived from the design matrix carries them: the fitted values from X %*% beta, the residuals formed from those, the row subset the cluster sort takes, and the drop() that turns each n-by-1 matrix into a vector. Only fitted.values keeps them in the returned object; estimatr returns residuals unnamed, as 1.0.6 did, so on that branch they were built only to be discarded.

They are now taken off the design matrix and the outcome once, in clean_model_data(), where both are still unshared and stripping them is free, and put back on fitted.values once in lm_return(). Nothing in between has to carry them. weights is named to match, as before.

The rest are local. Gaps in fixed-effect codes are detected with tabulate() rather than by hashing every observation with sort(unique()). Unweighted group sizes in the leverage computation are tabulate() rather than rowsum(), which rediscovers the groups it was already given, and the unit weight vector is no longer materialised on the one-way path that does not index it. The sums of squares behind the R-squared of a fixed-effects fit go through crossprod() instead of squaring the vector first. The absorbed group effects find a representative row per group by scatter assignment over the codes rather than by hashing all n of them. sweep() is gone from the total sum of squares, and the design matrix is no longer copied to itself when no column was dropped for collinearity.

The fixed-effects kernels

Two things dominated a multi-way fixed-effects fit, and both were doing work they already had the answer to.

fe_leverage() cross-tabulates each pair of factors. It built a composite index and handed it to rowsum(), which hashes all n observations to rediscover groups that are known to run 1..g, and then labels its result rows with those group values as strings, so the caller converted them back with as.integer(rownames(.)). At n = 100,000 over 1,000 by 30 groups that pair cost 7.3 ms, against 0.04 ms for a single counting pass, and there is one per factor pair. A dense weighted cross-tabulation in C++ (xtab_cpp) replaces it and serves the weighted case too, so the unit weight vector is no longer built for an unweighted fit. The masks and column offsets the leverage vector needs are now built only when a leverage vector is wanted, which is HC2 and HC3; every other se_type returns the rank without them.

demean_cpp() rebuilt the group weight sums on every sweep of the alternating projections, though they depend only on the codes and the weights and are identical each time, and then made a second full pass over the matrix to find its largest absolute value for the convergence test. The sums are now computed once, the maximum is accumulated by the subtraction pass that already touches every cell, and an unweighted fit no longer multiplies each element by a 1.0 it had just read out of a vector of ones. Every result is bit-identical; there is simply less of it.

Measured at n = 100,000 with two regressors, against 2.0.0 before this work: a plain fit 15.8 to 3.2 ms, a clustered fit 18.5 to 5.1 ms, a one-way fixed-effects fit 14.7 to 5.8 ms, a two-way fixed-effects fit 32.7 to 9.3 ms. On the cases that carry the most fixed-effects machinery: two-way HC2 36.7 to 12.8 ms, two-way weighted 35.0 to 12.6 ms.

All of it is in clean_model_data(), lm_robust_fit(), lm_return() and the fixed-effects helpers, so iv_robust() and lm_lin() get it for nothing: 2SLS 20.0 to 5.5 ms, clustered 2SLS 30.0 to 7.6 ms, 2SLS with two-way fixed effects 37.8 to 14.4 ms, lm_lin() 26.0 to 6.1 ms.

Every returned value is unchanged. Across 56 specifications covering each se_type, weights, clusters, multivariate outcomes, IV, and one to two fixed-effect factors, the largest relative difference from the previous version is 9e-15, and the names on every field are identical apart from the three naming bugs fixed below.

Fixed effects

Fixed effects absorption is implemented using the Frisch-Waugh-Lovell theorem with alternating projections (the algorithm used by fixest::feols). All FE variables are demeaned out of both the outcome and the design matrix before estimation; the intercept is absorbed into the group means and dropped.

SE types with fixed effects. With a single FE factor there is no restriction and no cost: the hat value of the full [dummies | X] design splits as h_ii = h_ii(demeaned X) + w_i / sum(w over i's group), so HC2 and HC3 are exact from the demeaned fit plus one vector. The default is therefore "HC2", the same as with no FE, and the same as 1.x returns.

lm_robust(y ~ z, data = dat, fixed_effects = ~block)
# se_type = "HC2", bit-identical to 1.x, 4.1 ms where 1.0.6 takes 41.0 s
# at n = 40,000 with 2,000 blocks

lm_robust(y ~ z, data = dat, fixed_effects = ~block, clusters = cl)
# se_type = "CR0", with a warning: 1.0.6 defaulted to "CR2" here

lm_robust(y ~ z, data = dat, fixed_effects = ~ block + year)
# se_type = "HC2", as in 1.0.6, and no warning: the identity covers any
# number of factors, so there is nothing to trade away

lm_robust(y ~ z, data = dat, fixed_effects = ~ block + year, se_type = "CR2",
          clusters = cl)
# se_type = "CR2", exact, by expanding the dummies. No warning: you asked.

With two or more FE factors the naive version of the identity fails, by roughly 3e-3 rather than by rounding: summing each factor’s own within-group share drops the cross terms. The projection decomposition P_[X | D] = P_D + P_{M_D X} is exact for any number of factors, and that is what is used. Only "CR2" with fixed_effects still expands the dummies, its adjustment being built from cluster-level blocks of the hat matrix rather than from h_ii. Every answer equals the explicit-dummy fit and equals 1.0.6; that equivalence is tested and recorded in the fixture.

What changes is only what a default will reach for, and only in the clustered case. Expanding the dummies is roughly O(g^3) in the number of levels, so no default expands them; "CR2" is the only estimator that still needs the expansion, so fixed_effects with clusters defaults to "CR0" with a warning, shown once per session, that names what 1.0.6 returned and names the se_type that accepts the new default. Writing se_type = "CR0" removes it. Unclustered fits keep the "HC2" default at any number of factors, as in 1.0.6.

Multi-way FE. Two-way (and higher) FE are supported via alternating projections that iterate to convergence (tolerance 1e-8, maximum 100 sweeps, after which it warns).

Return object additions. FE models add: - fes: logical flag indicating whether fixed effects were used - proj_r.squared: within-group R² from the demeaned regression - r.squared: full-model R² using the original (un-demeaned) outcome - fitted.values: reconstructed to the original scale via FWL (projected residuals equal full residuals)

Degrees of freedom. For a model with one FE variable having B levels, fe_rank = B degrees of freedom are consumed (B − 1 dummy contrasts plus the absorbed intercept), so df.residual = n − k − B where k is the number of non-FE regressors.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.