Regression Adjustment, Cross-Fitting, and Randomized Experiments with Many Controls

Regression adjustment in completely randomized experiments with many covariates

\address{Department of Economics, University of Wisconsin-Madison, 1180 Observatory Drive Madison, WI 53706-1393, USA.} \email{[email removed]} \address{Graduate School of Economics, Hitotsubashi University, 2-1 Naka, Kunitachi, Tokyo 186-8601, Japan.} \email{[email removed]} \address{Department of Economics, London School of Economics, Houghton Street, London, WC2A 2AE, UK.} \email{[email removed]}

Introduction

Randomized controlled trials (RCTs) remain among the most fundamental and influential tools for causal inference, widely employed by empirical researchers across the natural, social, and biomedical sciences. See, for instance, \citet*{fisher1925statistical,fisher1935}, \citet*{neyman1923application}, and \citet*{kempthorne1952design} for foundational developments, and \citet*{imbens2015causal} and \citet*{rosenberger2015randomization} for modern textbook treatments. Statistical inference in RCTs is typically approached from one of two distinct perspectives: the finite population and the superpopulation frameworks. First introduced by \citet*{neyman1923application}, the finite population perspective treats potential outcomes as fixed, with randomness arising solely from the treatment assignment mechanism. In contrast, the superpopulation approach assumes that observed units are independently drawn from a hypothetical infinite population. While both perspectives are influential and widely adopted, theoretical developments under the finite population framework---especially in more complex settings---remain relatively underexplored. This paper adopts the finite population perspective and examines causal inference using regression adjustment methods for RCTs.\footnote{It is not our intention to advocate for either perspective; see \citet*{reichardt1999justifying} for a philosophical comparison of the two.}

In various RCTs, researchers usually collect covariates that are predetermined characteristics of the experimental subjects and conduct regression adjustments to estimate treatment effects of interest since regression adjustments can potentially reduce variability of the estimates (see, for example, Section 7 in Imbens and Rubin, 2015). However, different opinions exist on whether to adjust for covariates; in an influential work, \citet*{freedman2008regressiona, freedman2008regression} criticized the practice of using regression adjustment for RCTs with three critiques: (i) lack of efficiency guarantee of ad hoc regression adjustment over the unadjusted estimator, (ii) inconsistency of the classical regression variance estimator, and (iii) presence of a bias term of order $O_{p}(n^{-1})$. When the number of covariates is treated as fixed, the first two critiques have been addressed by \citet*{lin2013agnostic}, in which the author suggested running a regression of the observed outcomes on the treatment variable, covariates, and their interactions. { \citet*{lin2013agnostic} showed that (i) this regression adjusted estimator is consistent and asymptotically more efficient than the simple difference in means estimator without regression adjustment, and (ii) this efficiency improvement is ensured despite arbitrary misspecification of the conditional mean function (called the no-harm property; see also \citet*{negi2021} for analogous results in the superpopulation framework).} In addition, \citet*{lin2013agnostic} showed that the heteroskedasticity robust variance estimators for linear regression is asymptotically conservative and thus provides valid size control. Recently, \citet*{chang2021exact} address the remaining criticism by providing analytic exact bias correction formulae for the regression adjustment estimators in \citet*{freedman2008regressiona, freedman2008regression} and \citet*{lin2013agnostic}. Thus far, at least under the asymptotic framework where the number of covariates held fixed, Freedman's critiques on regression adjustment for RCTs have been addressed.

In addition to these remarkable progresses, attempts have been made to study asymptotic regimes that allow the number of covariates to grow with the population size. Such analyses are empirically important because in many RCT studies, researchers record a sizable set of covariates whose dimension is often not negligible compared to the number of experimental subjects. Indeed, in such scenarios, theoretical guarantees derived under fixed dimensionality may be far less than compelling; with a sizable number of covariates, the bias, oftentimes non-negligible, becomes even more problematic. In such asymptotic environments, an important recent contribution came from \citet*{lei2021regression}; under fairly mild conditions, they establish asymptotic normality permitting growing number of covariates, and characterize the leading term of the bias for the regression adjustment estimator of \citet*{lin2013agnostic}. They go one step further by providing an analytic bias-correction estimator. Despite its promising theoretical guarantees, their proposed bias-corrected estimator does not appear to be nearly bias free in their simulation studies when the DGPs contain more nonlinearity as well as larger numbers of covariates. As a practical solution, they further recommend a trimming procedure for covariates to get around the unreliable finite sample bias performances of their bias-corrected estimator. Nevertheless, the means to effectively tackle the bias problem without resorting to artificial modification of the covariates remain unclear.

