Unbiased Regression-Adjusted Estimation of Average Treatment Effects in Randomized Controlled Trials

Unbiased Regression-Adjusted Estimation of Average Treatment Effects in Randomized Controlled Trials

Introduction

Establishing causal relationships is a central goal in economics and the social sciences. In observational studies, selection bias often obscures the effects of interventions. Randomized controlled trials (RCTs) address this problem by eliminating systematic pre-randomization differences between treatment and control groups. For this reason, RCTs occupy a central place in modern empirical research on causal inference.

As experimental practice has matured, attention has shifted toward improving efficiency and statistical power. Researchers now collect extensive pretreatment information to predict outcomes and use regression adjustment to explain part of the outcome variance. This approach reduces noise and increases the precision of treatment effect estimation without requiring complex assignment mechanisms.

Most formal guarantees for regression-adjusted estimators rely on asymptotic arguments. However, early-stage clinical trials, marketing geo-experiments, and studies of firm- or country-level policies frequently operate in small sample regimes, where standard regression-adjusted estimators become biased, with undesirable consequences for treatment effect estimation and policy design freedman2008regression. Moreover, in small samples, high-leverage observations (that is, observations whose covariate values exert a disproportionate influence on the regression fit) further degrade the performance of regression adjustment young2019channeling.

To address these challenges, we propose a regression adjustment framework that is unbiased in finite samples and robust to influential observations. A key component of our analysis is a leave-one-out regression adjustment (LOORA) procedure, which removes the bias of classical regression adjustment. Within this framework, we develop two estimators of the average treatment effect (ATE): LOORA-HT, a Horvitz–Thompson estimator for simple random assignment; and LOORA-DM, a difference-in-means estimator for complete random assignment. LOORA-HT and LOORA-DM are easy to implement and deliver a performance that previously required complex assignment mechanisms harshaw2022design.

LOORA-HT and LOORA-DM simultaneously address three core problems in regression-adjusted estimation of treatment effects. First, they answer the critique of freedman_2008_b by providing unbiased, regression-adjusted estimators for RCTs. Second, they attain the asymptotic variance of the estimator in lin2013agnostic, which is the efficient variance among all linearly adjusted estimators under simple and complete random assignment. Third, by stabilizing estimation in the presence of influential observations, they reduce the sensitivity to leverage documented in young2019channeling.

Related Work

The practice of regression adjustment in randomized experiments dates at least to the use of analysis of covariance (ANCOVA) in classical experimental design, where a linear adjustment removes the effects of pretreatment covariates to improve precision.\footnote{See, for example, discussions in fisher1971design.} Despite its pervasiveness, the formal large-sample properties of regression adjustment after randomization have been repeatedly questioned. In two influential critiques, freedman2008regression,freedman_2008_b argue that ordinary least squares (OLS) adjustments can introduce bias when the regression model is misspecified, and can be less efficient than the unadjusted difference-in-means estimator, even asymptotically.

In response, lin2013agnostic shows that, under complete random assignment, an OLS regression of outcomes on treatment, covariates, and all treatment-by-covariate interactions delivers an estimator that is consistent for the average treatment effect and asymptotically efficient among the class of linearly adjusted estimators. While these results clarify the large sample properties of regression adjustment after randomization, two important gaps remain: (i) lin2013agnostic does not provide finite-sample unbiasedness guarantees for regression adjustment, and (ii) the interacted specification in lin2013agnostic can be unstable in moderate samples because it increases the number of regressors by adding interaction terms and potentially inflating leverage scores and sensitivity to influential observations.

This article proposes two estimators, LOORA-HT under simple random assignment and LOORA-DM under complete random assignment, that (i) are exactly unbiased for the finite-population average treatment effect, (ii) admit closed-form variance expressions in finite samples, and (iii) are asymptotically efficient among linearly adjusted estimators, attaining the variance bound of lin2013agnostic. We restrict our attention to linear regression adjustments because of their ubiquity in empirical research and the potential instability of nonlinear adjustments in small samples. Linear adjustments, however, can accommodate nonlinear relationships through covariate transformations.

LOORA combines leave-one-out regression adjustment and ridge regression to remove the finite-sample bias of regression-adjusted estimators and to curb the influence of high-leverage observations. For the difference-in-means setting, we show that LOORA-DM admits a representation that is equivalent to a leave-two-out construction in the sense of spiess2025optimal, thereby satisfying necessary conditions for unbiasedness under complete random assignment. In this way, our estimators provide constructive finite-sample analogues of the asymptotic guarantees in lin2013agnostic, while remaining valid without the assumption of correctly-specified parametric outcome models.

In two within-subject experiments that provide a realistic ground truth for the joint distribution of potential outcomes, LOORA-HT and LOORA-DM substantially reduce estimation bias and yield confidence intervals with close-to-nominal coverage.

A related line of work examines precision gains in randomized studies through experimental design rather than regression adjustment. harshaw2022design proposes a Gram--Schmidt walk design, which chooses treatment assignments in a way that approximately balances covariates. Their goal is to optimize the randomization design to deliver low-variance Horvitz--Thompson estimators in finite samples. In contrast, we consider simple or complete random assignment and instead adjust outcomes through leave-one-out regressions. This distinction matters in practice, because our methods can be applied to commonly used experimental designs without the need to change the assignment rule.

ghadiri2023finite introduces a cross-fitted, regression-adjusted Horvitz--Thompson estimator under simple random assignment. It splits the sample into two folds, estimates regression coefficients on one fold, and applies those coefficients to adjust outcomes in the other. This construction yields an exactly unbiased estimator with an explicit nonasymptotic variance bound. However, the bound is looser than the lin2013agnostic asymptotic efficiency bound, and the article does not provide a procedure for constructing confidence intervals. lei2021regression proposes a bias-correction procedure for regression adjustment when the number of covariates grows with the sample size. The method removes higher-order bias terms but does not yield exact unbiasedness in finite samples. chang2024exact develops an exact bias-correction approach that eliminates the bias, but does not derive the finite-sample variance of the resulting estimator. gu2025assumption studies regression adjustment with high-dimensional covariates under a superpopulation framework and develops consistent estimators in this setting. zhao2024covariate considers a related framework and proposes a debiasing approach based on higher-order influence functions that corrects own-observation bias arising from regression adjustment. Their theoretical development primarily focuses on unregularized linear adjustment, or regimes in which the ridge regularization parameter vanishes asymptotically, which is equivalent to using the Moore--Penrose pseudoinverse of the Gram matrix in the associated regression problem. In contrast, our analysis accommodates arbitrary regularization parameters. This flexibility enables explicit damping of high-leverage observations, which can otherwise induce unstable or unreliable inference, especially in finite-population or small-sample settings.

