Mathematical notes and guarantees for estimatr

Written with AI, and checked accordingly

estimatr 2.0.0 was written by Alexander Coppock working with Claude (Anthropic), across design, implementation, tests, benchmarks and documentation.

While the code base has been reviewed, it was not written by hand. The guarantee offered here is therefore not that every line has been vouched for. It is narrower and it is checkable: estimatr implements these estimators correctly, and the way that is shown is validation to machine precision against the definitions themselves.

Which is what this document is. Every estimator gets its definition, stated in mathematics with the paper it comes from, and then, immediately underneath, the same quantity computed twice: once by calling estimatr, and once from the definition transcribed into a few lines of base R.

Two further layers are checked outside this document. Every estimator reproduces estimatr 1.0.6’s numbers wherever both versions answer, checked in tests/testthat/test_vs_estimatr.R against 695 values recorded from an installed 1.0.6. A further 808 assertions compare against implementations that share no lineage with this one: sandwich, clubSandwich, ivreg, Stata, fixest, plm and blkvar, in the five tests/testthat/test_vs_*.R files. vignette("estimatr2.0") sets out both layers under “How this was checked”; the suite holds 5,635 assertions in total.

How the checking works

An identity holds to machine precision or it is broken. Nothing below is random and nothing is replicated, so there is no sampling error to allow for and no tolerance to argue about. The two quantities being compared are the same number, and what gets reported is the largest relative gap between them, expected to sit near the floor of double-precision arithmetic.

The reference side is written in this document rather than borrowed from another package, on purpose. A comparison against sandwich shows that two implementations agree. A comparison against the formula shows what the estimator is, which is the question a reader of mathematical notes is actually asking. It also leaves the document depending on nothing but estimatr, so no check can vanish because a suggested package is missing.

CHECKS <- list()

check <- function(label, ours, theirs, tol = 1e-10) {
  gap <- max(abs(ours - theirs) / pmax(abs(theirs), 1))
  # Two jobs: record the gap in the running list for the final table, and
  # return a one-row data frame so the calling chunk prints its own result.
  CHECKS[[label]] <<- gap
  data.frame(gap = sprintf("%.1e", gap), holds = gap < tol)
}

Each section calls check() once, prints its own result, and adds it to a running list. Every promise in one table collects them at the end and the document refuses to build if any of them fails.

Notation

Throughout, $\mathbf{X}$ is the $N \times K$ design matrix, $\mathbf{y}$ the outcome, $\mathbf{e} = \mathbf{y} - \mathbf{X}\widehat{\beta}$ the residuals, and $\mathbf{x}_i$ the $i$th row of $\mathbf{X}$. $\mathbf{W}$ is a diagonal matrix of weights scaled to sum to one, and $\mathrm{diag}[\cdot]$ builds a diagonal matrix from a vector. For clustered designs, $S$ is the number of clusters and $\mathbf{X}_s$ and $\mathbf{e}_s$ are the rows belonging to cluster $s$. For blocked designs, $J$ is the number of blocks and $N_j$ the size of block $j$.

The data used throughout

One hundred units, a binary treatment, a covariate, weights, ten groups for the fixed-effects section, twenty clusters, and an instrument with the endogenous regressor it shifts. Drawn once, at a fixed seed, and reused by every check below.

set.seed(20260826)
N <- 100

d <- data.frame(
  z  = rbinom(N, 1, 0.5),
  x  = rnorm(N),
  w  = runif(N, 0.5, 2),
  g  = rep(1:10, each = 10),
  cl = rep(1:20, each = 5),
  inst = rnorm(N)
)
d$en <- d$inst + rnorm(N, 0, 0.5)
d$y <- 1 + 0.5 * d$z + 2 * d$x + d$en + rnorm(N)

lm_robust

Coefficients