In this paper, we contribute to the endeavor of understanding regression adjustment in multiple fronts. { First, we study higher-order properties of Lin's (2013) regression adjustment estimator, Lei and Ding's (2021) bias-corrected estimator, and the cross-fitted estimator proposed by aronow-middleton2013 when the number of covariates grows at a slower rate than the sample size under completely randomized experiments, and show that the cross-fitted estimator possesses improved asymptotic properties compared to the others.} Second, we derive a finer asymptotic variance expression for the estimators that takes into account of the higher-order term. As pointed out in Lei and Ding (2021, Section 4.3), the asymptotic variance of the regression adjustment estimators can deviate significantly from the theoretical ones in finite samples, especially when the dimensionality and/or nonlinearity in the DGPs is non-negligible. This further motivates us to propose an alternative bias-corrected version of the HC3 standard error. The simulation studies unveil supporting evidences that the cross-fitted estimator has favorable performances robustly over a variety of scenarios. Coupled with our bias-corrected HC3, it delivers more precise inference results than existing alternative estimation and inference methods when researchers utilize a modest or large number of covariates for causal inference in RCTs. Our methodology is also extended to cover stratified experiments with large strata.

In both social and natural sciences, researchers often find RCTs involve a sizable number of available covariates in their empirical applications. To formally cope with such settings, \citet*{bloniarz2016lasso} and \citet*{wager2016high} studied regression adjustments by machine learning techniques in a high-dimensional setup where the dimensionality $p$ may be larger than the population size $n$. On the other hand, \citet*{lei2021regression} investigated the situation where $p\ll n$ but $p$ may grow with $n$, and developed a bias correction method for the regression adjustment estimator; as eloquently reasoned by \citet*{lei2021regression}, this moderately growing $p$ asymptotics is of particular importance in a wide range of applications that involve RCTs. { This paper employs the same setup as \citet*{lei2021regression} and focuses on the case of $p\ll n$.}

Relationship to the literature

This paper is built upon a growing body of the important recent forays into innovating theory of RCTs under finite population asymptotics; these include but are not limited to, \citet*{freedman2008regressiona, freedman2008regression}, \citet*{lin2013agnostic}, \citet*{tan2014second}, \citet*{aronow2014sharp}, \citet*{dasgupta2015causal}, \citet*{bloniarz2016lasso}, \citet*{wager2016high}, \citet*{fogarty2018regression}, \citet*{li2018asymptotic}, \citet*{abadie2020sampling}, \citet*{li2020rerandomization}, \citet*{chang2021exact}, \citet*{imbens2021causal}, and \citet*{lei2021regression}. In particular, \citet*{tan2014second} studies the higher-order asymptotics for regression in various design-based setups assuming fixed covariate dimensionality. { In a recent independent work, \citet*{lu2025} proposed an alternative bias-corrected estimator and studied asymptotic properties when $p\ll n$ and $p\sim n$. For the case of $p\ll n$ considered in this paper, we investigate asymptotic properties of different estimators up to second-order terms.} It is also closely related to the studies of regression models with many regressors under superpopulation setups such as, e.g. \citet*{cattaneo2018alternative,cattaneo2018inference,cattaneo2019two}, to list a few. The idea of cross-fitting or sample splitting has been widely applied in causal inference literature; in fact, it is a common strategy to reduce bias terms in many semiparametric and high-dimensional models, see, e.g., \citet*{schick1986asymptotically}, \citet*{zheng2011cross}, \citet*{chernozhukov2018double}, \citet*{newey2018cross}, \citet*{spiess2018optimal}, \citet*{bradic2019sparsity}, to list a few.

{ The idea of cross-fitting is well established in RCT contexts. For instance, aronow-middleton2013 and wu2018loop discuss unbiased estimation of the average treatment effect using robust moment conditions and leave-one-out procedures (see also williams1961 and robins1999 for related antecedents). In contrast to aronow-middleton2013, who derive unbiased estimators under independent treatment assignments, this paper treats the number of treated units as deterministic and accommodates potential dependence among treatment assignments. By incorporating the treatment-covariate interactions highlighted by lin2013agnostic, the cross-fitted estimator considered in this paper not only retains the low-bias advantages of cross-fitting but also performs robustly in settings with a large number of covariates.}

Our work sheds new light on these literatures by providing a bias-corrected estimation procedure that combines the idea of cross-fitting and efficient regression-assisted estimation for RCTs, and further establishes formal theoretical justification for its advantages in performances for models in RCTs with large numbers of covariates under design-based finite population asymptotics.\footnote{ In a recent working paper by matsushita2025, two of the authors of this paper studied analogous issues under the sampling-based superpopulation setup, and developed an optimal covariates selection criterion and higher-order accurate standard error.}

Methodology