A large concurrent literature studies regression adjustment and inference in more structured designs: stratified or blocked experiments, matched pairs, cluster randomization, or covariate-adaptive randomization. For example, bai2024primer survey modern design-based analysis of randomized experiments, emphasizing the role of stratification, regression adjustment and cluster randomized experiments. bai2025new,bai2025inference develop design-based variance estimators and inference procedures for finely stratified and matched-pair experiments, including settings with imperfect compliance. cytrynbaum2024covariate characterizes asymptotically optimal linear covariate adjustment for general stratified randomization schemes.

Our approach provides an appealing alternative in settings where stratification becomes impractical or infeasible. In particular, stratification requires covariate information at the design stage, which is often unavailable in online experiments and field experiments with asynchronous enrollment. By contrast, simple and complete randomization designs apply naturally to these settings.

(ref) provides a more detailed comparison of our estimators with related methods.

Notation

\paragraph{Vectors and matrices.} We denote matrices and vectors by bold uppercase and lowercase letters, respectively. The $i^{\text{th}}$ entry of a vector ${\bm{u}}$ is denoted by $u_i$. The transposed $i^{\text{th}}$ row of a matrix ${\bm{X}}$ is denoted by ${\bm{x}}_i$, and its $(i,j)^{\text{th}}$ entry by $x_{ij}$. For a constant $c$ and a vector ${\bm{u}}$, expressions such as $c + {\bm{u}}$ and ${\bm{u}}^{-1}$ are interpreted entrywise. When ${\bm{u}}$ and ${\bm{v}}$ are vectors of equal dimension, ${\bm{u}} {\bm{v}}$ and ${\bm{u}} / {\bm{v}}$ denote entrywise (Hadamard) product and division, respectively. A vector’s associated diagonal matrix appears in uppercase; for instance, $\bT$ denotes the diagonal matrix ${\bm{t}}$ as its main diagonal. For any vector ${\bm{v}} \in \R^n$, the notation ${\bm{v}}_{-i}$ refers to the vector in $\R^{n-1}$ obtained by removing the $i^{\text{th}}$ entry of ${\bm{v}}$. Likewise, ${\bm{X}}_{-i}$ denotes the matrix obtained by deleting the $i^{\text{th}}$ row of ${\bm{X}}$. We denote by $\boldsymbol{1}$ the vector with all entries equal to one. For a matrix ${\bm{X}} \in \R^{n\times k}$, let $\overline{{\bm{X}}}$ denote the matrix whose rows are all equal to the average row of ${\bm{X}}$, that is, each row of $\overline{{\bm{X}}}$ equals $ \boldsymbol{1}^\top {\bm{X}}/n$. Similarly, for a vector ${\bm{y}} \in \R^n$, let $\overline{{\bm{y}}}$ denote the vector whose entries are all equal to the average of the entries in ${\bm{y}}$.

{\em Norms and asymptotic notation.} For a vector ${\bm{u}} \in \R^n$, $\left\|{\bm{u}}\right\|_{2} = \sqrt{{\bm{u}}^\top {\bm{u}}}$ denotes the Euclidean norm, and $\left\|{\bm{u}}\right\|_{\infty} = \max_{i \in [n]} |u_i|$ denotes the $\ell_\infty$ norm. The operator $\tr(\cdot)$ denotes the trace of a square matrix. For a matrix ${\bm{X}} \in \R^{n \times k}$, $\left\|{\bm{X}}\right\|_{F} = \sqrt{\tr({\bm{X}}^\top {\bm{X}})}$ denotes the Frobenius norm, and $\left\|{\bm{X}}\right\|_{2,\infty} = \max_{i \in [n]} \left\|{\bm{x}}_i\right\|_{2}$ denotes its $(2,\infty)$-operator norm. Let $\mathbb{N}$ be the set of natural numbers. For any two functions $f, g \colon \mathbb{N} \to \R$, we write $f(n) = O(g(n))$ if there exist constants $C > 0$ and $n_0 \in \mathbb{N}$ such that for all $n \ge n_0$, $\lvert f(n) \rvert \le C \lvert g(n) \rvert$. Similarly, we write $f(n) = \Omega(g(n))$ if there exist constants $C > 0$ and $n_0 \in \mathbb{N}$ such that for all $n \ge n_0$, $\lvert f(n) \rvert \ge C \lvert g(n) \rvert$.

{\em Projection matrices.} Consider a full-rank matrix ${\bm{X}} \in \R^{n \times k}$ with $n \ge k$. We denote the projection matrix by $\bH = {\bm{X}}({\bm{X}}^\top {\bm{X}})^{-1}{\bm{X}}^\top$, and its diagonal entries by $h_{ii} = {\bm{x}}_i^\top ({\bm{X}}^\top {\bm{X}})^{-1} {\bm{x}}_i$, which represent the leverage scores of rows $i \in [n]$. The identity matrix is denoted by ${\bm{I}}$. For $\lambda \ge 0$, the ridge projection matrix is $\bH_\lambda = {\bm{X}}({\bm{X}}^\top {\bm{X}} + \lambda {\bm{I}})^{-1}{\bm{X}}^\top$, whose diagonal entries $h_{\lambda ii} = {\bm{x}}_i^\top ({\bm{X}}^\top {\bm{X}} + \lambda {\bm{I}})^{-1} {\bm{x}}_i$ are the ridge leverage scores.

Finite Population Framework

We operate within the Neyman--Rubin potential outcomes framework neyman1923application,rubin1974estimating for a randomized controlled trial over a set of units indexed by $[n] = \{1, \ldots, n\}$. Each unit is assigned to either treatment or control. The vector ${\bm{d}} \in \R^{n}$ denotes the treatment assignment, with $d_i = 1$ if unit $i$ is in treatment and $d_i = 0$ if it is in control. We also define ${\bm{z}} = 2({\bm{d}} - 1/2)$, so that $z_i = 1$ for treated units and $z_i = -1$ otherwise. Finally, we define $\bq$ as the vector with $i^{\text{th}}$ element equal to $q_i = p_i d_i + (1-p_i)(1-d_i)$, where $p_i$ is the treatment-assignment probability for unit $i$.

For each unit $i$, the potential outcomes under treatment and control are denoted by $y^{(1)}_i$ and $y^{(0)}_i$, respectively. We treat potential outcomes as fixed quantities; the only source of randomness in our setting is the treatment assignment. The observed outcome for unit $i$ is \[ y_i = y^{(1)}_i d_i + y^{(0)}_i (1 - d_i). \] We collect the observed and potential outcomes across all units in the vectors ${\bm{y}} = (y_1, \ldots, y_n)^\top$, ${\bm{y}}^{(1)} = (y^{(1)}_1, \ldots, y^{(1)}_n)^\top$, and ${\bm{y}}^{(0)} = (y^{(0)}_1, \ldots, y^{(0)}_n)^\top$.

