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 7.8e-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 1.2e-15  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.

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 6.9e-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 5.6e-17  TRUE
#> 2     HC1 9.7e-17  TRUE
#> 3     HC2 8.3e-17  TRUE
#> 4     HC3 6.9e-17  TRUE

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.

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.2e-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 3.4e-16  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.

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 1.3e-15  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.

Correspondence with Stata

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

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

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 et al. (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.

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 2.5e-16  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.