Consider a treatment-control RCT, where $y_{i}(1)$ and $y_{i}(0)$ are potential outcomes of unit $i=1,\ldots,n$ for treatment and control, respectively, and $T_{i}$ is an indicator for assignment ($T_{i}=1$ corresponds to the treatment, and $T_{i}=0$ corresponds to the control). { This paper focuses on the completely randomized experiment, where a researcher randomly assigns $n_1$ units to the treatment group and $n_0=n-n_1$ units to the control group (see, Chapter 4.4 of Imbens and Rubin, 2015). More precisely, the experimenter deterministically chooses $n_1$ and $n_0$, and treatment units are randomly drawn size-$n_1$ subset of $\{1,\ldots,n\}$ uniformly over all possible $n!/(n_1!n_0!)$ subsets. We note that this setup is commonly applied in the RCT literature using the design-based analysis (e.g., freedman2008regressiona, lin2013agnostic, and lei2021regression discussed above), and that the treatment variables $(T_1,\ldots,T_n)$ are dependent in contrast to the iid sampling.}\footnote{ aronow-middleton2013 adopted the iid sampling to study unbiasedness of their estimation method that covers the cross-fitted estimator defined below. The focus of this paper is the higher-order asymptotic property of the cross-fitted estimator with the growing number of covariates under the completely randomized experiment.}

We wish to conduct estimation and inference on the average treatment effect \[ \tau=\mu_{1}-\mu_{0},\quad \mbox{where }\mu_{t} = \frac{1}{n}\sum_{i=1}^{n}y_{i}(t)\mbox{ for }t=0,1, \] based on the observed outcome \[ Y_{i}=y_{i}(T_{i}), \] and $(p-1)$-dimensional pretreatment covariates $x_{i}$. In this paper, we employ the finite population perspective neyman1923application, where the potential outcomes $y_{i}(1)$ and $y_{i}(0)$ are non-random and randomness comes solely from the treatment indicator $T_{i}$ (see, e.g., \citet*{imbens2015causal}, for an overview).

The simplest estimator of $\tau$ is the difference in means \[ \hat{\tau}^{\mathrm{dif}}=\frac{1}{n_{1}}\sum_{i=1}^{n}T_{i}Y_{i} - \frac{1}{n_{0}}\sum_{i=1}^{n}(1-T_{i})Y_{i}, \] where $n_{1}$ and $n_{0}$ are the sizes of the treatment and control groups, respectively. Although this estimator is unbiased and asymptotically normal, \citet*{lin2013agnostic} showed that a regression adjustment using $x_{i}$ yields a more efficient estimator than $\hat{\tau}^{\mathrm{dif}}$. This regression adjustment estimator $\hat{\tau}^{\mathrm{adj}}$ is obtained as the OLS coefficient on $T_{i}$ from the regression of $Y_{i}$ on $(1,T_{i},(x_{i}-\bar{x})^{\prime},T_{i}(x_{i}-\bar{x})^{\prime})$, where $\bar{x}=n^{-1}\sum_{i=1}^{n}x_i$. To facilitate our discussion on bias correction below, we present an alternative expression for $\hat{\tau}^{\mathrm{adj}}$. Let $z_{i}=(1,x_{i}^{\prime})^{\prime}$, where $\bar{x}$ is normalized to be zero for each coordinate, and $\hat{\beta}_{1}$ and $\hat{\beta}_{0}$ be the OLS estimators for the regression of $Y_{i}$ on $z_{i}$ by the treatment ($T_{i}=1$) and control ($T_{i}=0$) groups, respectively. Then the regression adjustment estimator can be written as \[ \hat{\tau}^{\mathrm{adj}}=\hat{\mu}_{1}^{\mathrm{adj}}-\hat{\mu}_{0}^{\mathrm{adj}},\quad \mbox{where }\hat{\mu}_{t}^{\mathrm{adj}}=\frac{1}{n}\sum_{i=1}^{n}z_{i}^{\prime}\hat{\beta}_{t}\mbox{ for }t=0,1. \] \citet*{lin2013agnostic} showed that $\hat{\tau}^{\mathrm{adj}}$ is consistent, asymptotically normal, and more efficient than the difference in means $\hat{\tau}^{\mathrm{dif}}$. It should be noted that these results hold true under the finite population setup with fixed $p$ without assuming correct specification of the linear model.

In practice, it is often the case that researchers observe many covariates. \citet*{lei2021regression} studied asymptotic properties of the regression adjustment estimator when the number of covariates $p$ grows with the sample size, and developed a bias-corrected estimator. To define Lei and Ding's (2021) approach, we introduce some notation. Let $Z=(z_{1},\ldots,z_{n})^{\prime}$, $P_{ij}$ be the $(i,j)$-th element of the projection matrix $P=Z(Z^{\prime}Z)^{-1}Z^{\prime}$, and $\hat{e}_{i}$ be the OLS residual, that is \[ \hat{e}_{i} =