The average treatment effect (ATE) is defined as \[ \tau = \frac{1}{n} \sum_{i=1}^n \bigl(y^{(1)}_i - y^{(0)}_i\bigr). \]

Each unit $i \in [n]$ comes with a vector of fixed pretreatment covariates ${\bm{x}}_i \in \R^k$. The matrix ${\bm{X}} = ({\bm{x}}_1, \ldots, {\bm{x}}_n)^\top \in \R^{n \times k}$ collects covariate values for all units. For simplicity of exposition, we assume throughout the article that ${\bm{X}}$ has full column rank. However, since our estimators are based on ridge regression, the results extend directly to cases in which ${\bm{X}}$ is rank-deficient.

We consider two standard treatment assignment mechanisms. Under simple random assignment, each unit is independently assigned to treatment with probability $p_i \in (0,1)$, with $\bp = (p_1, \ldots, p_n)^\top$. Under complete random assignment, the number of treated units is fixed at $n_T \in [n - 1]$, and the number of control units is $n_C = n - n_T$. Treatment is assigned by selecting uniformly at random from all possible subsets of $n_T$ units in the sample.

Under simple random assignment, the Horvitz--Thompson (HT) estimator horvitz1952generalization is defined as \[ \widehat{\tau}_{\text{HT}} = \frac{1}{n} \sum_{i=1}^n \frac{d_i y_i}{p_i} \;-\; \frac{1}{n} \sum_{i=1}^n \frac{(1 - d_i) y_i}{1 - p_i}. \]

Under complete random assignment, the difference-in-means (DM) estimator is given by \[ \widehat{\tau}_{\text{DM}} = \frac{1}{n_T} \sum_{i=1}^n d_i y_i \;-\; \frac{1}{n_C} \sum_{i=1}^n (1 - d_i) y_i. \]

Regression Adjustment for Horvitz-Thompson Estimator

In this section, we first review some known results on the variance of the Horvitz--Thompson (HT) estimator. We then introduce LOORA-HT, our leave-one-out regression-adjusted version of the HT estimator ((ref)). We establish its unbiasedness and derive an exact expression for its variance in the finite-population setting ((ref)). Next, we analyze the asymptotic behavior of LOORA-HT ((ref)) and prove its asymptotic efficiency ((ref)). Finally, we describe our approach to variance estimation and confidence interval construction for LOORA-HT ((ref)).

It is well known that, under simple random assignment, the classical Horvitz--Thompson (HT) estimator is unbiased for $\tau$. The variance of the HT estimator under simple random assignment is given by

where ${\boldsymbol \mu} = (\mu_1, \ldots, \mu_n)^\top$, and

Let $\breve{y}_i = y_i - {\bm{x}}_i^\top {\bm{b}}$ denote the covariate-adjusted outcomes, where ${\bm{b}}$ is a fixed vector of coefficients. The corresponding potential outcomes are

Because $\breve{y}^{(1)}_i - \breve{y}^{(0)}_i = y^{(1)}_i - y^{(0)}_i$, it follows that the HT estimator applied to the adjusted outcomes $\breve{y}_i$ is also unbiased for $\tau$. Since ${\bm{b}}$ is fixed, applying (ref) to the adjusted potential outcomes implies that the variance of the covariate-adjusted HT estimator is

where ${\bm{R}}$ is a diagonal matrix associated with the vector ${\bm{r}} \in \R^{n}$ defined by $r_i = \sqrt{p_i (1 - p_i)}$. Let $\widetilde{{\bm{X}}} = {\bm{R}}^{-1} {\bm{X}}$. Then the variance in (ref) is minimized by the choice of coefficients

Equation (ref) implies that regression adjustment can reduce the variance of the HT estimator. In practice, however, only one potential outcome is observed for each unit, so ${\boldsymbol \mu}$ cannot be computed directly, and the minimization problem in (ref) cannot be solved exactly. We address this challenge in LOORA-HT by replacing ${\boldsymbol \mu}$ with a proxy vector that is equal to ${\boldsymbol \mu}$ in expectation.

Leave-One-Out Regression-Adjusted Horvitz-Thompson

(ref) describes LOORA-HT. Rather than using a single value for ${\bm{b}}$, LOORA-HT applies a different vector of regression coefficients, $\widehat{\bbeta}_{\lambda}^{(-i)}$, to each unit $i$ in the sample. LOORA-HT constructs its adjustment vectors using a proxy vector $\widetilde{{\bm{y}}}$ that coincides with ${\boldsymbol \mu}$ in expectation, \[ \widetilde y_i =

\] LOORA-HT is defined as \[ \widehat{\tau}_{\textup{LHT}} = \frac{1}{n}\sum_{i=1}^n \frac{z_i}{q_i} (y_i - {\bm{x}}_{i}^\top \widehat{\bbeta}_{\lambda}^{(-i)}), \] where $\widehat{\bbeta}_{\lambda}^{(-i)}$ is obtained from a regression on the leave-one-out data, $(\widetilde{{\bm{y}}}_{-i}, \widetilde{{\bm{X}}}_{-i})$. In general, leave-one-out coefficients can be unstable, particularly for rows with high leverage scores (see (ref)). To address this challenge, we use ridge regression to mitigate the impact of high leverage scores. In (ref) below, we provide an upper bound on the ridge leverage scores as a function of the ridge regularization parameter.

The Variance of LOORA-HT

Let the ridge projection matrix of $\widetilde{{\bm{X}}}$ be \[ \widetilde{\bH}_\lambda = \widetilde{{\bm{X}}} (\widetilde{{\bm{X}}}^\top \widetilde{{\bm{X}}} + \lambda {\bm{I}})^{-1} \widetilde{{\bm{X}}}^\top. \] We denote its $(i,j)^{\text{th}}$ element by $\widetilde{h}_{\lambda ij}$ and its $i^{\text{th}}$ diagonal element by $\widetilde{h}_{\lambda ii}$. Next, we define the vector \[ {\bm{t}} = \frac{(1 - \bp)^2 {\bm{y}}^{(1)} - \bp^2 {\bm{y}}^{(0)}}{{\bm{r}}}, \] which characterizes how far $\widetilde{{\bm{y}}}$ deviates from ${\boldsymbol \mu}$, \[ \widetilde{{\bm{y}}} - {\boldsymbol \mu} = \frac{{\bm{z}} {\bm{t}}}{\bq}. \] The ridge regression coefficient on $({\boldsymbol \mu}, \widetilde{{\bm{X}}})$ is