$$\widehat{\beta} = (\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{y}$$

The solver is a rank-revealing column-pivoting QR factorization from the Eigen C++ library, reached through RcppEigen, so $(\mathbf{X}^{\top}\mathbf{X})^{-1}$ is never formed explicitly. On a rank-deficient design the pivoting can drop a different column than lm() drops; the fitted values and the variance are the same either way, but which coefficient comes back NA may differ. Unlike 1.x, estimatr names the dropped terms in a warning rather than leaving them to be noticed in the output. Setting try_cholesky = TRUE substitutes a Cholesky factorization, which is faster and is guaranteed only when $\mathbf{X}$ has full rank.

The promise: the point estimates are least squares. lm() is the reference, since it solves the same problem by a different factorization.

check("lm_robust()",
      coef(lm_robust(y ~ z + x, data = d)),
      coef(lm(y ~ z + x, data = d)))
#>       gap holds
#> 1 4.6e-16  TRUE

Weights

Weights are scaled to sum to one, then each row of the design matrix and each outcome are multiplied by $\sqrt{w_i}$. Estimation proceeds on the transformed data, which gives

$$\widehat{\beta} = (\mathbf{X}^{\top}\mathbf{W}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{W}\mathbf{y}.$$

Romano and Wolf (2017) set out the properties that recommend the estimator. Everything below applies to the transformed data, so $(\mathbf{X}^{\top}\mathbf{X})^{-1}$ should be read as $(\mathbf{X}^{\top}\mathbf{W}\mathbf{X})^{-1}$ and $\mathbf{X}$ as $\mathbf{W}^{1/2}\mathbf{X}$ wherever weights are in play.

A row with weight zero contributes nothing to the fit and is not counted as an observation in the residual degrees of freedom or in the HC1 and "stata" scale factors, which is how lm() counts it too. The row is still returned in residuals and fitted.values.

The promise: the weighted fit is weighted least squares.

check("lm_robust(weights = )",
      coef(lm_robust(y ~ z + x, data = d, weights = w)),
      coef(lm(y ~ z + x, data = d, weights = w)))
#>       gap holds
#> 1 2.2e-16  TRUE

Heteroskedasticity-robust variance

The default is HC2, from MacKinnon and White (1985). It is the choice that lines up with design-based inference: under complete randomization the HC2 variance of a treatment coefficient equals the conservative Neyman estimator (Samii and Aronow 2012). Against the HC1 variance that Stata defaults to it gives up a little efficiency in large samples and is better behaved in small ones, which is the reason for the default.

se_type	$\widehat{\mathbb{V}}[\widehat{\beta}]$	Degrees of freedom
"classical"	$\frac{\mathbf{e}^\top\mathbf{e}}{N-K}(\mathbf{X}^{\top}\mathbf{X})^{-1}$	$N-K$
"HC0"	$(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathrm{diag}\left[e_i^2\right]\mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1}$	$N-K$
"HC1", "stata"	$\frac{N}{N-K}(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathrm{diag}\left[e_i^2\right]\mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1}$	$N-K$
"HC2" (default)	$(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathrm{diag}\left[\frac{e_i^2}{1-h_{ii}}\right]\mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1}$	$N-K$
"HC3"	$(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathrm{diag}\left[\frac{e_i^2}{(1-h_{ii})^2}\right]\mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1}$	$N-K$

where $h_{ii} = \mathbf{x}_i(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{x}_i^{\top}$ is the $i$th leverage value. Long and Ervin (2000) review the family and its small-sample behaviour.

Transcribed, the four robust members of that column are one function. bread is $(\mathbf{X}^{\top}\mathbf{X})^{-1}$, h the leverage diagonal, and adj the bracketed term that distinguishes them.

hc_vcov <- function(fit, type) {
  X <- model.matrix(fit)
  e <- residuals(fit)
  bread <- solve(crossprod(X))
  h <- rowSums((X %*% bread) * X)
  n <- nrow(X)
  k <- ncol(X)
  adj <- switch(type,
    HC0 = e^2,
    HC1 = e^2 * n / (n - k),
    HC2 = e^2 / (1 - h),
    HC3 = e^2 / (1 - h)^2
  )
  bread %*% crossprod(X * sqrt(adj)) %*% bread
}

The promise: the classical variance is the textbook one, and each robust variance is its own row of that table.

fit_lm <- lm(y ~ z + x, data = d)

check("lm_robust(se_type = 'classical')",
      lm_robust(y ~ z + x, data = d, se_type = "classical")$vcov,
      vcov(fit_lm))
#>       gap holds
#> 1 2.8e-17  TRUE

do.call(rbind, lapply(c("HC0", "HC1", "HC2", "HC3"), function(ty) {
  cbind(se_type = ty,
        check(paste0("lm_robust(se_type = '", ty, "')"),
              lm_robust(y ~ z + x, data = d, se_type = ty)$vcov,
              hc_vcov(fit_lm, ty)))
}))
#>   se_type     gap holds
#> 1     HC0 6.9e-17  TRUE
#> 2     HC1 1.1e-16  TRUE
#> 3     HC2 5.6e-17  TRUE
#> 4     HC3 8.3e-17  TRUE

Leverage at and above one

HC2 and HC3 divide by $1 - h_{ii}$, so a leverage value at or above one is a special case rather than an ordinary one.

Leverage exactly equal to one is benign. The residual is exactly zero, the contribution is a $0/0$ that resolves to zero, and the standard error is finite. Leverage marginally above one is not benign, and it happens: a near-saturated design can compute $h_{ii} = 1 + 10^{-16}$, at which point $1 - h_{ii}$ is negative. Under HC3 that row contributes a negative term to a variance. Under HC2 the implementation takes a square root of it, so a single such row turns every standard error in the fit into NaN, however small the offending quantity.

estimatr 2.0 sets the contribution of any row with $1 - h_{ii} \le 0$ to zero and warns, naming how many rows were affected. estimatr 1.0.6 returned NaN for HC2 and a silently inflated number for HC3 on the same designs. The CR2 estimator below has no analogous hole: it never forms $1 - h_{ii}$, and the eigenvalue clamp described there covers the degenerate case.

Note what the check above does and does not cover. d is well conditioned, with a hundred observations and three parameters, so no leverage in it comes near one. The degenerate designs are checked in the suite, not here, which is a limit this document shares with any table built on rnorm().

Cluster-robust variance

The cluster-robust estimators are the analogues of the heteroskedasticity-consistent ones. The default is CR2, from Bell and McCaffrey (2002), in the generalized form of Pustejovsky and Tipton (2018), whose clubSandwich package applies the same correction across a wider range of models. Imbens and Kolesár (2016) compare the alternatives in small samples.

se_type	$\widehat{\mathbb{V}}[\widehat{\beta}]$	Degrees of freedom
"CR0"	$(\mathbf{X}^{\top}\mathbf{X})^{-1}\sum_{s=1}^{S}\left[\mathbf{X}_s^\top\mathbf{e}_s\mathbf{e}_s^\top\mathbf{X}_s\right](\mathbf{X}^{\top}\mathbf{X})^{-1}$	$S-1$
"stata"	$\frac{N-1}{N-K}\frac{S}{S-1}\times$ the CR0 expression	$S-1$
"CR2" (default)	$(\mathbf{X}^{\top}\mathbf{X})^{-1}\sum_{s=1}^{S}\left[\mathbf{X}_s^\top\mathbf{A}_s\mathbf{e}_s\mathbf{e}_s^\top\mathbf{A}_s^\top\mathbf{X}_s\right](\mathbf{X}^{\top}\mathbf{X})^{-1}$	Satterthwaite, below

Transcribed, CR0 is the same bread with the meat summed over clusters instead of over observations, and "stata" is CR0 times two finite-sample corrections.

cr_vcov <- function(fit, cluster, stata = FALSE) {
  X <- model.matrix(fit)
  e <- residuals(fit)
  bread <- solve(crossprod(X))
  meat <- Reduce(`+`, lapply(split(seq_len(nrow(X)), cluster), function(i) {
    tcrossprod(crossprod(X[i, , drop = FALSE], e[i]))
  }))
  v <- bread %*% meat %*% bread
  if (!stata) return(v)
  S <- length(unique(cluster))
  v * (S / (S - 1)) * ((nrow(X) - 1) / (nrow(X) - ncol(X)))
}

The promise: CR0 is the cluster sandwich, and "stata" is CR0 times Stata’s two corrections.

check("lm_robust(clusters = )",
      lm_robust(y ~ z + x, data = d, clusters = cl, se_type = "CR0")$vcov,
      cr_vcov(fit_lm, d$cl))
#>       gap holds
#> 1 1.1e-16  TRUE

check("lm_robust(se_type = 'stata')",
      lm_robust(y ~ z + x, data = d, clusters = cl, se_type = "stata")$vcov,
      cr_vcov(fit_lm, d$cl, stata = TRUE))
#>       gap holds
#> 1 1.1e-16  TRUE

CR2 is the one member of the family whose reference is not a few lines of base R, so it is checked in the suite against clubSandwich instead, live and at $10^{-10}$. The adjustment matrices come from

$$\begin{aligned} \mathbf{H} &= \mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^\top \\ \mathbf{B}_s &= (\mathbf{I}_N - \mathbf{H})_s (\mathbf{I}_N - \mathbf{H})_s^\top \\ \mathbf{A}_s &= \mathbf{B}_s^{+1/2} \end{aligned}$$

where $(\mathbf{I}_N - \mathbf{H})_s$ are the $N_s$ columns belonging to cluster $s$ and $\mathbf{B}_s^{+1/2}$ is the symmetric square root of the Moore-Penrose inverse. estimatr reaches that inverse through an eigendecomposition with eigenvalues clamped below $10^{-12}$, which is what lets a rank-deficient cluster (fixed effects that coincide with the clusters, for instance) return an answer where the Bell and McCaffrey (2002) form could not be computed at all. The two forms agree whenever $\mathbf{B}_s$ has full rank.

The degrees of freedom are computed per coefficient:

$$\mathrm{df}_k = \frac{\left(\sum_{s=1}^{S}\mathbf{p}_s^\top\mathbf{p}_s\right)^2}{\sum_{s=1}^{S}\sum_{t=1}^{S}\left(\mathbf{p}_s^\top\mathbf{p}_t\right)^2}, \qquad \mathbf{p}_s = (\mathbf{I}_N - \mathbf{H})_s^\top\mathbf{A}_s\mathbf{X}_s(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{z}_k$$

with $\mathbf{z}_k$ the $k$th standard basis vector. Different coefficients in one fit can therefore carry different degrees of freedom.

Under weights, CR2 and HC2 follow different conventions, and the difference is invisible at the call site. CR2’s small-sample adjustment is built against a working model with identity covariance, $\mathbf{\Phi} = \mathbf{I}$, where the weighted HC2 adjustment is built against precision weights. In clubSandwich’s terms the weighted CR2 here is vcovCR(..., inverse_var = FALSE) and the weighted HC2 is inverse_var = TRUE. Each is internally consistent; they are not the same convention as one another, and the choice is inherited from estimatr 1.0.6 rather than made here. Both halves are pinned explicitly in tests/testthat/test_vs_clubsandwich.R, with inverse_var named on each side so that a change in clubSandwich’s default fails the test rather than quietly asserting the other convention.

One cluster is refused. A cluster-robust variance needs variation across clusters. Given a single cluster, estimatr 1.0.6 returned a standard error of about $6 \times 10^{-17}$ and a confidence interval of zero width, in silence. estimatr 2.0 raises an error.

Absorbed fixed effects

fixed_effects = ~ g partials the dummies for g out of the outcome and the covariates rather than adding them as columns. Point estimates are identical to the dummy regression by the Frisch-Waugh-Lovell theorem. The variance is where the work is, because HC2 and HC3 are built from the leverage values of the full design, the one with every dummy in it, which absorbing is precisely the decision not to build.

The way out is an identity. Write $\mathbf{D}$ for the matrix of fixed-effect dummies and $\mathbf{M_D} = \mathbf{I} - \mathbf{P_D}$ for the residual-maker that demeans. The projection onto the full design splits exactly:

$$\mathbf{P}_{[\mathbf{X}\,|\,\mathbf{D}]} = \mathbf{P_D} + \mathbf{P}_{\mathbf{M_D}\mathbf{X}}$$

so each leverage value of the full design is the leverage value of the demeaned covariates, which the fitter already has, plus the $i$th diagonal element of $\mathbf{P_D}$, which is cheap:

• One factor. $\mathbf{P_D}$ is diagonal and the term is the unit’s share of its own group’s weight, $w_i / \sum_{j \in g(i)} w_j$. Unweighted, that is one over the group’s size.
• Several factors. Write $\mathbf{D}$ with the widest factor in full dummies and the rest contrast-coded. The leading block of $\mathbf{D}^\top\mathbf{W}\mathbf{D}$ is then diagonal, so block inversion needs only the Schur complement, a matrix of side $\sum_{k>1}(g_k - 1)$: the design’s narrowest dimension rather than its widest. No dummy matrix is built at any number of factors.

The identity holds for any number of factors, so HC2 and HC3 carry no restriction under fixed_effects. CR2 is the exception: its correction comes from cluster-level blocks of the hat matrix rather than from the diagonal, and blocks do not decompose this way, so CR2 still expands the dummies. That cost, roughly cubic in the number of levels, is why fixed_effects combined with clusters defaults to CR0 in 2.0 where 1.x defaulted to CR2. It is the only default that moved in the release, it warns once per session, and naming se_type = "CR2" still gets the 1.x number exactly.

The promise: absorbing a factor changes the speed, not the answer. The reference is the dummy regression the absorption is supposed to reproduce, coefficients and standard errors alike.

absorbed <- lm_robust(y ~ z + x, data = d, fixed_effects = ~ g)
dummies  <- lm_robust(y ~ z + x + factor(g), data = d)
keep <- c("z", "x")

check("lm_robust(fixed_effects = )",
      c(coef(absorbed)[keep], absorbed$std.error[keep]),
      c(coef(dummies)[keep], dummies$std.error[keep]))
#>       gap holds
#> 1 7.8e-16  TRUE

Rank

The Schur complement is inverted through its eigendecomposition rather than a solve, which does two things at once. A disconnected or nested fixed-effect design makes $\mathbf{D}$ rank deficient, and the pseudo-inverse returns the right projection anyway. The eigenvalues also give the exact rank of the fixed-effect design for free,

$$\mathrm{rank}(\mathbf{D}) = g_1 + \mathrm{rank}(\mathbf{S}),$$

which is what estimatr uses for the residual degrees of freedom. The nominal count $\sum_k g_k - K + 1$ overstates the rank whenever one factor is partly spanned by the others, and 1.x used the nominal count, so its absorbed fit disagreed with its own explicit-dummy fit on such designs. lm() and plm report the exact rank; fixest reports the nominal one unless asked for ssc(K.exact = TRUE).

Weights and the hat matrix: where Stata differs

With weights, lm_robust()’s HC2 and HC3 standard errors do not match Stata’s vce(hc2) and vce(hc3). The cause is a difference in how the hat matrix is defined, and the choice is a convention rather than an error on either side. It is the one place in this document where a definition is contested, so there is no identity to check and the section reports a disagreement instead.

Stata uses

$$\mathbf{H}_{\text{Stata}} = \mathbf{X}(\mathbf{X}^{\top}\mathbf{W}\mathbf{X})^{-1}\mathbf{X}^\top$$

while estimatr, sandwich, and Python’s statsmodels all use

$$\mathbf{H}_{R} = \mathbf{X}(\mathbf{X}^{\top}\mathbf{W}\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{W}.$$

Only HC2 and HC3 depend on the hat matrix, so the divergence is confined to those two. Weighted classical, HC0, HC1 and the clustered "stata" estimator all agree with Stata exactly, and the test suite pins both facts: the weighted HC2 and HC3 variances differ from Stata’s by a bounded amount, under 2 percent on the reference fits, while the weighted HC1 and clustered fits match to the precision Stata printed.

Two arguments favour $\mathbf{H}_R$. It is what you get by rescaling the data by $\sqrt{w_i}$ and running ordinary least squares, so it follows if you regard the weighted model as a rescaling of the unweighted one. Its diagonal elements are also the weighted leverages in the sense of Li and Valliant (2009), where $\mathbf{H}_{\text{Stata}}$ would have to be weighted a second time to recover them. Against that, Stata’s convention has the weight of Stata behind it, and the differences are small. The choice is genuinely open, which is why estimatr pins it from both sides in the suite rather than treating either answer as the error.

lm_robust(mpg ~ hp, data = mtcars, weights = wt, se_type = "HC2")$std.error
#> (Intercept)          hp 
#>     2.16282     0.01446

Stata 13 reports 0.0143083 on hp for the same fit, about one percent below the number above. Python’s statsmodels returns estimatr’s. Change se_type to "HC1" and Stata and estimatr agree exactly.

Confidence intervals and testing

With $\widehat{\mathbb{V}}_k$ the $k$th diagonal element of $\widehat{\mathbb{V}}$,

$$\mathrm{CI}^{1-\alpha} = \left(\widehat{\beta}_k + t^{\mathrm{df}}_{\alpha/2}\sqrt{\widehat{\mathbb{V}}_k},\; \widehat{\beta}_k + t^{\mathrm{df}}_{1-\alpha/2}\sqrt{\widehat{\mathbb{V}}_k}\right)$$

and two-sided p-values come from the same $t$ distribution. Under CR2 the degrees of freedom vary by coefficient, so the multiplier does too.

lm_lin

lm_lin() is a pre-processor for lm_robust() implementing the covariate adjustment of Lin (2013), which answers Freedman (2008)’s demonstration that regression adjustment can reduce precision. Rather than

$$y_i = \tau z_i + \mathbf{\beta}^\top\mathbf{x}_i + \epsilon_i,$$

it centers every covariate at its sample mean and interacts the centered covariates with treatment:

$$y_i = \tau z_i + \mathbf{\beta}^\top\mathbf{x}^c_i + \mathbf{\gamma}^\top\mathbf{x}^c_i z_i + \epsilon_i.$$

Centering is what makes $\tau$ the estimate of the average treatment effect: at $\mathbf{x}^c = \mathbf{0}$ the interaction terms drop out. Centering happens after any function in the covariates formula is evaluated, so ~ log(x) centers $\log(x)$ rather than the log of the centered $x$. The centers are returned in scaled_center.

Multi-valued treatments are handled by building a full set of dummies and interacting each with the centered covariates. Everything else, weights, clusters, se_type, is lm_robust()’s.

The promise: it is the Lin specification a user could write by hand. The reference is that specification, written by hand.

dd <- d
dd$x_c <- dd$x - mean(dd$x)
lin  <- lm_lin(y ~ z, covariates = ~ x, data = d)
byhand <- lm_robust(y ~ z * x_c, data = dd)

check("lm_lin()",
      c(coef(lin)[["z"]], lin$std.error[["z"]]),
      c(coef(byhand)[["z"]], byhand$std.error[["z"]]))
#>       gap holds
#> 1 0.0e+00  TRUE

iv_robust

Coefficients

$$\widehat{\beta}_{2SLS} = (\mathbf{X}^{\top}\mathbf{P_Z}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{P_Z}\mathbf{y}, \qquad \mathbf{P_Z} = \mathbf{Z}(\mathbf{Z}^{\top}\mathbf{Z})^{-1}\mathbf{Z}^\top$$

with $\mathbf{X}$ the regressors, endogenous ones included, and $\mathbf{Z}$ the instruments. Equivalently: regress $\mathbf{X}$ on $\mathbf{Z}$ to get $\widehat{\mathbf{X}} = \mathbf{Z}\widehat{\beta}_{FS}$, then regress $\mathbf{y}$ on $\widehat{\mathbf{X}}$. Weights are handled as in lm_robust(), by rescaling before estimation.

The promise: the point estimates are two-stage least squares. The reference is the two stages, run as two stages.

tsls_coef <- function(y, X, Z) {
  xhat <- Z %*% solve(crossprod(Z), crossprod(Z, X))
  as.vector(solve(crossprod(xhat), crossprod(xhat, y)))
}

check("iv_robust()",
      unname(coef(iv_robust(y ~ en + x | inst + x, data = d))),
      tsls_coef(d$y, model.matrix(~ en + x, d), model.matrix(~ inst + x, d)))
#>       gap holds
#> 1 8.5e-16  TRUE

Variance

The variance estimators are lm_robust()’s with two substitutions. The second-stage regressors $\widehat{\mathbf{X}}$ replace $\mathbf{X}$, and the residuals are $\mathbf{y} - \mathbf{X}\widehat{\beta}_{2SLS}$, formed from the endogenous, uninstrumented regressors rather than from the fitted ones. residuals() returns those structural residuals, not the first-stage ones.

Which leverage

HC2 and HC3 need leverage values, and 2SLS admits two candidates. estimatr uses the second-stage hat values,

$$h_i = \widehat{\mathbf{x}}_i(\widehat{\mathbf{X}}^{\top}\widehat{\mathbf{X}})^{-1}\widehat{\mathbf{x}}_i^{\top},$$

the diagonal of an orthogonal projection. The alternative is the diagonal of $\mathbf{H}^{*}$, the matrix carrying $\mathbf{y}$ to its fitted values. Belsley, Kuh, and Welsch (1980) considered it, observed that $\mathbf{H}^{*}$ is idempotent but not symmetric, and recommended the second-stage hat values on the ground that the diagonal of an asymmetric matrix is not a leverage.

The choice has consequences. A projection diagonal lies in $[0,1]$, so HC2 is always defined. The diagonal of $\mathbf{H}^{*}$ is already negative for one row of mtcars, and across 3,000 weak-first-stage designs it exceeded one in 10.8 percent of them, reaching 309.

The ivreg package makes the second-stage convention its default, and sandwich::vcovHC() applied to an ivreg::ivreg() fit returns estimatr’s standard errors to machine precision. sandwich has no leverage convention of its own; it calls hatvalues() on whatever fit it is given. AER::ivreg()’s hatvalues method predates ivreg and returns $\mathrm{diag}(\mathbf{H}^{*})$, so estimatr differs from AER, by up to 8.6 percent at HC2 and 18.5 percent at HC3 on mtcars, and agrees with the successor package that deprecates AER’s method. The numbers are bit-identical to estimatr 1.0.6.

Correspondence with Stata

Stata’s ivregress 2sls applies no finite-sample correction and uses z-tests unless told otherwise.

estimatr	Stata
no equivalent	ivregress 2sls y (x = z)
se_type = "classical"	ivregress 2sls y (x = z), small
se_type = "HC0"	ivregress 2sls y (x = z), rob
se_type = "HC1"	ivregress 2sls y (x = z), rob small
clusters = cl, se_type = "CR0"	ivregress 2sls y (x = z), vce(cl cl)
clusters = cl, se_type = "stata"	ivregress 2sls y (x = z), vce(cl cl) small
se_type = "HC2" (default), "HC3", "CR2"	no equivalent

lh_robust

lh_robust() fits a model with lm_robust() and then tests linear restrictions on it through car::linearHypothesis(), keeping the robust variance and the degrees of freedom of the fit rather than recomputing them classically.

For a restriction vector $\mathbf{a}$, the estimate and its standard error are the delta method applied to a linear function of the coefficients:

$$\widehat{\theta} = \mathbf{a}^\top\widehat{\beta}, \qquad \mathrm{se}(\widehat{\theta}) = \sqrt{\mathbf{a}^\top\widehat{\mathbb{V}}\mathbf{a}}$$

with $\widehat{\mathbb{V}}$ whichever variance the fit was asked for. Several restrictions at once, stacked into a matrix $\mathbf{R}$ against targets $\mathbf{q}$, additionally give a Wald statistic

$$F = \frac{(\mathbf{R}\widehat{\beta} - \mathbf{q})^\top\left[\mathbf{R}\widehat{\mathbb{V}}\mathbf{R}^\top\right]^{-1}(\mathbf{R}\widehat{\beta} - \mathbf{q})}{\mathrm{rank}(\mathbf{R})}$$

on $\mathrm{rank}(\mathbf{R})$ and the fit’s residual degrees of freedom, returned in the joint_hypothesis element. estimatr 1.x declines to compute it.

The promise: a linear hypothesis is the delta method on the fit.

fit <- lm_robust(y ~ z + x, data = d)
lh <- lh_robust(y ~ z + x, data = d, linear_hypothesis = "z + x = 0")
a <- c(0, 1, 1)

check("lh_robust()",
      c(lh$lh$coefficients[[1]], lh$lh$std.error[[1]]),
      c(sum(a * coef(fit)), sqrt(drop(t(a) %*% fit$vcov %*% a))))
#>       gap holds
#> 1 0.0e+00  TRUE

difference_in_means

difference_in_means() picks the point estimate, variance, and degrees of freedom that match the design, and reports which one it used in the design element of the fitted object. The design is inferred from which of blocks and clusters are supplied, and from the shape of the blocks.

Estimates

Unblocked.

$$\widehat{\tau} = \frac{1}{N_1}\sum_{i: z_i = 1} y_i \;-\; \frac{1}{N_0}\sum_{i: z_i = 0} y_i$$

Blocked. The sample-weighted average of the within-block estimates,

$$\widehat{\tau} = \sum_{j=1}^{J}\frac{N_j}{N}\widehat{\tau}_j.$$

With weights, the estimate and its variance are handed to lm_robust() with HC2 standard errors, within each block if the design is blocked.

Variance for unblocked and clustered designs

Design	$\widehat{\mathbb{V}}[\widehat{\tau}]$	Degrees of freedom
No blocks, no clusters	$\frac{\widehat{\mathbb{V}}[y_{i,0}]}{N_0} + \frac{\widehat{\mathbb{V}}[y_{i,1}]}{N_1}$	Welch-Satterthwaite
Clusters, no blocks	the CR2 estimator of lm_robust()	as CR2
Blocked and clustered	$\sum_j \left(\frac{N_j}{N}\right)^2\widehat{\mathbb{V}}[\widehat{\tau}_j]$	$S - 2J$
Matched-pair clustered	$\frac{J}{(J-1)N^2}\sum_j\left(N_j\widehat{\tau}_j - \frac{N\widehat{\tau}}{J}\right)^2$	$J-1$

The unblocked variance and its degrees of freedom are what R’s t.test() computes. The clustered variance is the one Gerber and Green (2012) recommend in their equation 3.23 when clusters are of even size. The matched-pair clustered variance is the SATE variance of Imai, King, and Nall (2009), their equation 6, with the degrees of freedom they suggest.

That first row has a second description: the Neyman variance of a two-arm experiment is exactly what HC2 returns on a regression of the outcome on the treatment indicator, which is Samii and Aronow (2012)’s equivalence and the reason HC2 is the package default.

The promise: for a two-arm design it is lm_robust() at HC2.

dim_fit <- difference_in_means(y ~ z, data = d)
ols <- lm_robust(y ~ z, data = d, se_type = "HC2")

check("difference_in_means()",
      c(dim_fit$coefficients[["z"]], dim_fit$std.error[["z"]]),
      c(ols$coefficients[["z"]], ols$std.error[["z"]]))
#>       gap holds
#> 1 4.9e-16  TRUE

Variance for blocked designs

Blocked designs are where estimatr 2.0 departs most from 1.x, and the estimators come from Pashley and Miratrix (2021).

The classification is by arm counts, not by block size. A block with at least two treated and at least two control units has an estimable within-block variance and carries its own Neyman variance. A block with a singleton arm, one treated unit or one control unit, does not: with a single observation in an arm there is nothing to take a variance of. The variation across such blocks stands in for the variance they cannot each supply, which is the logic that makes the matched-pairs estimator work.

Write $\mathcal{B}$ for the blocks with both arms of size two or more, $\mathcal{S}$ for the blocks with a singleton arm, $n_{\mathcal{B}} = \sum_{j \in \mathcal{B}} N_j$ and $n_{\mathcal{S}} = \sum_{j \in \mathcal{S}} N_j$.

The estimable part is the usual blocked variance over $\mathcal{B}$ alone (their equation 4):

$$\widehat{\mathbb{V}}_{\mathcal{B}} = \frac{1}{n_{\mathcal{B}}^2}\sum_{j \in \mathcal{B}} N_j^2\,\widehat{\mathbb{V}}[\widehat{\tau}_j], \qquad \mathrm{df}_{\mathcal{B}} = n_{\mathcal{B}} - 2|\mathcal{B}|.$$

The singleton part is estimated across blocks. If every block in $\mathcal{S}$ is the same size, the estimator is the familiar matched-pairs one (their equation 5),

$$\widehat{\mathbb{V}}_{\mathcal{S}} = \frac{1}{|\mathcal{S}|(|\mathcal{S}|-1)}\sum_{j \in \mathcal{S}}\left(\widehat{\tau}_j - \bar{\tau}_{\mathcal{S}}\right)^2 ,$$

with $\bar{\tau}_{\mathcal{S}} = \sum_{j \in \mathcal{S}} N_j\widehat{\tau}_j / n_{\mathcal{S}}$. If the blocks differ in size, their equation 8 handles it without requiring any two blocks to match:

$$\widehat{\mathbb{V}}_{\mathcal{S}} = \frac{\sum_{j \in \mathcal{S}} \omega_j \left(\widehat{\tau}_j - \bar{\tau}_{\mathcal{S}}\right)^2}{n_{\mathcal{S}} + \sum_{j \in \mathcal{S}}\omega_j}, \qquad \omega_j = \frac{N_j^2}{n_{\mathcal{S}} - 2N_j},$$

with $\mathrm{df}_{\mathcal{S}} = |\mathcal{S}| - 1$ in both cases. The equal-size form is kept separate because equation 8 is undefined at two equal-sized blocks.

Combining. A design holding both kinds of block is the hybrid of their section 3.3, and the two parts combine by squared share of the sample:

$$\widehat{\mathbb{V}}[\widehat{\tau}] = \left(\frac{n_{\mathcal{B}}}{N}\right)^2\widehat{\mathbb{V}}_{\mathcal{B}} + \left(\frac{n_{\mathcal{S}}}{N}\right)^2\widehat{\mathbb{V}}_{\mathcal{S}}.$$

The paper stops at the variance. estimatr combines the two degrees-of-freedom components by Welch-Satterthwaite, which reduces to $N - 2J$ when every block is estimable and to $J - 1$ when every block has a singleton arm, matching what each literature uses on its own.

design reports which case applied.

blocked <- data.frame(bl = rep(1:10, each = 10),
                      z = rep(rep(0:1, each = 5), times = 10))
blocked$y <- rnorm(100) + 0.3 * blocked$z
difference_in_means(y ~ z, data = blocked, blocks = bl)$design
#> [1] "Blocked"

pairs <- data.frame(bl = rep(1:50, each = 2), z = rep(c(0, 1), 50))
pairs$y <- rnorm(100) + 0.3 * pairs$z
difference_in_means(y ~ z, data = pairs, blocks = bl)$design
#> [1] "Matched-pair"

# Both kinds of block in one design: 1.x applied the matched-pairs estimator
# to all of it, after a warning.
hybrid <- rbind(blocked, transform(pairs, bl = bl + 100))
difference_in_means(y ~ z, data = hybrid, blocks = bl)$design
#> [1] "Hybrid blocked"

What is refused

Two blocked designs are errors rather than estimates, because the variance genuinely cannot be estimated.

Exactly one block with a singleton arm. The variation across such blocks is what stands in for their within-block variance, and one block has no variation to offer.

Singleton-arm blocks of different sizes where one holds half or more of their units. Equation 8’s weights $\omega_j = N_j^2/(n_{\mathcal{S}} - 2N_j)$ require $N_j < n_{\mathcal{S}}/2$, which is what keeps them positive and the estimator conservative.

Both errors suggest merging blocks or using lm_robust() with block fixed effects.

Blocks of clusters are separate. Pashley and Miratrix (2021) treat treatment assigned to units within blocks, not to clusters within blocks, so blocked designs that also specify clusters use the earlier estimators, and every block must hold at least two treated and two control clusters unless the design is matched-pair clustered. A block with a single treated or control cluster is refused: its within-block variance is not estimable, and estimating it anyway understates the standard error by roughly the block’s cluster count.

horvitz_thompson

horvitz_thompson() estimates the average treatment effect by inverse probability weighting, which is unbiased when the assignment probabilities are known. Aronow and Middleton (2013), Middleton and Aronow (2015) and Aronow and Samii (2017) develop the estimator and its variance.

Let $\pi_{zi}$ be the marginal probability that unit $i$ is assigned to condition $z$, and $\pi_{zi,wj}$ the joint probability that unit $i$ is in condition $z$ and unit $j$ in condition $w$. Write

$$\widetilde{Y}_{zi} = \frac{y_i}{\pi_{zi}}$$

for the inverse-probability-weighted outcome of a unit observed in condition $z$.

Estimates

$$\widehat{\tau} = \frac{1}{N}\left(\sum_{i: z_i = 1}\widetilde{Y}_{1i} - \sum_{i: z_i = 0}\widetilde{Y}_{0i}\right)$$

$N$ is the number of units the design covers, which matters with more than two arms. condition1 and condition2 select the contrast, but the estimand remains the average treatment effect over every unit of the design, so the estimator divides by $N$ rather than by the number of units landing in the two selected conditions, and data must carry one row per unit including the arms outside the contrast. A declaration whose size does not match nrow(data) is an error rather than a silent misalignment.

The promise: the estimate is the Horvitz-Thompson estimator. Two lines is the whole definition.

ht_estimate <- function(y, z, pr) {
  mean(y * z / pr) - mean(y * (1 - z) / (1 - pr))
}

pr <- rep(0.5, N)
check("horvitz_thompson()",
      horvitz_thompson(y ~ z, data = d, condition_prs = pr)$coefficients[[1]],
      ht_estimate(d$y, d$z, pr))
#>       gap holds
#> 1 0.0e+00  TRUE

Variance

The variance estimator is the conservative bound of Aronow and Middleton (2013), built from Young’s inequality. In its general form,

$$\widehat{\mathbb{V}}[\widehat{\tau}] = \frac{1}{N^2}\left[ \sum_{i: z_i = 1}\widetilde{Y}_{1i}^2 + \sum_{i: z_i = 0}\widetilde{Y}_{0i}^2 + \sum_{i \neq j} A_{ij}\,\widetilde{Y}_i\widetilde{Y}_j \right]$$

where the cross terms enter with a minus sign when $i$ and $j$ are in opposite conditions, and

$$A_{ij} = 1 - \frac{\pi_i\pi_j}{\pi_{ij}}.$$

Everything below is that expression with $A_{ij}$ worked out for a particular design.

Simple (Bernoulli) randomization. Assignments are independent, so $\pi_{ij} = \pi_i\pi_j$, every $A_{ij}$ is zero, and the bound collapses to

$$\widehat{\mathbb{V}}[\widehat{\tau}] = \frac{1}{N^2}\left[\sum_{i: z_i = 1}\widetilde{Y}_{1i}^2 + \sum_{i: z_i = 0}\widetilde{Y}_{0i}^2\right].$$

Which is short enough to check directly, and it is the variance the estimate above was reported with, since a bare probability vector says nothing about dependence between units.

Y1 <- d$y[d$z == 1] / 0.5
Y0 <- d$y[d$z == 0] / 0.5

check("horvitz_thompson() variance, simple randomization",
      horvitz_thompson(y ~ z, data = d, condition_prs = pr)$std.error[[1]],
      sqrt((sum(Y1^2) + sum(Y0^2)) / N^2))
#>       gap holds
#> 1 0.0e+00  TRUE

Complete randomization. With $n$ units of which $m_1$ go to condition 1 and $m_0$ to condition 0, exchangeability gives the joint probabilities in closed form, $\pi_{11} = m_1(m_1-1)/(n(n-1))$ and so on, so $A_{ij}$ takes only three values:

$$A^{11} = 1 - \frac{m_1(n-1)}{n(m_1-1)}, \qquad A^{00} = 1 - \frac{m_0(n-1)}{n(m_0-1)}, \qquad A^{10} = \frac{1}{n}.$$

The cross coefficient collapses to $1/n$ for any complete design. Because the three coefficients are constant within pair type, the double sum needs no matrix: $\sum_{i \neq j}\widetilde{Y}_{1i}\widetilde{Y}_{1j}$ is $\left(\sum_i \widetilde{Y}_{1i}\right)^2 - \sum_i \widetilde{Y}_{1i}^2$. The whole variance is therefore four sums over the data (the total and the sum of squares of the weighted outcomes, in each condition) plus the design’s $n$ and $m_1$. Where 1.x built an $N \times N$ matrix of joint probabilities, 2.0 evaluates a scalar formula.

When the design implies a non-integer $m_1 = \pi_1 n$, the realized count is $\lfloor m_1 \rfloor$ or $\lfloor m_1 \rfloor + 1$, and the joint probabilities average over that mixture.

Blocked. Randomization is complete and independent within each block, so the contributions add:

$$\widehat{\mathbb{V}}[\widehat{\tau}] = \frac{1}{N^2}\sum_{j=1}^{J} C_j$$

with $C_j$ the complete-randomization expression evaluated on block $j$’s units, at that block’s $N_j$ and $m_{1j}$.

Clustered. Assignment is at the cluster level, so the weighted outcomes are summed within cluster first, and the same expression is applied to the $S$ cluster totals: complete randomization at the cluster level if the clusters were completely randomized, the simple form if they were not. Blocked and clustered designs aggregate within cluster and then sum over blocks.

Arbitrary designs. Given a permutation matrix, the joint probabilities come from one tcrossprod() and $A_{ij}$ is evaluated directly, at $O(n^2)$. A pair of units that can never appear together in the observed conditions has $\pi_{ij} = 0$, and its term is not identified at any sample size. Those terms are dropped and replaced by the Young’s inequality bound: the unidentified quantity is at most $(y_i^2 + y_j^2)/2$ within a condition and $y_i^2 + y_j^2$ across conditions, and $y_i^2$ is estimated from the single observation of it. Left as $1 - x/0$, as in an earlier implementation, the whole variance became $-\infty$ and then a silent NA.

What the declaration buys

condition_prs takes an ra_declaration from randomizr, a named vector of marginal probabilities, or a matrix of per-unit probabilities. The choice is visible at the call site and it determines which variance you get.

A declaration carries the block structure, the cluster structure, the per-unit marginals, and whether the randomization was simple or complete, which is exactly what the design-aware expressions above need. A bare probability vector carries only the marginals, so estimatr falls back to the simple-randomization bound, which is valid for any design and exact only for Bernoulli assignment. For a complete or blocked design it overstates the uncertainty.

In 1.x the same distinction existed but was buried in which combination of five arguments happened to be supplied.

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

# The design-aware variance
horvitz_thompson(Y ~ Z, data = dat_ht, condition_prs = decl)$std.error
#>      1 
#> 0.1414

# The conservative bound, from the marginals alone
horvitz_thompson(Y ~ Z, data = dat_ht,
                 condition_prs = c("0" = 0.6, "1" = 0.4))$std.error
#>      1 
#> 0.1506

Confidence intervals and testing

Inference for the Horvitz-Thompson estimator rests on a normal approximation:

$$\mathrm{CI}^{1-\alpha} = \left(\widehat{\tau} + z_{\alpha/2}\sqrt{\widehat{\mathbb{V}}[\widehat{\tau}]},\; \widehat{\tau} + z_{1-\alpha/2}\sqrt{\widehat{\mathbb{V}}[\widehat{\tau}]}\right)$$

with two-sided p-values from the same distribution.

Every promise in one table

Promise	Largest relative gap	Holds
lm_robust()	4.6e-16	TRUE
lm_robust(weights = )	2.2e-16	TRUE
lm_robust(se_type = ‘classical’)	2.8e-17	TRUE
lm_robust(se_type = ‘HC0’)	6.9e-17	TRUE
lm_robust(se_type = ‘HC1’)	1.1e-16	TRUE
lm_robust(se_type = ‘HC2’)	5.6e-17	TRUE
lm_robust(se_type = ‘HC3’)	8.3e-17	TRUE
lm_robust(clusters = )	1.1e-16	TRUE
lm_robust(se_type = ‘stata’)	1.1e-16	TRUE
lm_robust(fixed_effects = )	7.8e-16	TRUE
lm_lin()	0.0e+00	TRUE
iv_robust()	8.5e-16	TRUE
lh_robust()	0.0e+00	TRUE
difference_in_means()	4.9e-16	TRUE
horvitz_thompson()	0.0e+00	TRUE
horvitz_thompson() variance, simple randomization	0.0e+00	TRUE

stopifnot(all(unlist(CHECKS) < 1e-10))

Every promise above is met. The numbers are computed when the vignette is built, so they are what your installed copy produces rather than values recorded from a run somewhere else.

That last line is stopifnot() rather than a printed TRUE on purpose. A vignette that computes its own table can report FALSE in a cell and still build, which would leave a broken promise sitting inside a clean R CMD check. Written this way the document refuses to build, so the check fails and the table cannot quietly disagree with the sentences above it. The margin is wide enough for that to be safe: the gaps sit at 1e-15 or below against a tolerance of 1e-10, so the linear algebra library on your machine would have to be five orders of magnitude worse than the one this was written on before the build broke.

What this does not cover

The checks above are a demonstration, not a proof, and they are deliberately a small set.

They say nothing about what these estimators are good for. That difference_in_means() computes the difference in means is a fact about this package. Whether that quantity is unbiased for your estimand, whether its interval covers, whether covariate adjustment helps you: none of that is estimatr’s to guarantee, and none of it is checked here. The papers cited throughout are where those questions are answered. The guarantee is implementation, and the demonstration is arithmetic.

Each promise is checked in one configuration. The test suite checks many: the same identities across weighted and unweighted fits, single and multivariate outcomes, one and two absorbed factors, instrumental variables with and without clusters, and the rank-deficient and near-saturated designs that have caused bugs. That is where a guarantee is enforced. What this document adds is that the promises are stated in words a reader can disagree with, next to the mathematics they are supposed to implement, and measured where a reader can watch.

Agreement with a formula written here is not agreement with the literature. Transcribing HC2 into this document and matching it shows estimatr computes what it says. It cannot show that the definition is the one the field settled on, which is why every definition above carries its citation, and why the suite compares against sandwich, clubSandwich, ivreg, Stata’s regress, areg and ivregress, fixest, plm and blkvar, none of which shares any lineage with this package. Two known divergences are pinned from both sides there rather than dropped: weighted HC2 and HC3 differ from Stata by a bounded amount, and iv_robust() uses second-stage leverage, which agrees with ivreg exactly and departs from AER::ivreg()’s deprecated hatvalues() method by up to 18.5 percent.

The data here are well conditioned. Every fit above is full rank with far more observations than parameters. Leverage exactly equal to one is benign; leverage marginally above one is not, and HC2 and HC3 are guarded there rather than answered, warning and contributing zero for the offending rows. A cluster-robust variance on a single cluster is refused outright, where 1.0.6 returned a standard error of 5.9e-17 and a zero-width interval in silence. Those are documented in NEWS.md, and none of them is visible in a table built on rnorm().

References

Aronow, Peter M., and Joel A. Middleton. 2013. “A Class of Unbiased Estimators of the Average Treatment Effect in Randomized Experiments.” Journal of Causal Inference 1 (1): 135–54. https://doi.org/10.1515/jci-2012-0009.

Aronow, Peter M., and Cyrus Samii. 2017. “Estimating Average Causal Effects Under General Interference, with Application to a Social Network Experiment.” The Annals of Applied Statistics 11 (4). https://doi.org/10.1214/16-AOAS1005.

Bell, Robert M., and Daniel F. McCaffrey. 2002. “Bias Reduction in Standard Errors for Linear Regression with Multi-Stage Samples.” Survey Methodology 28 (2): 169–82.

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

Freedman, David A. 2008. “On Regression Adjustments in Experiments with Several Treatments.” The Annals of Applied Statistics 2 (1): 176–96. https://doi.org/10.1214/07-AOAS143.

Gerber, Alan S., and Donald P. Green. 2012. Field Experiments: Design, Analysis, and Interpretation. New York: W.W. Norton.

Imai, Kosuke, Gary King, and Clayton Nall. 2009. “The Essential Role of Pair Matching in Cluster-Randomized Experiments, with Application to the Mexican Universal Health Insurance Evaluation.” Statistical Science 24 (1): 29–53. https://doi.org/10.1214/08-STS274.

Imbens, Guido W., and Michal Kolesár. 2016. “Robust Standard Errors in Small Samples: Some Practical Advice.” Review of Economics and Statistics 98 (4): 701–12. https://doi.org/10.1162/REST_a_00552.

Li, Jianzhu, and Richard Valliant. 2009. “Survey Weighted Hat Matrix and Leverages.” Survey Methodology 35 (1): 15–24.

Lin, Winston. 2013. “Agnostic Notes on Regression Adjustments to Experimental Data: Reexamining Freedman’s Critique.” The Annals of Applied Statistics 7 (1): 295–318. https://doi.org/10.1214/12-AOAS583.

Long, J. Scott, and Laurie H. Ervin. 2000. “Using Heteroscedasticity Consistent Standard Errors in the Linear Regression Model.” The American Statistician 54 (3): 217–24. https://doi.org/10.1080/00031305.2000.10474549.

MacKinnon, James G., and Halbert White. 1985. “Some Heteroskedasticity-Consistent Covariance Matrix Estimators with Improved Finite Sample Properties.” Journal of Econometrics 29 (3): 305–25. https://doi.org/10.1016/0304-4076(85)90158-7.

Middleton, Joel A., and Peter M. Aronow. 2015. “Unbiased Estimation of the Average Treatment Effect in Cluster-Randomized Experiments.” Statistics, Politics and Policy 6 (1-2): 39–75. https://doi.org/10.1515/spp-2013-0002.

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

Pustejovsky, James E., and Elizabeth Tipton. 2018. “Small-Sample Methods for Cluster-Robust Variance Estimation and Hypothesis Testing in Fixed Effects Models.” Journal of Business & Economic Statistics 36 (4): 672–83. https://doi.org/10.1080/07350015.2016.1247004.

Romano, Joseph P., and Michael Wolf. 2017. “Resurrecting Weighted Least Squares.” Journal of Econometrics 197 (1): 1–19. https://doi.org/10.1016/j.jeconom.2016.10.003.

Samii, Cyrus, and Peter M. Aronow. 2012. “On Equivalencies Between Design-Based and Regression-Based Variance Estimators for Randomized Experiments.” Statistics & Probability Letters 82 (2): 365–70. https://doi.org/10.1016/j.spl.2011.10.024.