. \] Lei and Ding's (2021) bias-corrected estimator for $\tau$ is defined as

where \[ \hat{\Delta}_{1}=n_{1}^{-1}\sum_{i=1}^{n}T_{i}P_{ii}\hat{e}_{i}, \quad \hat{\Delta}_{0}=n_{0}^{-1}\sum_{i=1}^{n}(1-T_{i})P_{ii}\hat{e}_{i}. \] Note that $\frac{n_{0}}{n_{1}}\hat{\Delta}_{1}$ and $\frac{n_{1}}{n_{0}}\hat{\Delta}_{0}$ are correction terms to estimate the higher-order bias terms of $\hat{\mu}_{1}^{\mathrm{adj}}$ and $\hat{\mu}_{0}^{\mathrm{adj}}$ under the moderate-$p$ asymptotics, respectively. The terms involving $\hat{\Delta}_{t}$ are analytic bias estimates that replace the unknown bias terms in their asymptotic theory. Although this bias correction method works in theory, the quality of these bias estimates may not be ideal, as illustrated in Section 4.4 of \citet*{lei2021regression}.

{ This paper studies an alternative bias correction approach via cross-fitting adapted from aronow-middleton2013. We first note that the regression adjustment estimators for $\mu_{1}$ and $\mu_{0}$ can be alternatively written as

where $\pi=n_{1}/n$ is the fraction of treated units treated as non-random in our setup.} Albeit the implementation differences, the estimation based on ((ref)) is equivalent to the full-sample regression adjustment estimation with treatment-covariate interactions first proposed by \citet*{lin2013agnostic} and the regression adjustment estimator in \citet*{lei2021regression}. The key idea of the bias correction is to replace the OLS estimators $\hat{\beta}_{1}$ and $\hat{\beta}_{0}$ with their leave-one-out counterparts

Then the cross-fitted estimator of the average treatment effect $\tau$ is defined as \[ \hat{\tau}^{\mathrm{cf}}=\hat{\mu}_{1}^{\mathrm{cf}}-\hat{\mu}_{0}^{\mathrm{cf}}, \] where

Although this estimator may appear to be computationally demanding, in practice, one may utilize the identity for leave-one-out OLS estimation (see, e.g., Theorem 3.7 in \citet*{hansen2022econometrics}):