The next theorem establishes the unbiasedness of LOORA-HT and provides an exact expression for its variance.

The first term in (ref) mirrors the variance of the infeasible regression-adjusted version HT estimator in (ref), and quantifies the variance reduction achieved by regression adjustment. The factor $(1 - \widetilde{h}_{\lambda ii})^{-2}$ arises from the removal of a single row in the leave-one-out regressions. More specifically, by (ref), the first term of (ref) can be written as \[ \frac{1}{n^2}\sum_{i=1}^n (\widetilde{{\bm{x}}}_i^\top \bbeta_\lambda^{(-i)} - \mu_i)^2, \] where \[ \bbeta_{\lambda}^{(-i)} = \mathop{\mathrm{argmin}\,}\limits_{{\bm{b}}\in \mathbb R^k} \left\|{\boldsymbol \mu}_{-i}-\widetilde{{\bm{X}}}_{-i} {\bm{b}}\right\|_{2}^2 + \lambda\left\|{\bm{b}}\right\|_{2}^2. \]

When a leverage score $\widetilde{h}_{\lambda ii}$ approaches one, the variance becomes unbounded. In such high-leverage cases, it is essential to employ ridge regression with a sufficiently large regularization parameter $\lambda$. The following lemma provides an upper bound on the ridge leverage scores.

High-leverage observations distort regression-adjusted estimates. LOORA applies ridge regularization, which stabilizes the adjustment and delivers more reliable estimates.

The second term of (ref) arises from using the random vector $\widetilde{{\bm{y}}}$, with $\E[\widetilde{{\bm{y}}}]={\boldsymbol \mu}$, in place of ${\boldsymbol \mu}$ in the regression adjustment. In other words, it captures the error introduced by estimating $\bbeta_{\lambda}^{(-i)}$ with $\widehat{\bbeta}_{\lambda}^{(-i)}$.

We next provide a loose upper bound for this term to show that it scales as $k / n^2$. Let \[ \widetilde{{\bm{X}}}_{\lambda} =

. \] It follows that

where the second inequality holds because $\widetilde{\bH}_\lambda$ is an $(n \times n)$ block of the projection matrix associated with $\widetilde{{\bm{X}}}_{\lambda}$, and the final equality holds since the rank of $\widetilde{{\bm{X}}}_{\lambda}$ is $k$, and the trace of a projection matrix equals its rank.

Therefore, the second term of (ref) is bounded by

where $\widetilde{\bh}_\lambda$ denotes the vector of diagonal entries of $\widetilde{\bH}_{\lambda}$. That is, for a fixed number of covariates, the second term of (ref) scales as $1 / n^2$ as $n$ increases, and is typically much smaller than the first term, which scales as $1 / n$.

In particular, for $\lambda = \left\|\widetilde{{\bm{X}}}\right\|_{2,\infty}^2$, by (ref), the variance of LOORA-HT is bounded by \[ \frac{4}{n^2}\sum_{i=1}^n (\widetilde{{\bm{x}}}_{i}^\top \bbeta_\lambda - \mu_i)^2 + \frac{8k}{n^2} \left\|{\bm{t}}/{\bm{r}}\right\|_{\infty}^2\,. \]

We next study the asymptotic behavior of the LOORA-HT estimator.

Asymptotic Normality of LOORA-HT

This section derives a large-sample approximation to the distribution of LOORA-HT under the following assumptions.

\Crefname{enumi}{Assumption}{Assumptions}

Assumptions (ref) and (ref) provide bounds on potential outcomes and covariates, but allow the means of $(y_i^{(1)}, y_i^{(0)}, x_i)$ to fluctuate, without necessarily converging to a limit as $n$ increases. Assumption (ref) requires that each experimental unit be assigned to treatment and control with strictly positive probabilities. It rules out degenerate experimental designs in which causal effects cannot be identified.

Let \[ \bbeta_{\lambda} = (\widetilde{{\bm{X}}}^\top \widetilde{{\bm{X}}} + \lambda {\bm{I}})^{-1} \widetilde{{\bm{X}}}^\top {\boldsymbol \mu} \] be the true solution of the ridge regression problem on ${\boldsymbol \mu}$. The next theorem establishes the large-sample distribution of the LOORA-HT estimator and identifies its asymptotic variance in closed form.

In the next subsection, we build on Theorem (ref) to analyze the asymptotic efficiency of LOORA-HT.

Asymptotic Efficiency of LOORA-HT

We now examine the asymptotic efficiency of the LOORA-HT estimator. Theorem (ref) established the asymptotic normality of LOORA-HT. In this subsection, we show that under simple random assignment, LOORA-HT is asymptotically efficient among the class of Horvitz-Thompson linearly adjusted estimators, defined as estimators of the form, \[ {\widehat{\tau}}_{\textup{HLA}} = \frac{1}{n} \sum_{i=1}^n \frac{z_i}{q_i} \bigl(y_i - {\bm{x}}_i^\top \widehat{\bgamma}_{n}^{(-i)}\bigr), \] where there exists a bounded deterministic sequence of vectors $\bgamma_1,\bgamma_2,\ldots$ and a sequence of vectors $\widehat{\bgamma}_1,\widehat{\bgamma}_2,\ldots$ such that

The class of HT linearly adjusted estimators is larger than the class of estimators usually considered for the purpose of determining efficiency (e.g., see cytrynbaum2024covariate), because it allows for unit-dependent adjustment vectors, $\widehat{\bgamma}_{n}^{(-i)}$.

LOORA-HT is contained in the class of HT linearly adjusted estimators since

where \[ \widehat{\bbeta}_{\lambda} = (\widetilde{{\bm{X}}}^\top \widetilde{{\bm{X}}} + \lambda {\bm{I}})^{-1} \widetilde{{\bm{X}}}^\top \widetilde{{\bm{y}}}, \] as shown in the proof of (ref). The next result describes the asymptotic behavior of HT linearly adjusted estimators.

Using this result, we can establish the asymptotic efficiency of LOORA-HT in the case where $\left\|\widetilde{{\bm{X}}} \bbeta_{0} - {\boldsymbol \mu}\right\|_{2} = \Omega(\sqrt{n})$.

For the case of $\lambda=0$ and provided that the regression adjustment includes an intercept, by (ref), the variance of LOORA-HT reduces to \[ \min_{\bbeta \in \R^k} \frac{1}{n^2} \left\|(\widetilde{{\bm{X}}} - \overline{\widetilde{{\bm{X}}}})\bbeta - ({\boldsymbol \mu} - \overline{{\boldsymbol \mu}})\right\|_{2}^2, \] This expression is similar to the variance of the estimator of lin2013agnostic.

