estimatr

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 groupe