for $i\in\{1,...,n:T_{i}=t\}$, where $\tilde{e}_{i}=\hat{e}_{i}/(1-P_{t,ii})$ for $i$ with $T_i=t$, $P_{t,ij}$ is the $(i,j)$-th entry of the matrix $P_{t}=Z_{t}(Z_{t}'Z_{t})^{-1}Z_{t}'$, and $Z_{t}$ is the $n_{t}\times p$ submatrix that consists of $n_{t}$-rows of matrix $Z$ with $T_{j}=t$. This identity significantly lessens the computational burden to implement the cross-fitted estimator.

Asymptotic theory

In this section, we study asymptotic properties of the cross-fitted estimator $\hat{\tau}^{\mathrm{cf}}$ to compare with the existing ones, $\hat{\tau}^{\mathrm{adj}}$ and $\hat{\tau}^{\mathrm{bc}}$, and associated variance estimators.

Stochastic expansion

We first establish stochastic expansions for the estimators of $\tau$. To this end, we consider the setup employed by Lei and Ding (2021), where the number of covariates $p$ is allowed to grow with the sample size $n$. Let \[ e_{i}(t)=y_{i}(t)-z_{i}^{\prime}\beta_{t},\quad \mbox{where }\beta_{t}=(Z^{\prime}Z)^{-1}Z^{\prime}Y(t) \mbox{ for }t=0,1, \] and \[ \kappa=\max_{1\le i\le n}P_{ii},\quad \mathcal{E}_{2}=\max_{t\in\{0,1\}}\frac{1}{n}\sum_{i=1}^{n}e_{i}(t)^{2},\quad \mathcal{E}_{\infty}=\max_{t\in\{0,1\}}\max_{1\le i\le n}|e_{i}(t)|. \] We impose the following assumptions.

Assumptions (i)-(iv) are identical to Assumptions 1-4 in \citet*{lei2021regression}, respectively. { It should be noted that we impose no assumption on the functional forms of the outcome regression functions.} Assumption (i) holds if the proportions of treatment and control groups are fixed.

{ Assumption (ii) restricts the growth rate of $p$. Since $\kappa \in [p/n,1]$, this assumption implies $p = o(n)$, i.e., $p$ should grow slower than $n$. Also this assumptions allows influential observations as long as their leverages are of smaller orders than $1/\log p$.} { Note that $0 \le P_{ii} \le 1$ and $\sum^{n}_{i=1}P_{ii}=p$. Thus, in the favorable case where all leverage values are close to their average $p/n$, this condition holds if $\frac{p \log p}{n}=o(1)$.} Assumption (iii) imposes a mild restriction on the correlation between the potential residuals from the population ordinary least squares. It rules out perfectly negative correlation between the treatment and control potential residuals. Finally, Assumption (iv) imposes a Lindeberg-Feller type condition that none of potential residual dominates the others, while permitting heavy-tailed outcomes with $\mathcal{E}_{2}$ growing with $n$.

Let

Under the above assumptions, higher-order asymptotic properties of the estimators for the ATE $\tau$ are obtained as follows. {

}

Theorem (ref) (i) decomposes the estimation errors into a first-order dominant linear term $\mathcal{L}$, a second-order quadratic term $\mathcal{W}$, and the bias terms $B^{\mathrm{adj}}$ for $\hat{\tau}^{\mathrm{adj}}$. Note that the linear and quadratic terms are identical for all the estimators, and the differences are attributed to the bias term and stochastic orders of the remainder terms.

The bias term $B^{\mathrm{adj}}$ for the conventional regression adjustment estimator is studied by \citet*{lei2021regression}. In contrast, Lei and Ding's (2021) bias-corrected estimator $\hat{\tau}^{\mathrm{bc}}$ and the cross-fitted estimator $\hat{\tau}^{\mathrm{cf}}$ do not involve such a bias term and have better higher-order bias properties. Since the stochastic term $\mathcal{L}+\mathcal{W}$ is identical for all estimators, such bias reducing features of $\hat{\tau}^{\mathrm{bc}}$ and $\hat{\tau}^{\mathrm{cf}}$ do not inflate the variance compared to the conventional regression adjustment estimator. { However, our theorem does not necessarily mean $\hat{\tau}^{\mathrm{cf}}$ is exactly unbiased: the remainder term in ((ref)) is typically small but non-zero.} Furthermore, compared to \citet*{lei2021regression}, our expansions also characterize the second order quadratic term $\mathcal{W}$, which will be useful to investigate higher-order properties of the variance estimators in the next subsection.\footnote{ One may also consider a leave-one-out version of the regression adjustment estimator $\tilde{\tau}^{\mathrm{adj}} = \tilde{\mu}_{1} - \tilde{\mu}_{0}$, where $\tilde{\mu}_{t} = n^{-1}\sum_{i=1}^{n}z^{\prime}_{i}\hat{\beta}^{(i)}_{t}$. However, this estimator possesses an analogous bias term to the regression adjustment estimator. To see this, the identity in ((ref)) and an analogous argument in the proof of Theorem (ref) yield \[ \tilde{\mu}_{1} = \hat{\mu}^{\mathrm{cf}}_{1} - \frac{1}{n}\sum_{i=1}^n \frac{T_i}{\pi}z_i^{\prime}(Z_1^{\prime}Z_1)^{-1}z_i\tilde{e}_i = \left\{\hat{\mu}^{\mathrm{cf}}_{1} - \frac{1}{n\pi}\sum_{i=1}^n \frac{T_i}{\pi}P_{ii}e_i(1)\right\}(1+o_p(1)). \] Thus, the bias term of $\tilde{\tau}^{\mathrm{adj}}$ will be of same order as $B^{\mathrm{adj}}$.}

We now compare $\hat{\tau}^{\mathrm{bc}}$ and $\hat{\tau}^{\mathrm{cf}}$. As shown in Theorem (ref) (i), the remainder terms of these exhibit different orders. Also Theorem (ref) (ii) characterizes the stochastic components $\mathcal{L}$ and $\mathcal{W}$. Combining these results, the stochastic orders of the estimation errors are

{ We note that $\mathcal{L}$ is asymptotically normal with mean $0$ and variance $\sigma_L^2/n$; see Lei and Ding (2021, Theorem 3). Hence, for $a\in\{\mathrm{adj},\mathrm{bc},\mathrm{cf}\}$, the estimator $\hat{\tau}^{a}$ has the same asymptotic normality as $\mathcal{L}$ if the remainder terms (as well as the bias term $B^{\mathrm{adj}}$ for $\hat{\tau}^{\mathrm{adj}}$) vanish when multiplied by $\sqrt{n}/\sigma_L$. Specifically, the convergence $\sqrt{n}(\hat{\tau}^{a}-\tau)/\sigma_{L}\overset{d}{\to}N(0,1)$ holds if $\kappa p=o(1)$ for $a=\mathrm{adj}$, if $\kappa^2 p=o(1)$ for $a=\mathrm{bc}$, and if $\kappa^2 p^{1/2}+\kappa^3p(\log p)^2=o(1)$ for $a=\mathrm{cf}$, respectively. In the favorable case where all leverage scores are close to their average $p/n$, the regression adjustment estimator and Lei and Ding's bias corrected estimator are asymptotically normal when $p=o(n^{1/2})$ and $p=o(n^{2/3})$, respectively, but the cross-fitted estimator is asymptotically normal when $p=o(n^{3/4}/(\log n)^{1/2})$.}

In addition to the linear component $\mathcal{L}$, Theorem (ref) (ii) characterizes the variance of the second-order quadratic term $\mathcal{W}$. The term $\sigma_{L}^{2}$ is identical to the conventional variance term for the regression adjustment estimator as in Lin (2013). { Therefore, similar to the regression adjustment and Lei and Ding's (2021) bias-corrected estimators, the cross-fitted estimator is guaranteed to be asymptotically more efficient than the difference in means estimator despite arbitrary misspecification of the conditional mean function.} Note that the third component in the expression of $\sigma_{L}^{2}$, $(n-1)^{-1}\sum_{i=1}^{n}(e_{i}(1)-e_{i}(0))^{2}$, has no consistent estimator in general. The additional term $\sigma_{W}^{2}$ also contains a component which cannot be consistently estimated (i.e., the third term of $\sigma_{W}^{2}$).

Compared to the existing results such as \citet*{lei2021regression}, the results on the second-order term $\mathcal{W}$ and its variance $\sigma_{W}^{2}$ are new. Indeed, in their simulation study, Lei and Ding (2021) reported that $\sigma_{L}^{2}$ tends to be lower than the Monte Carlo variance of the point estimator for $\tau$ for larger values of $p$. Based on our higher-order analysis, we argue that this discrepancy can be attributed to the second-order component $\sigma_{W}^{2}$ whose order increases with $p$.

Variance estimation

We next consider variance estimation of the treatment effect estimator, particularly the HC0 and HC3 variance estimators

{ Lei and Ding (2021, Theorem 5) showed that under Assumptions (i)-(iv), both estimators are asymptotically conservative, i.e., $\hat{\sigma}_{\text{HC}j}^{2}/\sigma^{2}_{L}\ge 1-a_{jn}$ with a non-negative sequence $a_{jn}=o_p(1)$ for $j=0$ and $3$.} Under our setup, the properties of these variance estimators are characterized as follows. {

}

This theorem depicts the means of the HC0 and HC3 variance estimators, taking into account of the higher-order terms. First, the first two terms of $\mathbb{E}[\hat{\sigma}_{\text{HC0}}^{2}]$ and $\mathbb{E}[\hat{\sigma}_{\text{HC3}}^{2}]$ are the exact match to the first two terms of $\sigma_{L}^{2}$. However, the third term of $\sigma_{L}^{2}$ is not consistently estimable. Thus, as far as we are concerned with the first-order dominant terms, HC0 and HC3 are conservative estimators of the asymptotic variance of the treatment effect estimators. Second, the third and fourth terms of $\mathbb{E}[\hat{\sigma}_{\text{HC3}}^{2}]$ closely match to the first and second terms of $\sigma_{W}^{2}$, except for the factors $n_{0}/n$ and $n_{1}/n$, respectively. It is interesting to note that the HC3 estimator is interpreted as a jackknife variance estimator. So these multiplicative discrepancies can be understood as emergence of Efron and Stein's (1981) bias for the jackknife variance in higher-order terms in the context of the design-based asymptotic analysis. Third, it should be noted that the signs of the third and fourth terms of $\mathbb{E}[\hat{\sigma}_{\text{HC0}}^{2}]$ are opposite to the corresponding ones in the first and second terms of $\sigma_{W}^{2}$ (or the signs of the third and fourth terms of $\mathbb{E}[\hat{\sigma}_{\text{HC3}}^{2}]$). Therefore, the higher-order term of HC3 slightly overestimates $\sigma_{W}^{2}$, while HC0 severely underestimates $\sigma_{W}^{2}$. This explains relatively poor performances of HC0 in finite samples, as observed in the literature (e.g., simulation studies in \citet*{lei2021regression}).

{ Note that all the terms of $\sigma_W^2$ except the third term are estimable. It is of interest whether one can construct a variance estimator that is guaranteed to be asymptotically conservative for the variance component $\sigma^2_L + \sigma^2_W$ up to the second-order. To this end, we propose a modified version of the HC3 variance estimator:

The asymptotic conservativeness of $\hat{\sigma}_{\text{mHC3 }}^{2}$ is obtained as follows.

} We investigate its finite sample performance in the simulation study below.\footnote{ By estimating the estimable components in $\sigma_W^2$, we can propose alternative bias-corrected versions of the HC0 and HC3 variance estimators defined as follows:

However, analogous arguments to the proof of Corollary (ref) show that these asymptotic variance estimators do not achieve asymptotic conservativeness up to the second-order.}

Stratified randomized experiments with large strata

Stratified randomized experiments using regression adjustment estimators have been considered in liu2020regression under a fixed dimensional asymptotic regime. Our methodology can be extended to the stratified randomized experiments with a finite number of large strata. Consider stratified randomized experiments with $N$ units in the population grouped into strata $s=1,\ldots,S$ for a finite $S$. For each stratum, a randomized experiment is then conducted independently from other strata. The size of the $s$-th stratum is denoted as $n_s\ge 2$. Within the stratum, $n_{1s}$ of them are sampled without replacement and receive treatment while the remaining $n_{0s}=n_{s}-n_{1s}$ are assigned to control. Let $\pi_s=n_{1s}/(n_{1s}+n_{0s})$. Denote the potential outcomes of unit $i$ in stratum $s$ as $(y_{is}(1),y_{is}(0))$, and $T_{is}$, $Y_{is}$, and $z_{is}$ are the corresponding observed outcome, treatment indicator variable, and vector of covariates, respectively. The population average treatment effect is defined as \[ \tau=\sum_{s=1}^S\sum_{i=1}^{n_s} \frac{\tau_{is}}{N}=\sum_{s=1}^S c_s\tau_s, \] where $\tau_{is}=y_{is}(1)-y_{is}(0)$, $\tau_s=\sum_{i=1}^{n_s}\tau_{is}/n_s$, and $c_s=n_s/N$. For each unit $i$ in stratum $s$, define the leave-one-out estimators $\hat \beta_{1s}^{(j)}=\left(\sum_{i\in I_s \setminus \{j\}}T_{is}z_{is}z_{is}^{\prime}\right)^{-1}\left(\sum_{i\in I_s \setminus \{j\}}T_{is}z_{is}Y_{is}\right)$, and $\hat \beta_{0s}^{(j)}=\left(\sum_{i\in I_s \setminus \{j\}}(1-T_{is})z_{is}z_{is}^{\prime}\right)^{-1}\left(\sum_{i\in I_s \setminus \{j\}}(1-T_{is})z_{is}Y_{is}\right)$, where $I_s$ means the units in stratum $s$. The cross-fitted estimator for $\tau$ is then defined as $\hat\tau = \sum_{s=1}^S \sum_{i=1}^{n_s} \hat\tau_s/N$, where $\hat\tau_s=\hat \mu_{1s}^{\text{cf}}-\hat \mu_{0s}^{\text{cf}}$, \[ \hat{\mu}_{1s}^{\mathrm{cf}}=\frac{1}{n_s}\sum_{j\in I_s}\left\{ \frac{T_{js}}{\pi_s}Y_{js}-\left(\frac{T_{js}}{\pi_s}-1\right)(z_{js}^{\prime}\hat{\beta}_{1s}^{(j)})\right\},\quad\hat{\mu}_{0s}^{\mathrm{cf}}=\frac{1}{n_s}\sum_{j\in I_s}\left\{ \frac{1-T_{js}}{1-\pi_s}Y_{js}-\left(\frac{1-T_{js}}{1-\pi_s}-1\right)(z_{js}^{\prime}\hat{\beta}_{0s}^{(j)})\right\}. \] Then, as long as Assumptions (ref) (i)-(iv) hold for each stratum $s$, as $\min_{s=1,\ldots,S} n_s$ diverges to infinity, the conclusions of Theorems (ref)-(ref) continue to hold for each stratum. As the estimators $\hat \tau_1,\ldots,\hat \tau_S$ are mutually independent, the variance of $\hat\tau$ can be estimated by {

} where $\tilde e_{is}$ and $P_{ij,s}$ are defined in the same way as $\tilde e_{i}$, $P_{ij}$, respectively, with observations restricted to the $s$-th stratum. Asymptotically valid inference under this stratified setup can be conducted based on these variance estimators.

Simulation

In this section, we illustrate our theoretical findings through a series of simulation studies. The simulation designs closely follow those in \citet*{lei2021regression}. Specifically, we set the sample size to \( n = 500 \) and define the treatment group size as \( n_{1} = n\pi_{1} \) with \( \pi_{1} = 0.2 \). The covariate matrix \( \mathcal{X} \in \mathbb{R}^{n \times n} \) consists of independent and identically distributed entries drawn from the $t(3)$ distribution. To mimic the design-based asymptotic framework, the matrix \( \mathcal{X} \) is generated once and held fixed across all Monte Carlo replications. Similarly, a vector \( b \in \mathbb{R}^n \), with entries iid from a standard normal distribution, is generated once at the beginning and subsequently kept fixed.

For each dimension \( p \in \{5, 10, \dots, 75\} \), we construct the covariate matrix \( X \in \mathbb{R}^{n \times p} \) by taking the first \( p \) columns of \( \mathcal{X} \), and define parameter vectors \( \beta_{1}^{*} = \beta_{0}^{*} = (b_1, \dots, b_p)^{\prime} \) using the first \( p \) entries of \( b \). The potential outcomes are specified as \[ Y(t) = X \beta_t^* + \epsilon(t), \quad t \in \{0, 1\}, \] where the error vectors \( \epsilon(t) \in \mathbb{R}^n \) are generated under two distinct designs: a worst-case configuration and a normal-error design.

For the worst-case errors, we define \( \epsilon(0) = \epsilon \) and \( \epsilon(1) = 2\epsilon \), where the vector \( \epsilon \in \mathbb{R}^n \) solves the constrained optimization problem: \[ \max_{\epsilon \in \mathbb{R}^n} \left| \frac{n_1}{n_0} \Delta_0 - \frac{n_0}{n_1} \Delta_1 \right| \quad \text{subject to } \frac{1}{n} \epsilon^{\prime} \epsilon = 1, \quad X^{\prime} \epsilon = (1,...,1)^{\prime} \epsilon = 0, \] with \( \Delta_t = n^{-1} \sum_{i=1}^n e_i(t) P_{ii} \). This construction maximizes the first-order bias of the regression-adjusted estimator \( \hat{\tau}^{\mathrm{adj}} \) under the increasing-dimensional asymptotics developed by \citet*{lei2021regression}. For the normal error design, we consider homoskedastic normal errors with \( \epsilon(0) = \epsilon(1) = \epsilon \), where \( \epsilon \sim N(0, I_n) \). In this setting, the potential outcome equations are linear, and biases from regression adjustment are generally small. As with the covariates \( \mathcal{X} \), the error vectors \( \epsilon(t) \) are generated once and fixed throughout the simulation replications. Each simulation design is evaluated with $10,000$ Monte Carlo repetitions.

We compare three estimators for the average treatment effect of interest: (i) the standard regression adjustment estimator based on equation ((ref)) (RA), (ii) the bias-corrected regression adjustment estimator from \citet*{lei2021regression}, given in equation ((ref)) (BC), and (iii) the cross-fitted regression adjustment estimator introduced in this paper, defined in equation ((ref)) (CF). For all estimators, inference is conducted using heteroskedasticity-robust standard errors from the Eicker-Huber-White family, specifically HC2 and HC3. As reported in \citet*{lei2021regression}, HC3 tends to yield the most reliable performance across simulation settings, particularly when the covariate dimension is high. For CF, we further consider inference based on the modified HC3 variance estimator (mHC3), proposed in this paper. The mHC3 estimator is designed to provide improved higher-order accuracy in the context of cross-fitting, potentially offering more robust inference than standard HC3.

Figure (ref) shows the root mean squared errors (RMSE) for the three estimators under normal and worst-case errors. Note that when the errors are normal, all three estimators performed decently. The differences become much more pronounced under the worst-case error structure. RA suffers severe degradation in performance as $p$ increases. Its RMSE rises sharply with $p$, underscoring its vulnerability to adversarial alignment between the error and the design. BC successfully controls this bias, maintaining relatively low RMSE throughout. However, its RMSE still shows mild growth in higher dimensions. CF delivers the best performance, and consistently outperforms RA and BC, and the performance gap widens with increasing $p$, affirming the theoretical advantages of cross-fitting in high-dimensional and misspecified environments.

Figure (ref) shows the coverage rates for the three estimators with different variance estimators. Observe that HC2 fails to maintain nominal coverage as dimension increases. Coverage drops steadily with $p$, indicating that HC2 underestimates the variability of the estimator in challenging settings. HC3 performs better than HC2 but still exhibits noticeable undercoverage for larger values of $p$. Figure (ref) focuses on the CF estimator and highlights the contrast in coverage rates between HC2, HC3, and mHC3. The performance gap between mHC3 and the conventional estimators grows with $p$, showcasing the benefit of incorporating bias correction in the standard error estimator when dealing with highly misspecified or adversarial settings.

{ In sum, we recommend to use the cross-fitted estimator for point estimation of the average treatment effect, and the HC3 or mHC3 variance estimator for interval estimation to practitioners particularly when the number of covariates $p$ is moderate or large. Theoretically, the cross-fitted estimator is asymptotically more efficient than the unbiased difference in means estimator, and admits asymptotic normality under weaker conditions on $p$ compared to the regression adjustment and Lei and Ding's (2021) bias-corrected estimators.}