In the next subsection, we develop valid and asymptotically correct confidence intervals for LOORA-HT.

Confidence Intervals for LOORA-HT

The next theorem presents a variance estimator for $\widehat{\tau}_{\textup{LHT}}$.

As is generally the case in all finite-population inference abadie2020sampling, $\widehat{V}_{\textup{LHT}}$ is consistent under homogeneous individual treatment effects and conservative when treatment effects are heterogeneous. We construct confidence intervals using the estimate of standard deviation as $\sqrt{\widehat{V}_{\textup{LHT}}/n}$ and the appropriate quantiles of the standard normal distribution.

Regression Adjustment for Difference-in-Means Estimator

This section considers the difference-in-means (DM) estimator under complete random assignment, where exactly $n_T$ units are assigned to treatment and $n_C = n - n_T$ to control.

The DM estimator \[ \widehat{\tau}_{\textup{DM}} \;=\; \frac{1}{n_T}\sum_{i=1}^n d_i y_i \;-\; \frac{1}{n_C}\sum_{i=1}^n (1-d_i) y_i \] is unbiased for the average treatment effect under complete random assignment. Its finite-population variance is given by \[ \Var(\widehat{\tau}_{\textup{DM}}) \;=\; \frac{S_n^2(y^{(1)})}{n_T} \;+\; \frac{S_n^2(y^{(0)})}{n_C} \;-\; \frac{S_n^2(\tau)}{n}, \] where $S_n^2(\cdot)$ denotes finite-population variances computed with an $n-1$ divisor (see, e.g., imbens2015causal). Equivalently—and most convenient for our regression-adjustment analysis—this variance admits the representation

where

For a fixed vector ${\bm{b}} \in \mathbb{R}^k$, the covariate-adjusted DM estimator \[ \widehat{\tau}_{\textup{DM}}({\bm{b}}) = \frac{1}{n_T}\sum_{i=1}^n d_i\,(y_i - \eta^{-1}{\bm{x}}_i^\top{\bm{b}}) \;-\; \frac{1}{n_C}\sum_{i=1}^n (1-d_i)\,(y_i - \eta^{-1}{\bm{x}}_i^\top{\bm{b}}), \] where $\eta=\sqrt{n_C/n_T} + \sqrt{n_T/n_C}$, is unbiased for $\tau$, and its variance is equal to

with $\overline{(\widetilde{{\boldsymbol \mu}}-{\bm{X}}{\bm{b}})} = (n^{-1}\sum_{i=1}^n(\widetilde{\mu}_i - {\bm{x}}_i^\top{\bm{b}}))\boldsymbol{1}$ . Therefore, the minimum variance of $\widehat{\tau}_{\textup{DM}}({\bm{b}})$ is equal to

The vector $\widetilde{{\boldsymbol \mu}}$ is not observed, so the variance in (ref) cannot be attained in practice. The next section introduces an estimable surrogate and uses it to define a leave-one-out regression-adjusted DM estimator (LOORA-DM).

Leave-One-Out Regression-Adjusted Difference-in-Means

Similar to LOORA-HT, our proposed leave-one-out regression-adjusted difference-in-means (LOORA-DM) estimator ((ref)) computes a distinct coefficient vector ${\bm{b}}$ for each unit $i$. A key difference, however, is the dependent variable used in the regression to estimate the adjustment coefficients. This regression must be modified to account for the dependence in treatment assignments induced by complete randomization.

For each unit $i$, we define \[ \widetilde{{\bm{y}}}^{(-i)} \;=\; \frac{n_T n_C (n-1)}{n}\,{\bm{f}}^{(-i)} {\bm{y}}, \] where the components of ${\bm{f}}^{(-i)}$ are given by \[ f^{(-i)}_j =

\] if $d_i=1$, and \[ f^{(-i)}_j =

\] if $d_i=0$.

The coefficient vector $\widehat{\bbeta}^{(-i)}_{\lambda}$ is then obtained by regressing $\widetilde{{\bm{y}}}^{(-i)}_{-i}$ on ${\bm{X}}_{-i}$ using ridge regression. Finally, the outcome of unit $i$ is adjusted by ${\bm{x}}_i^\top \widehat{\bbeta}^{(-i)}_{\lambda}$, \[ \widehat{\tau}_{\textup{LDM}}=\frac{1}{n_T}\sum_{i=1}^n d_i\,(y_i - {\bm{x}}_{i}^\top \widehat{\bbeta}^{(-i)}_\lambda) \;-\; \frac{1}{n_C}\sum_{i=1}^n (1-d_i)\,(y_i - {\bm{x}}_{i}^\top \widehat{\bbeta}^{(-i)}_\lambda). \]

The Bias and Variance of LOORA-DM

This section establishes that LOORA-DM is exactly unbiased for the average treatment effect and derives an exact expression for its finite-sample variance. The expression of the variance is more intricate for LOORA-DM than for LOORA-HT because complete randomization induces dependence between treatment assignments. Unbiasedness is also slightly more delicate to verify than in the LOORA-HT case, but it ultimately follows from the fact that $\widetilde{{\bm{y}}}^{(-i)}$ is constructed so that $\E[\widetilde{{\bm{y}}}^{(-i)}_{-i}]=(\sqrt{n_c n_T}/n)\widetilde{{\boldsymbol \mu}}_{-i}$ for each $i$.

The first line in the variance expression is the finite-population variance of $(\widetilde{\mu}_i-{\bm{x}}_i^\top\bbeta_{\lambda})/(1-h_{\lambda ii})$ rescaled by $1/n$. The $1-h_{\lambda ii}$ in the denominator provides a leverage correction that generates from leave-one-out adjustment. The second line of the variance expression in (ref) collects cross-unit terms that include the off-diagonal leverages $h_{\lambda ij}$. These cross-unit contributions did not appear in the LOORA-HT analysis because simple random assignment makes assignments independent across units, whereas complete random assignment induces dependence across units. The matrix $\bQ$ bundles additional corrections that depend only on $(n_T,n_C)$ and the leverage structure encoded in $\bH_{\lambda}$. For a fixed number of covariates, the first term scales as $O(1/n)$, while both the cross-unit term and the $\bQ$-term scale as $O(1/n^2)$.\footnote{To see why the cross-unit term in the result of Theorem (ref) is $O(1/n^2)$, apply Cauchy-Schwarz inequality and note that $(\sum_{\substack{i,j \in [n]: i \neq j}} h_{\lambda ij}^2 )^{0.5} \leq \|\bH\|_{F} \leq \sqrt{k}$.}

