lm_robust: Ordinary Least Squares with Robust Standard Errors

Description

Fits a linear model by ordinary least squares and returns heteroskedasticity-robust or cluster-robust standard errors, with the small-sample corrections used in design-based work. Fixed effects can be absorbed rather than expanded into dummy columns, at no cost in the available standard error types.

Usage

lm_robust(
  formula,
  data,
  weights,
  subset,
  clusters,
  fixed_effects,
  se_type = NULL,
  ci = TRUE,
  alpha = 0.05,
  return_vcov = TRUE,
  try_cholesky = FALSE
)

Value

An object of class "lm_robust", a list holding the estimate table in coefficients, std.error, df, statistic, p.value, conf.low, conf.high, term, and outcome; the fit in fitted.values, residuals, vcov, nobs, k, rank, df.residual, and res_var; the summary statistics r.squared, adj.r.squared, tss, and fstatistic; and se_type, weighted, clustered, fes, alpha, terms, xlevels, and call.

Absorbed fits add fixed_effects, felevels (the absorbed levels of each factor), and the within-projection summaries proj_r.squared, proj_adj.r.squared, proj_tss, and proj_fstatistic.

Arguments

formula

(required) An object of class formula, as in lm()

data

(optional) A data.frame

weights

(optional) The bare (unquoted) name of the weights variable

subset

(optional) A bare (unquoted) expression specifying a subset

clusters

(optional) A bare (unquoted) name of the cluster variable

fixed_effects

(optional) A one-sided formula of fixed effects to absorb rather than expand into dummy columns, such as ~ blockID or ~ block + year. Each variable is demeaned within the groups before OLS is run, so by the Frisch-Waugh-Lovell theorem the coefficients and residuals are the dummy regression's exactly.

Absorbing costs nothing in available standard error types. "HC2" and "HC3" are exact at any number of factors, because the leverage of the full design splits into the demeaned-X leverage plus a term that is cheap to compute, so no dummy hat matrix is built.

"CR2" is the exception: its adjustment is built from cluster-level blocks of the hat matrix rather than from the diagonal, and blocks do not split that way, so it expands the dummies and pays for the expansion. That is why fixed_effects with clusters defaults to "CR0". Asking for se_type = "CR2" still works and still gives the 1.0.6 number. Refused is the three together: "CR2" with both weights and fixed_effects, as in estimatr 1.0.6.

The projection identity, the several-factor case, the exact-rank calculation, and the weighted CR2 and HC2 conventions are derived in vignette("mathematical-notes").

se_type

(optional) The standard error type. Defaults depend on whether clusters and/or fixed effects are present:

No clusters, no FE: "HC2" (default), "HC0", "HC1", "HC3", "classical", "stata", "none".
Clusters, no FE: "CR2" (default), "CR0", "stata", "none".
No clusters, with FE (any number of factors): "HC2" (default), "HC0", "HC1", "HC3", "classical", "stata", "none". The same menu as with no FE at all.
Clusters, with FE: "CR0" (default), "CR2", "stata", "none". "CR2" expands the fixed effects into dummies, so it is not the default here; it is refused with weights.

"stata" means two different things. With no clusters it is exactly "HC1", and the fitted object reports se_type = "HC1". With clusters it is not an alias for "CR0": it is CR0 scaled by Stata's finite-sample factor, (J / (J - 1)) * ((N - 1) / (N - K)) on the variance, and the object reports se_type = "stata" to keep the distinction visible.

ci

(optional) Logical. Whether to compute p-values and confidence intervals.

alpha

(optional) The significance level, 0.05 by default.

return_vcov

(optional) Logical. Whether to return the vcov matrix.

try_cholesky

(optional) Logical. Whether to solve by Cholesky decomposition of X'X rather than by the default pivoted QR. FALSE by default, and worth turning on in most applied settings: about 1.4 times faster at n = 100,000 with two regressors, and 1.7 times faster at n = 200,000 with 60 regressors, where it is 0.15s against 0.25s. The saving is per fit, so it is worth most in a simulation that fits the same design thousands of times.

Rank deficiency is caught on either path. Redundant columns come back as NA exactly as they do from lm() whichever path ran, and a design that is rank deficient falls back to the QR.

Whether it is safe turns on one question, whether two regressors are nearly the same variable. Forming X'X squares the condition number, so the Cholesky path has about twice the rounding error of the QR, and only near-collinearity makes that visible. Differences of scale do not, because the columns are normalized before either decomposition, so a covariate in dollars beside one in years costs nothing. For a treatment indicator, a few covariates, block or cluster dummies, the centered interactions lm_lin() builds, or a factorial, the two paths agree to at least 10 significant digits, which is why difference_in_means() sets it to TRUE internally. Agreement falls to about 3 digits as the scaled condition index reaches 1e6, and the QR fallback takes over above roughly 1e8. Nothing interpretable lives in that range: a design at 1e6 returns a coefficient of 4.8e4 with a standard error of 4.6e4 on a regressor whose true effect is zero. To check a design directly, scale the columns first, since the unscaled condition number of a design in mixed units is large for a reason that does not affect the fit: kappa(sweep(X, 2, sqrt(colSums(X^2)), "/"), exact = TRUE).

Examples

Run this code

set.seed(15)
dat <- data.frame(
  y = rpois(40, lambda = 4),
  x = rnorm(40),
  z = rbinom(40, 1, prob = 0.4),
  cl = rep(1:10, each = 4),
  bl = rep(c("A", "B", "C", "D"), each = 10),
  w = runif(40)
)

# HC2 is the default
fit <- lm_robust(y ~ x + z, data = dat)
fit
tidy(fit)
summary(fit)
confint(fit, level = 0.8)

# Other variance estimators, including Stata's
lm_robust(y ~ x + z, data = dat, se_type = "classical")
lm_robust(y ~ x + z, data = dat, se_type = "stata")

# Clustered inference defaults to CR2
lm_robust(y ~ x + z, data = dat, clusters = cl)
lm_robust(y ~ x + z, data = dat, clusters = cl, se_type = "stata")

# Weights and subsets behave as they do in lm()
lm_robust(y ~ x + z, data = dat, weights = w, clusters = cl)
lm_robust(y ~ x, data = dat, subset = z == 1)

# Fixed effects are absorbed rather than expanded into dummies. With a
# single factor the HC2 default is exact and costs nothing extra.
lm_robust(y ~ z, data = dat, fixed_effects = ~ bl)

Run the code above in your browser using DataLab