Theorem (ref) shows that the variance of LOORA-DM is dominated by the centered dispersion of the leverage-adjusted signal and that dependence induced by complete random assignment contributes only second-order terms. This decomposition is useful to establish asymptotic distributional results and efficiency comparisons.

Asymptotic Normality of LOORA-DM

This section analyzes the large-sample behavior of LOORA-DM under complete random assignment. It shows that, in large samples, LOORA-DM satisfies a central limit theorem. Similar to (ref), we study the behavior of truncated sequences $({\bm{x}}_i, y_i^{(1)}, y_i^{(0)})_{i \in [n]}$ as $n \to \infty$, and assume that the ratio $n_T/n$ is bounded away from zero and one; see (ref).

Recall that for each unit $i$, the adjusted outcome is constructed using a coefficient vector estimated from all other units:

To study the limiting behavior of these regression adjustments, define the following population version of $\widehat{\bbeta}^{(-i)}_\lambda$

$\bbeta_{\lambda}$ represents the deterministic coefficient vector that best captures the relation between covariates and the signal component $\widetilde{{\boldsymbol \mu}}$ in the population.

The variance expression above parallels that of the LOORA-HT estimator but differs in two key aspects. First, the large sample variance of the difference-in-means estimator under complete random assignment automatically recenters ${\bm{X}}$ and $\widetilde{{\boldsymbol \mu}}$. In contrast, LOORA-HT requires the inclusion of an intercept term to obtain a similar variance expression. Second, the scaling of $\widetilde{{\boldsymbol \mu}}$ reflects the finite-sample proportions of treated and control units (equation (ref)).

Theorem (ref) shows that LOORA-DM achieves an asymptotically normal distribution with a variance determined by the projection of the population signal onto the null space of the covariate matrix. It implies that LOORA-DM is as efficient in large samples as the interacted adjustment estimator of lin2013agnostic.

Confidence Intervals for LOORA-DM

The next theorem provides a consistency result for an estimator of the variance of LOORA-DM.

As with LOORA-HT, and more generally with finite-population variance estimation for treatment-effect estimators, $\widehat{V}_{\textup{LDM}}$ is consistent under treatment-effect homogeneity and conservative otherwise. Large-sample confidence intervals for $\tau$ follow by combining this variance estimator with quantiles of the standard normal distribution.

Discussion

LOORA targets three properties that existing estimators do not achieve simultaneously: finite-sample unbiasedness, asymptotic efficiency, and robustness to high-leverage observations. Classical estimators such as Horvitz–Thompson and the difference-in-means are exactly unbiased, but they leave efficiency gains on the table. Regression-adjusted estimators recover those gains asymptotically, yet they introduce finite-sample bias and can become unstable in the presence of high-leverage observations. LOORA delivers all three properties simultaneously: it preserves exact unbiasedness, attains the lin2013agnostic efficiency bound in large samples, and remains stable under influential covariates through its leave-one-out construction.

Table (ref) summarizes the properties for LOORA and competing estimators under simple and complete randomization.

\noindentClassical estimators. Under simple randomization, the Horvitz--Thompson (HT) estimator is exactly unbiased but can exhibit high variance. Regression adjustment aims to improve efficiency by controlling for observed covariates. The most commonly used adjusted estimator, $\widehat\tau_{\text{ADJ}}$, regresses ${\bm{y}}$ on ${\bm{X}}$ and the treatment indicator ${\bm{d}}$ with intercept, producing a consistent estimator only when treatment effects are homogeneous. The interactive version, $\widehat\tau_{\text{INT}}$, adds treatment-by-covariate interactions and is asymptotically efficient under both simple (with equal treatment probability assignment for all units) and complete randomization lin2013agnostic. However, both estimators, $\widehat\tau_{\text{ADJ}}$ and $\widehat\tau_{\text{INT}}$, can become unstable when the sample contains high-leverage observations, and neither is unbiased.

\noindentLeave-one-out regression adjustment. The LOORA framework modifies regression adjustment to preserve finite-sample unbiasedness without sacrificing efficiency. By using leave-one-out fitted values, the LOORA-HT and LOORA-DM estimators remove the bias induced by regression adjustment while maintaining asymptotic optimality. Under simple randomization, LOORA-HT achieves the same variance bound as the lin2013agnostic estimator but is exactly unbiased and robust to leverage. Under complete randomization, LOORA-DM extends these results to dependent assignment structures by appropriately rescaling the leave-one-out predictions. Both estimators combine the finite-sample exactness of HT and DM with the asymptotic efficiency of regression adjustment.

\noindentRelation to recent work. Recent work has proposed innovative methods for covariate adjustment in randomized experiments. ghadiri2023finite develops a cross-fitting unbiased regression-adjusted estimator for simple random assignment. This approach, however, does not attain the lin2013agnostic variance in large samples. harshaw2022design, in contrast, takes a design-based perspective: rather than modifying the estimator, they optimize the randomization mechanism itself to minimize the variance of the unadjusted Horvitz--Thompson estimator. Finally, spiess2025optimal characterizes the necessary conditions under which design-based estimators such as Horvitz--Thompson or difference-in-means are unbiased. For the HT estimator, these conditions require a leave-one-out structure, which our LOORA-HT estimator satisfies by construction. For the DM estimator, unbiasedness requires a leave-two-out structure. We show below that LOORA-DM can be represented as a leave-two-out estimator, thereby satisfying Spiess’s condition. Hence, our estimators are not only consistent with the theoretical framework of spiess2025optimal, but also provide constructive, closed-form realizations of unbiased estimators that satisfy these necessary design-based criteria. The remainder of this section provides a detailed comparison with prior work.

Comparison with Cross-Fitted Regression-Adjusted HT

ghadiri2023finite proposes a finite-sample unbiased regression-adjusted Horvitz–Thomp-\allowbreak son estimator built on cross-fitting. The procedure splits the sample into two groups, learns a regression vector within each group, and then uses each group’s fitted vector to adjust the outcomes in the other. ghadiri2023finite considers only the case of constant $p_i=0.5$ for all $i$ and derives a high-probability upper bound on the variance of their estimator. Specifically, conditional on an event $\mathcal{E}$ that holds with probability at least $1-\delta$, the variance of the estimator in ghadiri2023finite is bounded by

where $\zeta_{{\bm{X}}}^2 = \left\|{\bm{X}}\right\|_{2, \infty}^2$ and $\varepsilon$ is an approximation error that depends on $\delta$.\footnote{$\mathcal{E}$ is the event under which the leverage-score sampling in ghadiri2023finite delivers a $(1+\epsilon)$-approximate solution to the full regression problem. Concretely, the method randomly subsamples rows of ${\bm{X}}$ and solves the regression on the subsample rather than on all $n$ observations. The resulting approximation error relative to the full-sample regression induces the $(1+\epsilon)$ factor in (ref)}. By (ref), ghadiri2023finite choice of regularization parameter $\lambda = 100 \log(n/\delta)\, \zeta_{{\bm{X}}}^2$, the ridge leverage scores satisfy $h_{\lambda i i}({\bm{X}}) \leq 1/(1 + 100 \log(n/\delta))$. Moreover, since all $p_i = 0.5$, it follows that $ \max_{i \in [n]} \left\|\widetilde{{\bm{x}}}_i\right\|_{2} = 2 \max_{i \in [n]} \left\|{\bm{x}}_i\right\|_{2}$. Consider

Note that for any ${\bm{b}}\in\R^k$, \[ \left\| {\bm{X}} {\bm{b}} - {\boldsymbol \mu} \right\|_{2}^2 + 100 \log(n/\delta)\, \zeta_{{\bm{X}}}^2 \left\| {\bm{b}} \right\|_{2}^2 = \left\| \widetilde{{\bm{X}}} ({\bm{b}}/2) - {\boldsymbol \mu} \right\|_{2}^2 + 100 \log(n/\delta)\, \zeta_{\widetilde{{\bm{X}}}}^2 \left\| {\bm{b}} / 2 \right\|_{2}^2. \] Therefore, \[ \mathop{\mathrm{argmin}\,}\limits_{{\bm{b}} \in \mathbb{R}^k} \Big( \left\| \widetilde{{\bm{X}}} {\bm{b}} - {\boldsymbol \mu} \right\|_{2}^2 + 100 \log(n/\delta)\, \zeta_{\widetilde{{\bm{X}}}}^2 \left\| {\bm{b}} \right\|_{2}^2 \Big) = \bbeta/2. \]

For $\lambda = 100 \log(n/\delta)\, \zeta_{\widetilde{{\bm{X}}}}^2$ with $\delta\leq \exp(-1)$, the variance of LOORA-HT from (ref), is bounded by

Relative to the bound of ghadiri2023finite in (ref), the bound above improves the first term by roughly a factor of eight. It also eliminates the additional ridge-regularization penalty term that appears in (ref) and improves the final term by about a factor of sixteen.

In addition, unlike the procedure in ghadiri2023finite, LOORA comes with consistent variance estimation and valid confidence intervals.

Comparison with Gram--Schmidt Random Walk Design

harshaw2022design proposes an ingenious design-based approach to variance reduction. Rather than modifying the estimator, they optimize the randomization mechanism itself through a Gram-Schmidt Walk (GSW) design that ensures covariate balance while controlling the marginal assignment probabilities. One potential limitation of GSW is that it requires covariate information on all sample units prior to treatment assignment. In contrast, simple or complete random assignment can be implemented online (sequentially) as additional experimental subjects are recruited for a study. Moreover, with simple or complete randomization, a researcher can freely add or remove covariates, analyze different subgroups, or draw additional samples even after the experiment has been conducted. This flexibility, together with the transparency and ease of implementation of standard randomization, makes simple and complete randomization attractive in many applied settings. Motivated by these considerations, LOORA-HT and LOORA-DM are explicitly designed for simple and complete randomization, respectively.

harshaw2022design shows that the GSW design achieves the following variance bound for any choice of $\phi \in (0,1)$:

We now compare this bound with the LOORA-HT variance. The second term in (ref) and the corresponding regularization term in (ref) are not directly comparable in general, although both scale on the order of $k/n^2$. For the first term, with a constant probability of assignment, \[ \min_{\bbeta \in \mathbb{R}^k} \Big[ \frac{1}{\phi}\left\|{\boldsymbol \mu} - {\bm{X}} \bbeta\right\|_{2}^2 + \frac{\zeta_{{\bm{X}}}^2}{1 - \phi}\left\|\bbeta\right\|_{2}^2\Big] = \min_{\bbeta \in \mathbb{R}^k} \Big[ \frac{1}{\phi}\left\|{\boldsymbol \mu} - \widetilde{{\bm{X}}} \bbeta\right\|_{2}^2 + \frac{\zeta_{\widetilde{{\bm{X}}}}^2}{1 - \phi}\left\|\bbeta\right\|_{2}^2\Big]. \] Therefore, setting $\lambda = \zeta_{\widetilde{{\bm{X}}}}^2 / (1 - \phi)$, (ref) implies that \[ \widetilde h_{\lambda ii} \;\leq\; \frac{1 - \phi}{2 - \phi}. \] Consequently, \[ \frac{1}{(1 - \widetilde h_{\lambda ii})^2} \;\leq\; (2 - \phi)^2. \] For $\phi \in [0, (3 - \sqrt{5})/2)$, we have $(2 - \phi)^2 \leq 1/\phi$. Hence, the scaling factor that appears in our variance bound is uniformly smaller than that of the GSW design. Moreover, the first term in (ref) diverges as $\phi \to 0$, whereas the corresponding constant in the variance of LOORA-HT remains bounded.

Finally, when the leverage scores of $\widetilde{{\bm{X}}}$ are small, LOORA can set $\lambda=0$. In contrast, under the GSW design the smallest feasible regularization parameter (i.e., coefficient of $\left\|\bbeta\right\|_{2}^2$ in (ref)) is $\lambda=\zeta_{{\bm{X}}}^2$.

Leave-One-Out and Leave-Two-Out Estimators

aronow2013class studies a general class of unbiased Horvitz-Thompson estimators. These estimators adjust unit $i$'s outcome by a function $f({\bm{x}}_i,\btheta_i)$, where $\btheta_i$ is a vector of parameters. The estimator is unbiased when $d_i$ is uncorrelated with $f({\bm{x}}_i,\btheta_i)$. A useful case arises when $\btheta_i$ depends only on outcomes $y_j$ whose treatment indicators $d_j$ are independent of $d_i$. Under independent treatment assignment, this condition suggests leave-one-out constructions. LOORA-HT follows this logic by setting $\btheta_i=\widehat{\bbeta}^{(-i)}_{\lambda}$, as defined in (ref).

wu2018loop studies a special case of aronow2013class in which the adjustment term for observation $i$ is constructed as a linear combination of two counterfactual outcomes for unit $i$ (obtained by any method of our choosing, such as linear regression or random forests)---a predicted outcome for unit $i$ under treatment and a predicted outcome under control---via a leave-one-out fit, which satisfies the sufficient condition for $f$ stated in aronow2013class.

Building on these contributions, spiess2025optimal provides a general characterization of unbiased regression-adjusted estimators for average treatment effects under arbitrary randomization schemes when the treatment-assignment probability is the same for all units. He shows that for Horvitz--Thompson and difference-in-means type estimators, unbiasedness imposes specific structural constraints on how each unit’s outcome can depend on treatment assignments of other units. In particular, under simple random assignment with equal treatment assignment probability for all units, the Horwitz--Thompson estimator must have a leave-one-out structure—each adjusted outcome for unit $i$ may depend on the treatment assignments of all other units except $i$—whereas under complete random assignment, the estimator must satisfy a stronger leave-two-out condition, meaning that each pairwise contribution between units $i$ and $j$ can depend only on the treatment assignments of units other than $i$ and $j$.

We next show that LOORA-DM satisfies Spiess's leave-two-out requirement and therefore meets the necessary conditions for unbiasedness in spiess2025optimal. Specifically, we show that the LOORA--DM estimator admits the representation \[ \widehat{\tau}_{\textup{LDM}} = \frac{1}{n_T n_C} \sum_{i<j} (d_i - d_j) \big(y_i - y_j - \phi_{ij}({\bm{z}}_{-ij})\big), \] where $\phi_{ij}$ is an adjustment function and ${\bm{z}}_{-ij}$ collects the treatment assignments of all units other than $i$ and $j$. Note that \[ \frac{1}{n_T n_C} \sum_{i<j} (d_i - d_j)(y_i - y_j) = \frac{1}{n_T}\sum_{d_i=1} y_i - \frac{1}{n_C}\sum_{d_i=0} y_i, \] so it suffices to characterize the adjustment term $\phi_{ij}({\bm{z}}_{-ij})$.

\paragraph{Pairwise representation.} We define \[ \widetilde{y}_{\ell} =

\] and set \[ \phi_{ij}({\bm{z}}_{-ij}) = {\bm{x}}_i^\top ({\bm{X}}_{-i}^\top {\bm{X}}_{-i} + \lambda {\bm{I}})^{-1} {\bm{X}}_{-ij}^\top \widetilde{{\bm{y}}}_{-ij} - {\bm{x}}_j^\top ({\bm{X}}_{-j}^\top {\bm{X}}_{-j} + \lambda {\bm{I}})^{-1} {\bm{X}}_{-ij}^\top \widetilde{{\bm{y}}}_{-ij}. \] This form makes explicit that $\phi_{ij}$ depends only on ${\bm{z}}_{-ij}$—that is, on the treatment assignments of all units other than $i$ and $j$—and hence defines a leave-two-out structure.

Substituting this definition, we obtain

Some algebra shows that these expressions are exactly the adjustment values in (ref) for both treated and control units. Thus, $\widehat{\tau}_{\textup{LDM}}$ admits a leave-two-out representation. (ref) gives a detailed proof.

Evidence from Within-Subject Experiments

This section presents empirical evidence on the finite-sample performance of LOORA estimators. The analysis uses experimental data from two within-subject studies. In a within-subject experiment, each unit experiences both the treatment and control conditions, so both potential outcomes are observed for every unit. We use within-subject experiment data to simulate randomized evaluations with known ground truth and with a realistic joint distribution of potential outcomes.

The datasets are drawn from allcott2015evaluating and mcdonald2025evaluating. A brief description of each is provided below.

We simulate randomized assignments by masking one outcome per unit according to a prescribed randomization rule. This setup lets us assess each estimator’s finite-sample performance, bias, standard deviation, RMSE, and the average coverage of confidence intervals, under realistic distributions for potential outcomes and covariates. We consider three treatment-assignment mechanisms. In the first two designs, treatment assignments are independent across units. In the first design, treatment assignment probabilities depend on covariate values. Specifically, we draw a random Gaussian vector in $\R^k$, compute for each unit $i$ the cosine similarity $c_i$ between this vector and its covariate vector, and set \[ p_i = \max\!\left\{\min\!\left\{\frac{1+c_i}{2},\, 0.8\right\},\, 0.2\right\}. \] The second design is simple random assignment, with all treatment assignment probabilities equal to $0.5$. The third design is complete random assignment, with half the units assigned to treatment.

Simulation results are based on $100{,}000$ repetitions. (ref) report results for covariate-dependent assignment; (ref) for simple random assignment with equal probabilities; and (ref) for complete random assignment. For simple random assignment, we compare our estimator, LOORA--HT, with two benchmarks: (i) linear regression-adjustment, and (ii) the regression adjustment of lin2013agnostic with additional interacted terms. For complete random assignment, in addition to these, we also compare our estimator, LOORA--DM, with the exact bias correction approach of chang2024exact, which we refer to as “Debiased” in (ref).

We report the coverage probabilities of confidence intervals in (ref) using heteroskedasticity-consistent (HC) variance estimators. HC0 denotes White's original HC estimator white1980heteroskedasticity, while HC2 denotes the leverage-adjusted HC estimator of mackinnon1985some. The variance estimator in (ref) corresponds to HC0, whereas the variance estimator in (ref) corresponds to HC2.

All simulations use recentered covariates and an intercept. That is, regression adjustment uses the matrix of regressors \[ \bM =

. \] The value of $\zeta$ in the tables is equal to $\left\|\bM\right\|_{2,\infty}$.

Throughout Tables (ref)--(ref), LOORA estimators consistently achieve lower bias, higher precision, and lower RMSE than competing methods. Moreover, the confidence intervals for LOORA estimators based on heteroskedasticity-robust standard errors exhibit coverage rates remarkably close to the nominal level, whereas the alternative methods tend to undercover substantially. These results demonstrate the potential role of LOORA and ridge regularization to stabilize regression adjustment without sacrificing unbiasedness or inferential validity.

Conclusion

This article develops leave-one-out regression-adjusted estimators of average treatment effects under simple and complete random assignment. We show that these estimators are exactly unbiased in finite populations, achieve the asymptotic efficiency of the interacted regression of lin2013agnostic, and admit closed-form expressions for their variance. Regularization through leave-one-out ridge adjustment ensures robustness to high-leverage observations and stabilizes inference even in small samples. Our theoretical results establish a design-based foundation for regression adjustment that unifies classical unbiasedness, modern efficiency results, and practical inference procedures. Empirical evaluations confirm the estimators’ excellent finite-sample performance across different designs.

Natural next steps include extending the framework to stratified or covariate-adaptive designs and to clustered or networked experiments, and incorporating high-dimensional or nonparametric adjustments. Taken together, these directions define a broad agenda for design-based regression adjustment in modern experimental settings.