Assumption-lean covariate adjustment under covariate adaptive randomization when $p = o (n)$

bibunit[apalike] \begin{abstract} Adjusting for (baseline) covariates with working regression models becomes standard practice in the analysis of randomized clinical trials (RCT). When the dimension $p$ of the covariates is large relative to the sample size $n$, specifically $p = o (n)$, adjusting for covariates even in a linear working model by ordinary least squares can yield overly large bias, defeating the purpose of improving efficiency. This issue arises when no structural assumptions are imposed on the outcome model, a scenario that we refer to as the assumption-lean setting. Several new estimators have been proposed to address this issue. However, they focus mainly on simple randomization under the finite-population model, not covering covariate adaptive randomization (CAR) schemes under the superpopulation model. Due to improved covariate balance between treatment groups, CAR is more widely adopted in RCT; and the superpopulation model fits better when subjects are enrolled sequentially or when generalizing to a larger population is of interest. Thus, there is an urgent need to develop procedures in these settings, as the current regulatory guidance provides little concrete direction. In this paper, we fill this gap by demonstrating that an adjusted estimator based on second-order $U$-statistics can almost unbiasedly estimate the average treatment effect and enjoy a guaranteed efficiency gain if $p = o (n)$. In our analysis, we generalize the coupling technique commonly used in the CAR literature to $U$-statistics and also obtain several useful results for analyzing inverse sample Gram matrices by a delicate leave-$m$-out analysis, which may be of independent interest. Both synthetic and semi-synthetic experiments are conducted to demonstrate the superior finite-sample performance of our new estimator compared to popular benchmarks. \end{abstract} {\smallKeywords: covariate adaptive randomization, covariate adjustment, leave-$m$-out analysis, randomized experiments, $U$-statistics} \doublespacing \section{Introduction} Randomization is widely regarded as the gold standard for evaluating treatment effects from clinical trials and interventional studies in other disciplines. Although simple randomization is straightforward to implement in practice and is optimal in theory under certain criteria bai2023randomize, it can nevertheless lead to imbalances in the distributions of (baseline) covariates in different treatment groups and, in turn, to accidental bias and loss of efficiency in estimating treatment effects student1938comparison, efron1971forcing, kasy2016experimenters, banerjee2020theory. As a remedy, randomization schemes that leverage the information in covariates, in particular covariate adaptive randomization (CAR), are now routinely adopted in practice, especially in clinical settings Lin2015. When balancing stratified covariates, stratified block randomization zelen1974randomization and minimization Taves1974,Pocock1975 are commonly deployed in randomized clinical trials (RCT). Stratified block randomization defines sets of strata of a subset of covariates and allocates units within each stratum using block randomization. This method is widely used in the design of RCT, accounting for approximately 70% of all cases ciolino2019ideal, Lin2015. Minimization was initially designed to balance covariates over their marginal distributions Taves1974, Pocock1975, but it has also been generalized to control various types of imbalance measures, including marginal and/or within-stratum imbalances hu2012asymptotic. Other randomization strategies that fall under the purview of CAR include stratified biased coin design efron1971forcing and various model-based approaches Atkinson1982,Begg1980. A rather comprehensive discussion of these randomization schemes in RCT can be found in Rosenberger2015. We also refer those who are interested in the use of CAR in disciplines beyond medicine to duflo2007using. Beyond balancing covariates in the design stage, it is now well accepted that adjusting for covariates in the analysis stage can also improve the estimation efficiency leon2003semiparametric, zhang2008improving, tsiatis2008covariate, ma2024new, van2024covariate. In fact, such practice has been encouraged by two quite impactful regulatory guidelines EMA2015, FDA2023. The robustness of (covariate-)adjusted estimators is highly valued and critical for their acceptance into routine practice because the data-generating process in RCT is typically unknown beyond the treatment assignment mechanism. In general, it is required that the bias of any adjusted estimator, regardless of its source (e.g., model misspecification or curse of dimensionality), must be dominated by sampling variability of order $n^{-1 / 2}$. This is due to the availability of a $\sqrt{n}$-consistent and asymptotically normal ($\sqrt{n}$-CAN) estimator that adjusts for no covariates in the analysis stage, thus imposing no structural assumptions on the outcome model. We now review several existing covariate adjustment methods that are relevant to, but not restricted to, CAR. Consider a trial enrolling $n$ subjects. For simple randomization (without balancing any covariates in the design stage), several studies have explored the properties of regression-adjusted estimators freedman2008regressiona, Freedman2008, lin2013. However, under covariate-adaptive randomization, the robustness of these estimators becomes more challenging to probe due to dependencies in treatment assignments. Bugni2018,Bugni2019 proposed a model-assisted approach that adjusts for stratification covariates. The method proposed there ensures valid statistical inference that is robust to potential misspecification of the outcome model. To adjust for additional covariates and improve efficiency, recent work developed stratum-common and stratum-specific estimators for settings where the covariate dimension $p$ is fixed ma2022regression, ye2023toward, gu2023regression. For high-dimensional regimes ($p \gg n$), liu2022lasso introduced lasso-adjusted estimators that ensure both robustness and efficiency gains. Despite the above progress, there remains an important gap in the theoretical underpinnings of how to construct adjusted estimators of treatment effects, when the dimension $p$ of adjusted covariates is moderately high relative to $n$, or more precisely $p = o (n)$ and $p < n$, without imposing any structural assumptions on the outcome model, such as sparsity as in Bloniarz2016 or low complexity (quantified by the metric entropy integral) as in guo2023generalized. We coin the setting in which no structural assumptions are imposed on the outcome model as the assumption-lean setting, borrowing a terminology recently popularized in the causal inference literature vansteelandt2022assumption. Under simple randomization and the design-based framework (or equivalently, the finite-population model), various bias-corrected adjusted estimators with linear working models have been proposed to reduce bias due to large $p$ lei2021regression, chang2024exact, lu2025debiased, zhao2024hoif, chiang2025regression. As delineated in zhao2024hoif, all the above bias-corrected adjusted estimators can essentially be interpreted as estimators based on a particular type of $U$-statistics motivated by the theory of higher-order influence functions liu2017semiparametric. When it comes to the more general CAR designs that are the main focus of this paper, the relevant literature is scant. We also deviate from the above literature by adopting the superpopulation model to better fit RCT, in which subjects are often enrolled sequentially, and one is often interested in generalizing the conclusions drawn from an RCT to a larger population. To our knowledge, jiang2025adjustments is the only work investigating the statistical properties of adjusting for high-dimensional covariates up to $p = O (n)$. In particular, they demonstrate that the standard OLS adjusted estimator with linear working models ma2022regression, ye2023toward is $\sqrt{n}$-CAN and has a guaranteed efficiency gain compared to the unadjusted estimator as long as $p \ll \sqrt{n}$. However, when $p \gtrsim \sqrt{n}$, jiang2025adjustments assumed a correct linear relationship between the outcome and the covariates to preserve the same statistical guarantee using the OLS adjusted estimator. As will be shown in the sequel, the OLS estimator is not robust against model misspecification once $p \gtrsim \sqrt{n}$. A more detailed comparison of our work with the articles mentioned above is provided in Remark (ref). In summary, it remains elusive whether one can construct a $\sqrt{n}$-CAN adjusted estimator of the average treatment effect (ATE) under the assumption-lean setting and the superpopulation model for CAR, and this research gap is reflected in the following statement in the most recent FDA guideline on covariate adjustment FDA2023: \begingroup \begin{quote} The statistical properties of covariate adjustment are best understood when the number of covariates adjusted for in the study is small relative to the sample size (Tsiatis et al., 2008). Therefore, sponsors should discuss their proposal with the relevant review division if the number of covariates is large relative to the sample size or if proposing to adjust for a covariate with many levels (e.g., study site in a trial with many sites). \end{quote} \endgroup \paragraph{Main contributions and organization} Our main objective is to address the gap targeted by the above quote; we refer the interested readers to Section (ref) for our concrete recommendation. In particular, our paper advances the methodological and theoretical understanding of covariate adjustment on several fronts, while providing insights and guidance on its implementation in practice. Our contributions are summarized as follows. \begin{enumerate}[label = \arabic*)] • Methodologically, we propose a new adjusted estimator of the ATE under CAR and the assumption-lean setting. Our new estimator is essentially a second-order $U$-statistic. In particular, this new estimator is $\sqrt{n}$-CAN, and has guaranteed efficiency gain compared to the unadjusted estimator (Theorem (ref)), as long as $p = o (n)$ in the assumption-lean setting, only under some additional mild tail assumptions on the covariates $X$ and outcome $Y (a)$; see Assumptions (ref)--(ref). Thus, our result resolves an important gap in the literature on covariate adjustment in randomized experiments. Furthermore, a consistent variance estimator of our new adjusted estimator of the ATE is also constructed to facilitate downstream hypothesis testing and statistical inference (Theorem (ref)). • Theoretically, new technical results on the coupling technique adopted in Bugni2018, Bugni2019 for CAR are established to analyze $U$-statistics to prove the statistical properties of our new estimator; see Lemma (ref). Furthermore, compared to the randomization-based framework zhao2024hoif, the proof complexity under the superpopulation model is dramatically escalated due to the randomness in the covariates, demanding a fairly involved leave-two-out analysis of inverse sample Gram matrices; see Lemma (ref). These new technical results could be of independent interest. • Empirically, through extensive simulation studies and a semi-synthetic data analysis, we demonstrate that our new estimator exhibits highly competitive performance in a variety of data generating processes, compared to popular benchmark methods. Based on both theoretical and empirical results, we provide parallel practical guidance that complements the most recent FDA guideline, in particular targeting the paragraph quoted on page (ref). We have incorporated our new adjusted estimator under CAR into an R package available to download from \href{https://cran.r-project.org/web/packages/HOIFCar/index.html}{CRAN} to make the new method ready to be used in practice. \end{enumerate} The remainder of the paper is structured as follows. Section (ref) sets the stage by introducing our statistical problem and the mathematical framework. Both intuitive explanation and numerical evidence are presented in Section (ref) to elucidate why the existing adjusted estimators with linear working models under CAR may have an overly large bias. We also formally introduce our new adjusted estimator at the end of Section (ref). Section (ref) then presents the statistical properties of our new estimator, which are the main theoretical contribution of this paper. In Section (ref), we construct a consistent variance estimator for our new adjusted estimator to facilitate downstream hypothesis testing and valid inference. Synthetic and semi-synthetic data analyses are conducted in Section (ref), to empirically demonstrate the competitive performance of our new adjusted estimator compared to several benchmarks. Section (ref) concludes the paper with a discussion of practical recommendations and some extensions. Technical details are deferred to the Appendix. \section{Framework and Setup} Suppose that we conduct a randomized experiment with $n$ units (enrolled subjects). For each unit $i$, we use $A_i$ to denote the treatment assignment with $A_i=1$ for the treatment group and $A_i=0$ for the control group. In this paper, we only consider two treatment groups, but it is straightforward to extend all our results to more than two groups. We let $n_1 = \sum_{i=1}^n A_i$ and $n_0 = \sum_{i=1}^n (1-A_i)$ be the number of units in the treated group and the control group, respectively. Since we are mainly interested in CAR, we further introduce $\{B_i\}_{i = 1}^{n}$ to denote the strata labels that are balanced in the design stage and each $B_{i}$ takes values in $\mathcal{K} = \{1,2,\dots,K\}$, where $K$ represents the total number of strata, assumed to be fixed throughout. We also collect $p$-dimensional covariates $\{X_i\}_{i = 1}^{n}$ that are not used in the design stage. We let $p_{[k]} = \Pr(B_i=k)$ be the target proportion of stratum $k$ and $\pi_{[k]} = \Pr(A_i=1|B_i=k)$ be the target proportion of treatment group in stratum $k$. We use subscript $[k]$ to indicate statistics calculated in stratum $k$, e.g. $n_{[k]} = |[k]|$, $n_{[k]1} = \sum_{i\in [k]} A_i$ and $n_{[k]0} = \sum_{i\in [k]} (1-A_i)$ representing the numbers of all units, treated units and control units in stratum $k$, respectively. Here, $i \in [k]$ indexes units in stratum $k$. Let $p_{n[k]} = n_{[k]}/n$ and $\pi_{n[k]} = n_{[k]1}/n_{[k]}$ be the proportion of stratum size and the treated group in stratum $k$. We adopt the potential outcomes notation to define treatment effects. Let $Y_i(0)$ and $Y_i(1)$ be the potential outcomes. By the standard consistency assumption, the observed outcome is $Y_i = A_i Y_i(1) + (1-A_i)Y_i(0)$. For any random variable $V$, let $\sigma^2_{[k]V} = \mathrm{var}(V|B = k)$ denote the variance of $V$ in stratum $k$. Let $\Sigma_{[k]} = E_{[k]}\{X X^{\top}\}$ be the population-level Gram matrix of the covariates $X$ in stratum $k$, where $E_{[k]} (\cdot) = E(\cdot|B=k)$ denotes the conditional mean in stratum $k$. We use $\eta_{[k]}(a) = E_{[k]}\{X Y (a)\}$ to denote the inner product between $X$ and $Y (a)$, for $a\in\{0,1\}$. For a generic random variable $r (a)$ related to the treatment group $A = a$, let $\bar{r}_{[k]a} = n_{[k]a}^{-1} \sum_{i\in[k]} \mathbbm{1} \{A_i = a\} r_i(a)$, for $a \in \{0, 1\}$, be the corresponding sample mean in stratum $k$. For example, $r (a)$ may correspond to the potential outcomes $Y (a)$, the covariate $X$, or their product $X \cdot Y (a)$. Also, let $\widehat{\Sigma}_{[k]} = n_{[k]}^{-1} \sum_{i\in[k]}X_i X_i^{\top}$ be the sample Gram matrix of the covariates $X$ in stratum $k$. Let $W_i = \{Y_i(0),Y_i(1),B_i,X_i\}_{i=1}^n$ be independent and identically distributed (i.i.d.) samples from the population $W = \{Y(0),Y(1),B,X\}$. Under CAR, the ATE is identified as \begin{align*} \tau = E \{Y (1) - Y (0)\} = \sum\limits_{k =1}^K p_{[k]} E_{[k]}\{Y (1) - Y (0)\} = \sum\limits_{k=1}^K p_{[k]}\tau_{[k]}, \end{align*} where $\tau_{[k]}$ is the ATE in stratum $k$. Our goal is to construct point and interval estimators of $\tau$ based on the non-i.i.d.\ observed data $\{Y_i,B_i,X_i,A_i\}_{i=1}^n$, where the dependency is induced by CAR. Throughout, we impose the following on the treatment assignment mechanism. \begin{assumption}\ \begin{enumerate} • Conditional on $\{B_1,\dots,B_n\}$, $\{A_1,\dots,A_n\}$ is independent of $\{Y_i(0),Y_i(1),X_i\}_{i=1}^n$. • $\pi_{n[k]} \stackrel{P}\to \pi_{[k]}$ as $n \to \infty$, for all $k \in \mathcal{K}$. \end{enumerate} \end{assumption} Assumption (ref) is also imposed in Bugni2019, ma2022regression and liu2022lasso, which is satisfied by several commonly used CAR procedures, including simple randomization, stratified block randomization zelen1974randomization and stratified biased coin randomization efron1971forcing. When $\pi_{[k]}$'s are the same across different strata, Pocock and Simon's minimization method Pocock1975 also satisfies this assumption Ma2015. Before formally presenting our new covariate-adjusted ATE estimator (abbreviated as adjusted estimator for short), we first introduce the stratified difference-in-means estimator (referred to as the \emph{unadjusted estimator} here, as it does not adjust for additional covariates $X_i$), which is, to the best of our knowledge, first proposed in Bugni2019: \begin{equation} \widehat\tau_{\rm unadj} = \sum\limits_{k=1}^K p_{n[k]} \big(\bar Y_{[k]1} -\bar Y_{[k]0} \big). \end{equation} This unadjusted estimator $\widehat{\tau}_{\mathrm{unadj}}$ serves as a benchmark against which both existing adjusted estimators and our newly proposed adjusted estimator are compared. Importantly, Bugni2019 showed that $\widehat{\tau}_{\mathrm{unadj}}$ $\sqrt{n}$-CAN, without any assumptions on the outcome model. Despite the advantages mentioned above, $\widehat{\tau}_{\mathrm{unadj}}$ does not fully incorporate all available information from the data, therefore resulting in a loss of efficiency. We next introduce the OLS estimator $\widehat\tau_{\rm OLS}$ proposed by ma2022regression, which is commonly used in practice to improve efficiency compared to the unadjusted estimator when the covariate dimension $p$ is fixed, defined as \begin{equation} \begin{split} \widehat\tau_{\rm OLS} = & \ \sum\limits_{k=1}^K p_{n[k]} \Big( \Big\{\bar Y_{[k]1} - \dfrac{1}{n_{[k]1}}\sum\limits_{i\in[k]}(A_i-\pi_{n[k]})X_i^{\top}\widehat\beta_{[k]}(1) \Big\} \\ & - \Big\{\bar Y_{[k]0} - \dfrac{1}{n_{[k]0}}\sum\limits_{i\in[k]}(\pi_{n[k]}-A_i)X_i^{\top}\widehat\beta_{[k]}(0) \Big\} \Big), \end{split} \end{equation} where \[\widehat\beta_{[k]}(a) = \widehat{\Sigma}_{[k]}^{-1} \cdot \frac{1}{n_{[k]a}}\sum\limits_{i\in[k]} \mathbbm{1}\{A_i=a\}X_iY_i,\quad a\in\{0,1\}.\] A practical advantage of $\widehat{\tau}_{\rm OLS}$ is that it can be computed by straightforward post-processing of the output of standard, off-the-shelf statistical software packages. To see this, the OLS estimator can in fact be obtained by first regressing $Y$ against $A$, $X$ and their interactions within each stratum and then taking a weighted sum of the regression coefficients of $A$ across strata, with weights $p_{n[k]}$. Specifically, $\widehat\beta_{[k]}(a)$ corresponds to the coefficient of the interaction term between $A$ and $X$ in each stratum. We now recall the statistical properties of $\widehat\tau_{\rm OLS}$ in Proposition (ref) below, which is adapted from ma2022regression. \begin{proposition} Under Assumption (ref) and $E\{Y^2_i(a)\}<\infty$, let $\beta_{[k]}(a) = \Sigma_{[k]}^{-1} \eta_{[k]}(a)$, $\beta_{[k]} = (1-\pi_{n[k]})\beta_{[k]}(1) + \pi_{n[k]}\beta_{[k]}(0)$ and $r_i(a) = Y_i(a) - X_i^{\top}\beta_{[k]}$. Define $\sigma^2_{\rm OLS} = \zeta^2_{H} + \zeta^2_{\mathrm{I},r}(\pi_{[k]})$. When $p$ is fixed, we have \[\sqrt{n} (\widehat\tau_{\rm OLS} - \tau)/\sigma_{\rm OLS} \stackrel{d}\to \operatorname{N}(0,1),\] where \begin{equation} \zeta^2_{\mathrm{I},r}(\pi_{[k]}) = \sum\limits_{k=1}^K p_{[k]}\Big(\frac{\sigma^2_{[k]r(1)}}{\pi_{[k]}} + \frac{\sigma^2_{[k]r(0)}}{1-\pi_{[k]}}\Big),\ \zeta^2_{H} = E\{\check{r}_i(1) - \check{r}_i(0)\}^2, \end{equation} and $\check{r}_i(a) = E_{[k]}\{r_i(a)\} - E\{r_i(a)\}$ is the difference between the mean of $r (a)$ in stratum $k$ and the mean of $r (a)$ in the population. \end{proposition} Before moving forward, we unpack the seemingly cumbersome notation $\zeta_{\rm I, r}^{2} (\pi_{[k]})$ that appeared in Proposition (ref) above. The subscript $\mathrm{I}$ is used to represent sample means, distinguished from second-order $U$-statistics (demarcated by subscript $\mathrm{II}$) to appear in the next section. The other subscript $r$ indicates that $\zeta_{\rm I, r}^{2} (\pi_{[k]})$ corresponds to the asymptotic variance of the transformed outcome $r$. At this stage, it may seem that we are overloading the notation, but it will facilitate the presentation in Section (ref), when we discuss how to estimate the variance of our new estimator. \begin{remark} In this paper, we simply take $p$ as the dimension of all the covariates. In principle, we allow the dimension of all the covariates to be different from that of the covariates being adjusted. As we will see in Section (ref), $p = o (n)$ suffices for our proposed estimator to be asymptotic normal with asymptotic variance no greater than that of the unadjusted estimator. Therefore, a natural next step of our work is to develop estimators that can incorporate data-adaptive variable selection methods van2024automated that adaptively choose $o (n)$ many variables if the total number of covariates is greater than $n$. \end{remark} \section{Bias-Corrected Estimator Based on \texorpdfstring{$U$}-Statistics} In the previous section, in particular in Proposition (ref), we have seen that the OLS estimator $\widehat\tau_{\rm OLS}$ can be more efficient than the unadjusted estimator $\widehat{\tau}_{\mathrm{unadj}}$, when the dimension $p$ of adjusted covariates is fixed or, more generally, when $p \ll \sqrt{n}$. However, in the assumption-lean setting, even when $p$ is moderately high in the sense that $p = o (n)$, the conclusions of Proposition (ref) no longer hold and $\widehat{\tau}_{\rm OLS}$ can suffer from non-negligible bias. If the bias of $\widehat{\tau}_{\rm OLS}$ exceeds or even just equals $n^{-1 / 2}$ in order, it defeats the purpose of improving estimation efficiency by adjusting for covariates, because downstream hypothesis testing or statistical inference based on standard Wald confidence intervals will be invalid. We now explain the limitations of $\widehat\tau_{\rm OLS}$ when $p$ is close to $n$. Recall that the OLS estimator takes the form in (ref). Suppose that the population-level Gram matrix $\Sigma_{[k]}$ is known and we temporarily replace $\widehat{\Sigma}_{[k]}$ in $\widehat{\beta}_{[k]} (a)$ by $\Sigma_{[k]}$, which appeared in $\widehat\tau_{\rm OLS}$ defined in (ref). Using the treatment group as an example, by spelling out $\widehat{\beta}_{[k]} (a)$ for $a \in \{0, 1\}$, we can rewrite the linear adjustment term (or sometimes interpreted as the augmentation term) of $\widehat{\tau}_{\rm OLS}$ in stratum $k$ as a second-order $V$-statistic \[\frac{1}{n^2_{[k]1}}\sum\limits_{i\in[k]}\sum\limits_{j\in [k]}(A_i-\pi_{n[k]})X_i^{\top}\Sigma_{[k]}^{-1} A_j X_jY_j.\] Based on this $V$-statistic representation, it is not difficult to see that it introduces a bias of order $p / n$ through the diagonal component in the double summation: \begin{equation} \frac{1}{n_{[k]1}^2}\sum\limits_{i\in[k]}(A_i-\pi_{n[k]})X_i^{\top}\Sigma_{[k]}^{-1} A_i X_iY_i = \frac{1-\pi_{n[k]}}{n_{[k]1}^2}\sum\limits_{i\in[k]}A_iX_i^{\top}\Sigma_{[k]}^{-1}X_iY_i =O_P\left(\frac{p}{n}\right), \end{equation} where the last equality follows because each summand $|A_iX_i^{\top}\Sigma_{[k]}^{-1}X_iY_i| \lesssim \Vert X_i \Vert_{2}^{2} \Vert \Sigma_{[k]} \Vert_{\rm op} = O_{P} (p)$, if we assume that each dimension of $X_{i}$ is of order $O_{P} (1)$ and $\Sigma_{[k]}$ has bounded eigenvalues (see Assumptions (ref) and (ref) in Section (ref)). In Figure (ref), we numerically illustrate how the bias of $\widehat\tau_{\rm OLS}$, analytically characterized in (ref), scales with $p$. For comparison, we also report the bias of an oracle bias-corrected estimator to be introduced later, based on the oracle knowledge of $\Sigma_{[k]}$. The simulated data are generated according to a similar setting as Model 1 in Section (ref) later and we defer the details to Appendix (ref). In Figure (ref)(a), we gather both the (absolute) empirical bias of $\widehat{\tau}_{\rm OLS}$ calculated by taking the average of $\widehat{\tau}_{\rm OLS} - \tau$ over repeated Monte Carlo draws and the analytical bias calculated using (ref), confirming that the absolute bias of $\widehat{\tau}_{\rm OLS}$ indeed increases substantially with $p$ and tracks the analytic formula (ref) closely. In contrast, the oracle bias-corrected estimator is nearly unbiased, with an absolute bias near zero across all $p$'s considered. We then display in Figure (ref)(b) the histograms of both $\widehat{\tau}_{\rm OLS}$ and the oracle bias-corrected estimator respectively at $p=5$ and $p=60$. We rescaled both estimators by subtracting the true effect and dividing their corresponding Monte Carlo standard deviations. When $p=5$, the two estimators exhibit similar performance due to the small bias. However, when $p=60$, the OLS estimator exhibits a substantial bias, whereas the oracle bias-corrected estimator remains nearly unbiased. \begin{figure}[H] \caption{(a) Absolute bias of the estimators relative to the true effect as the covariate dimension $p$ increases, averaged over 2000 Monte Carlo replicates. The covariate dimension refers to the set of variables adjusted in the estimators. The oracle estimator $\widehat{\tau}_{\mathrm{ora}}$ will be introduced later in this section. The vertical dotted lines mark the $p$'s at which the histograms of estimators are displayed in panel (b). (b) Histograms of the scaled estimators for two representative cases, $p=5$ and $p=60$, where scaling is performed by dividing each estimator by its Monte Carlo standard deviation computed from 2000 replicates.} \end{figure} The above observation prompts the need to develop new estimators that have negligible bias compared to sampling variability when $p$ is large compared to $n$ under the assumption-lean setting. To introduce our new estimator, we start by still assuming the oracle knowledge of $\Sigma_{[k]}$. We then reveal the oracle estimator reported in Figure (ref), which is a second-order $U$-statistic and removes the bias of $\widehat{\tau}_{\rm OLS}$. First, for each stratum $k$, we compute two separate second-order $U$-statistics, respectively, in the treatment and control groups: \begin{equation} \begin{split} & \mathbf{U}_{n_{[k]},2}(\Sigma_{[k]}^{-1}; 1) = \frac{1}{n_{[k]}(n_{[k]}-1)}\frac{1}{\pi^2_{n[k]}} \sum_{\substack{1\leq i\neq j \leq n\\ i,j\in[k] }} (A_i - \pi_{n[k]})X_i^{\top} \Sigma_{[k]}^{-1} A_{j} X_j Y_j, \\ & \mathbf{U}_{n_{[k]},2}(\Sigma_{[k]}^{-1};0) = \frac{1}{n_{[k]}(n_{[k]}-1)}\frac{1}{(1-\pi_{n[k]})^2} \sum_{\substack{1\leq i\neq j \leq n\\ i,j\in[k] }} (\pi_{n[k]}-A_i)X_i^{\top} \Sigma_{[k]}^{-1} (1-A_j) X_j Y_j. \end{split} \end{equation} Here, we introduce the short-hand notation $\mathbf{U}_{n_{[k]},2}(\Sigma_{[k]}^{-1}; 1)$ and $\mathbf{U}_{n_{[k]},2}(\Sigma_{[k]}^{-1}; 0)$ for the ease of exposition. The oracle bias-corrected adjusted estimator of $\tau$ then reads as: \[\widehat\tau_{\mathrm{ora}} = \sum\limits_{k=1}^K p_{n[k]}\Big[\Big\{\bar Y_{[k]1} - \mathbf{U}_{n_{[k]},2}(\Sigma_{[k]};1) \Big\}-\Big\{\bar Y_{[k]0} - \mathbf{U}_{n_{[k]},2}(\Sigma_{[k]};0) \Big\}\Big].\] Viewing $\widehat{\tau}_{\rm OLS}$ as a $V$-statistic, by removing the diagonal components using a $U$-statistic, the oracle adjusted estimator reduces the bias due to large $p$. In reality, the population-level Gram matrix $\Sigma_{[k]}$ is in general unknown and needs to be estimated from the data. A natural approach is to use the sample Gram matrix estimator $\widehat{\Sigma}_{[k]}$ and plug it in. Replacing $\Sigma_{[k]}$ with $\widehat{\Sigma}_{[k]}$, we finally obtain the feasible adjusted estimator $\widehat\tau$ defined as: \begin{equation} \widehat\tau = \sum\limits_{k=1}^K p_{n[k]}\Big[\Big\{\bar Y_{[k]1} - \mathbf{U}_{n_{[k]},2}(\widehat{\Sigma}_{[k]};1) \Big\}-\Big\{\bar Y_{[k]0} - \mathbf{U}_{n_{[k]},2}(\widehat{\Sigma}_{[k]};0) \Big\}\Big]. \end{equation} The new adjusted estimator $\widehat{\tau}$ constitutes one of the main methodological contributions of our paper. As will be demonstrated in Section (ref), the bias of $\widehat{\tau}$ is negligible compared to its standard deviation, and it is guaranteed that the asymptotic variance of $\widehat{\tau}$ never exceeds that of $\widehat{\tau}_{\mathrm{unadj}}$ under the assumption-lean setting, as long as $p = o (n)$. We refer readers to Remark (ref) for further discussions on the bias of $\widehat{\tau}$ when $p > n$. \begin{remark} If treatments are independently assigned within each stratum, the oracle adjusted estimator is exactly unbiased, i.e., $E\{\widehat\tau_{\mathrm{ora}} - \tau\}=0$. In contrast, in CAR, treatment assignments $A_i$ and $A_j$ within the same block may be dependent, but, as we will show later in Section (ref), $\widehat\tau_{\mathrm{ora}}$ is still asymptotically unbiased even after being scaled by $\sqrt{n}$: \[\lim_{n \rightarrow \infty} E\{\sqrt{n} (\widehat\tau_{\mathrm{ora}}-\tau)\}\to 0,\] It is not difficult to completely remove the remaining bias of order $o (n^{-1 / 2})$, which, however, is a less important issue and will not be further examined in this paper. \end{remark} \begin{remark} It is noteworthy that our new estimator can also be interpreted as a modified OLS estimator simply by replacing the regression coefficients $\widehat{\beta}_{[k]} (a)$ by their leave-one-out (LOO) forms. To see this, define the oracle LOO regression coefficients as: for any $i = 1, \cdots, n$, \begin{align*} \widetilde{\beta}^{-i}_{[k]}(1) = \Sigma_{[k]}^{-1} \cdot \frac{1}{n_{[k]}-1} \sum_{j\neq i} \frac{A_j}{\pi_{n[k]}} X_j Y_j, \, \, \, \, \text{and} \, \, \, \, \widetilde{\beta}^{-i}_{[k]}(0) = \Sigma_{[k]}^{-1} \cdot \frac{1}{n_{[k]}-1} \sum_{j\neq i} \frac{1-A_j}{1-\pi_{n[k]}} X_j Y_j. \end{align*} We then obtain that, for each $a \in \{0, 1\}$, \[\mathbf{U}_{n_{[k]},2}(\Sigma_{[k]}^{-1};a) = \frac{1}{n_{[k]a}} \sum\limits_{i\in [k]}(\mathbbm{1} \{A_i=a\} - \pi_{n[k]})X_i^{\top} \widetilde{\beta}^{-i}_{[k]}(a).\] We in turn define the oracle adjusted estimator $\widehat{\tau}_{\rm ora}$ based on second-order $U$-statistics: \begin{equation*} \widehat\tau_{\mathrm{ora}} = \ \sum\limits_{k=1}^K p_{n[k]}\Big[\Big\{\bar Y_{[k]1} - \frac{1}{n_{[k]1}}\sum\limits_{i\in[k]}(A_i-\pi_{n[k]})X_i^{\top} \widetilde{\beta}^{-i}_{[k]}(1) \Big\}-\Big\{\bar Y_{[k]0} - \frac{1}{n_{[k]0}}\sum\limits_{i\in[k]}(\pi_{n[k]}-A_i)X_i^{\top} \widetilde{\beta}^{-i}_{[k]}(0) \Big\}\Big]. \end{equation*} Compared to $\widehat\tau_{\rm OLS}$, this LOO formulation of $\widehat{\tau}_{\mathrm{ora}}$ replaces the OLS regression coefficients $\widehat{\beta}_{[k]}(a)$ by the oracle LOO coefficients $\widetilde{\beta}^{-i}_{[k]}(a)$ for each $a\in\{0,1\}$. \end{remark} \section{Statistical Properties of the Bias-Corrected Estimator} In this section, for ease of our exposition, we first present the statistical properties of the oracle estimator $\widehat\tau_{\mathrm{ora}}$, before moving on to the properties of the feasible estimator $\widehat\tau$, which are the main theoretical results of this paper. The oracle estimator $\widehat\tau_{\mathrm{ora}}$ serves as an ideal technical device that bridges the feasible estimator $\widehat{\tau}$ and the true ATE $\tau$. Unlike $\widehat{\tau}$, the $U$-statistic kernel of $\widehat\tau_{\mathrm{ora}}$ does not depend on the entire sample through $\widehat{\Sigma}_{[k]}$, so it is relatively straightforward to establish its bias, variance, and asymptotic normality. Theoretical results of $\widehat{\tau}_{\mathrm{ora}}$ then conceptually simplify the analysis of the feasible estimator $\widehat{\tau}$: Once the results are in place in the oracle setting, deriving the properties of $\widehat{\tau}$ reduces to controlling the difference between $\widehat{\Sigma}_{[k]}$ and $\Sigma_{[k]}$. However, it should be noted that this final step turns out to be technically challenging and involves a delicate decoupling leave-out analysis of $\widehat{\Sigma}_{[k]}^{-1}$. \subsection{The oracle estimator} Before stating our theoretical results, we further introduce some mild regularity conditions required on the distributions of $X_i$ and $Y_i (a)$. \begin{assumption}[Distributions of $X$ and $Y (a)$] There exists an absolute constant $M > 0$ such that the following hold: \begin{enumerate}[label = (\arabic*)] • $X_i$ is uniformly bounded by $M$ in the sense that there exists constant $M$ such that $\max_{i,j} |X_{ij}| \leq M < \infty$, where $X_{ij}$ is the $j$-th covariate of the $i$-th unit. • Given $X_i$, the moments of $Y_i(a)$ up to order 4 are all bounded by $M$ almost surely: $\max_{i} E_{[k]} [|Y_i(a)|^4|X_i] \leq M < \infty$ almost surely for $a \in\{0,1\},\ k\in\mathcal{K}$. \end{enumerate} \end{assumption} \begin{assumption}[Eigenvalues of $\Sigma_{[k]}$] Let $\Lambda_{\min}$ and $\Lambda_{\max}$ be, respectively, the minimum and maximum eigenvalues of the stratum-specific population Gram matrix $\Sigma_{[k]}$. There exist two absolute constants $\kappa_l$ and $\kappa_{u}$ independent of $n$ such that $0 < \kappa_l \leq \Lambda_{\min} \leq \Lambda_{\max} \leq \kappa_u < \infty$. \end{assumption} Assumption (ref) imposes tail conditions on the covariates $X$ and the potential outcome $Y (a)$, which are mild in the context of RCT analysis. Specifically, it is standard to rescale $X$ in practice so that $X$ can be viewed as bounded. The bounded fourth-moment condition on $Y (a)$ is imposed in Freedman2008 and lin2013. We conjecture that it is possible to further relax Assumption (ref)(1) from bounded $X$ to sub-Gaussian tails, but proving this is beyond the scope of this paper. Assumption (ref) incurs almost no loss of generality since we only assume that the population Gram matrix has bounded eigenvalues from above and below. These two assumptions are needed to upper bound the variance and prove the asymptotic normality of $\widehat{\tau}_{\mathrm{ora}}$ or $\widehat{\tau}$. Before presenting the main theoretical result of this section, we first introduce the following generic notation that will be frequently encountered in the rest of our paper: for $a, b \in \{0, 1\}$ and $l, m \in \{1, 2\}$, \begin{equation} \zeta^2_{\mathrm{II},r[l,m]}(a,b,w(\pi_{[k]})) = \sum\limits_{k=1}^K w(\pi_{[k]})\frac{p_{[k]}}{n_{[k]}-1}E_{[k]}\{(X_1^{\top}\Sigma_{[k]}^{-1}X_2)^2 r_l(a)r_m(b)\}, \end{equation} where $w: [0, 1] \to \mathbb{R}$ is a function of the stratum-specific treatment-assignment probability $\pi_{[k]}$. Here, the new notation $\zeta^2_{\mathrm{II},r[l,m]}(a,b,w(\pi_{[k]}))$ resembles the notation $\zeta^2_{\mathrm{I},r}(\pi_{[k]})$ introduced in (ref) previously. Armed with (ref), we can introduce the following quantities that are useful in representing the variance of $\widehat{\tau}_{\rm ora}$. To this end, define \begin{equation*} \sigma^2 = \zeta^2_{H} + \zeta^2_{\mathrm{I},r}(\pi_{[k]}) + \zeta^2_{\mathrm{II}}, \end{equation*} where $\zeta^2_{\mathrm{II}} = \zeta^2_{\mathrm{II},Y(1)} + \zeta^2_{\mathrm{II},Y(0)} - 2 \zeta^2_{\mathrm{II},Y[1,2]}(1,0,1)$, and \begin{align*} \zeta^2_{\mathrm{II},Y(1)} &= \zeta^2_{\mathrm{II},Y[1,1]}\Big(1,1,\frac{1-\pi_{[k]}}{\pi^2_{[k]}}\Big)+ \zeta^2_{\mathrm{II},Y[1,2]}\Big(1,1,\frac{(1-\pi_{[k]})^2}{\pi^2_{[k]}}\Big),\\ \zeta^2_{\mathrm{II},Y(0)} &= \zeta^2_{\mathrm{II},Y[1,1]}\Big(0,0,\frac{\pi_{[k]}}{(1-\pi_{[k]})^2}\Big) + \zeta^2_{\mathrm{II},Y[1,2]}\Big(0,0,\frac{\pi^2_{[k]}}{(1-\pi_{[k]})^2}\Big). \end{align*} We are now ready to present the statistical properties of $\widehat\tau_{\mathrm{ora}}$ in Theorem (ref) below. \begin{theorem} Under Assumptions (ref)--(ref), when $p \lesssim n$, the following hold. \begin{enumerate}[label = (\alph*)] • The bias and variance of the oracle estimator $\widehat\tau_{\mathrm{ora}}$ satisfy: \begin{align*} & E\{\sqrt n (\widehat\tau_{\mathrm{ora}} - \tau)\} = o(1), \text{ and} \\ & \mathrm{var}\{\sqrt n (\widehat\tau_{\mathrm{ora}} - \tau)\} = \sigma^2 + o (1). \end{align*} • Furthermore, $\widehat\tau_{\mathrm{ora}}$ is $\sqrt{n}$-CAN, or more precisely: \[\sqrt{n}(\widehat\tau_{\mathrm{ora}} - \tau)/\sigma\stackrel{d}\to \operatorname{N}(0,1).\] \end{enumerate} \end{theorem} Since $\widehat{\tau}_{\mathrm{ora}}$ is a $U$-statistic, the proof of Theorem (ref) generalizes the coupling technique of Bugni2019 to $U$-statistics, which could be of independent interest; see Appendix (ref). In Theorem (ref), we only need $p \lesssim n$ for the statements to hold because $\Sigma_{[k]}$ is known and we do not need $p < n$ to ensure $\Sigma_{[k]}$ to be invertible (guaranteed by Assumption (ref)). In fact, the asymptotic normality of $\widehat{\tau}_{\mathrm{ora}}$ is maintained as long as $p = o (n^{2})$, if we scale $\widehat{\tau}_{\mathrm{ora}} - \tau$ by $\max \{\sqrt{n}, \sqrt{n^{2} / p}\}$ instead of $\sqrt{n}$ liu2020nearly. It is well known in the literature that $\sigma_{\rm OLS}^{2}$ is guaranteed to be smaller than or equal to the asymptotic variance of $\sqrt{n} \widehat{\tau}_{\mathrm{unadj}}$. In the proposition below, we also establish the relationship between $\sigma_{\rm OLS}^{2}$ and $\sigma^{2}$. In particular, for $\sigma^{2}$ to be of the same order as $\sigma_{\rm OLS}^{2}$, $p = o (n)$ is still needed. Therefore, if $p = o (n)$, $\widehat{\tau}_{\mathrm{ora}}$ is also guaranteed to be never less efficient than $\widehat{\tau}_{\mathrm{unadj}}$. \begin{proposition} The following hold. \begin{enumerate}[label = (\alph*)] • $\zeta^{2}_{\rm II} \geq 0$, or equivalently $\sigma_{\rm OLS}^{2} \leq \sigma^{2}$; • Under the assumptions of Theorem (ref), if $p = o (n)$, $\lim_{n \rightarrow \infty} \sigma^{2} / \sigma^{2}_{\rm OLS} = 1$. \end{enumerate} \end{proposition} \begin{remark} ma2022regression have shown that when $p$ is fixed, $\widehat{\tau}_{\mathrm{OLS}}$ has a guaranteed efficiency gain compared to $\widehat{\tau}_{\mathrm{unadj}}$. When $p = o(n)$, our oracle estimator $\widehat{\tau}_{\mathrm{ora}}$ achieves the same asymptotic variance as $\widehat{\tau}_{\mathrm{OLS}}$, and is thus also more efficient than or as efficient as $\widehat{\tau}_{\mathrm{unadj}}$. The extra term $\zeta^2_{\mathrm{II}}$ in $\mathrm{var} \{\widehat{\tau}_{\mathrm{ora}}\}$ compared to $\mathrm{var} \{\widehat{\tau}_{\rm OLS}\}$ is $O (p/n)$, and it can be shown that $\zeta^2_{\mathrm{II}}\geq0$ (see Appendix (ref)). As a result, when $p/n \to \gamma\in(0,1]$, the asymptotic variance of $\widehat\tau_{\mathrm{ora}}$ may exceed that of $\widehat\tau_{\mathrm{unadj}}$ in the proportional asymptotic regime, potentially reducing efficiency unless further conditions are assumed. See Section (ref) for some further discussions. \end{remark} \begin{remark} Recent work has also focused on covariate adjustment in randomized experiments in which $p$ diverges with $n$. Most closely related to our work, zhao2024hoif employed second-order $U$-statistics to estimate linear adjustments under the finite-population model with simple randomization. Separately, chang2024exact and lu2025debiased derived a debiased regression-adjusted estimator within the same framework; zhao2024hoif noted that these estimators are asymptotically equivalent and had been introduced in almost identical forms in observational studies liu2020nearly, liu2020rejoinder, liu2023new using the framework of higher-order influence functions liu2017semiparametric, leading to the same asymptotic normality. We focus on CAR under a superpopulation model, which differs from the finite-population model in both theoretical assumptions and randomization schemes. While jiang2025adjustments also consider the regime in which $\sqrt{n} \lesssim p \lesssim n$ in the superpopulation model, they assume a correct linear relationship between potential outcomes and covariates. Under this predicate, the OLS estimator $\widehat{\tau}_{\rm OLS}$, although a $V$-statistic instead of a $U$-statistic, is still $\sqrt{n}$-CAN. In contrast, Theorem (ref) allows for possible misspecification of the linear model. As we have illustrated in Section (ref), the regression-adjusted estimator based on $V$-statistics may have bias diverging to infinity after being standardized by $\sqrt{n}$ under model misspecification. \end{remark} \subsection{The feasible estimator} Theorem (ref) establishes the asymptotic normality of the oracle estimator $\widehat\tau_{\mathrm{ora}}$, but this theoretical result is derived under the often unrealistic assumption that the stratum-specific population Gram matrix $\Sigma_{[k]}$ is known. In practice, $\Sigma_{[k]}$ shall be estimated from the observed data. As indicated at the end of Section (ref), the feasible estimator $\widehat{\tau}$ simply replaces $\Sigma_{[k]}$ by the corresponding sample covariance estimator $\widehat{\Sigma}_{[k]}$. Recall from (ref) in Section (ref) that \[\widehat\tau = \sum\limits_{k=1}^K p_{n[k]}\Big[\Big\{\bar Y_{[k]1} - \mathbf{U}_{n_{[k]},2}(\widehat{\Sigma}_{[k]};1) \Big\}-\Big\{\bar Y_{[k]0} - \mathbf{U}_{n_{[k]},2}(\widehat{\Sigma}_{[k]};0) \Big\}\Big].\] We first present the most important theoretical result of this paper, concerning the asymptotic statistical properties of $\widehat{\tau}$. \begin{theorem} Under Assumptions (ref)--(ref), the following hold. \begin{itemize} • When $p<n$, the feasible estimator $\widehat\tau$ is asymptotic unbiased in the sense that \begin{align*} E\{\sqrt n (\widehat\tau - \tau)\} = o(1). \end{align*} \end{itemize} We then also assume $p = o (n)$. \begin{itemize} • The asymptotic variance of $\widehat\tau$ satisfies \[\mathrm{var}\{\sqrt{n} (\widehat\tau - \tau)\} = \sigma^2 + o (1).\] • Moreover, $\widehat\tau$ is $\sqrt n$-CAN in the sense that \[\sqrt{n}(\widehat\tau - \tau)/\sigma \stackrel{d}\to \operatorname{N}(0,1).\] In view of Proposition (ref), when $p = o (n)$, $\lim_{n \rightarrow \infty} \sigma^{2} / \sigma^{2}_{\rm OLS} = 1$, we also have \[\sqrt{n}(\widehat\tau - \tau)/\sigma_{\rm OLS} \stackrel{d}\to \operatorname{N}(0,1).\] \end{itemize} \end{theorem} In fact, to maintain the same statistical properties of $\widehat{\tau}$ as in Theorem (ref), in principle $\widehat{\Sigma}_{[k]}^{-1}$ can be replaced by any generic estimator $\widetilde{\Sigma}_{[k]}^{-1}$ of $\Sigma_{[k]}^{-1}$, as long as $\widetilde{\Sigma}_{[k]}$ satisfies the following: \begin{equation} \sqrt{n}\Big\{\mathbf{U}_{n_{[k]},2}(\widetilde{\Sigma}_{[k]}^{-1} - \Sigma_{[k]}^{-1};1) - \mathbf{U}_{n_{[k]},2}(\widetilde{\Sigma}_{[k]}^{-1} - \Sigma_{[k]}^{-1};0)\Big\} = o_P(1), \end{equation} for all $k \in \mathcal{K}$. Here, for $a = 0, 1$, the notation $\mathbf{U}_{n_{[k]},2}(\widetilde{\Sigma}_{[k]}^{-1} - \Sigma_{[k]}^{-1};a)$ is the same as $\mathbf{U}_{n_{[k]},2}(\Sigma_{[k]}^{-1};a)$ defined in (ref), replacing $\Sigma_{[k]}^{-1}$ by $\widetilde{\Sigma}_{[k]}^{-1} - \Sigma_{[k]}^{-1}$. Condition (ref) characterizes the relevant functional convergence rate of $\widetilde{\Sigma}_{[k]}^{-1}$ to $\Sigma_{[k]}^{-1}$ such that the difference between $\widehat\tau_{\mathrm{ora}}$ and $\widehat\tau$ is $o_P (n^{-1/2})$. In the proof, we show that the sample Gram matrix estimator $\widehat{\Sigma}_{[k]}$ satisfies Condition (ref) when $p = o(n)$. However, it is possible to consider other estimators of $\Sigma_{[k]}^{-1}$, such as the ridge inverse Gram matrix estimator liu2025augmented, abadie2025unbiased. We leave the justification whether other estimators of $\Sigma_{[k]}^{-1}$ also satisfy Condition (ref) to future work. In the proof of Theorem (ref) (see Appendix (ref)), the use of the sample Gram matrix $\widehat{\Sigma}_{[k]}$ introduces two technical challenges. First, since $\widehat{\Sigma}_{[k]}$ involves all the data in stratum $k$, it induces complicated dependencies into the corresponding $U$-statistic kernel. A standard kernel $h$ of a second-order $U$-statistic is data-independent, while the corresponding kernel of $\widehat{\tau}$ depends on all $X_{i}$'s in stratum $k$ through $\widehat{\Sigma}_{[k]}$. To decouple such a complicated dependency structure, we obtain a technical lemma (Lemma (ref)) that repeatedly applies the Sherman–Morrison formula to carry out a delicate leave-two-out analysis of $\widehat{\Sigma}_{[k]}^{-1}$. This result could be of independent interest in other problems dealing with inverse sample Gram matrices with large $p$ liu2017semiparametric, bao2025leave. The second challenge lies in characterizing the asymptotic variance of $\widehat{\tau}$. To show that $\widehat{\tau}$ has the same asymptotic variance $\sigma^{2}$ as $\widehat\tau_{\mathrm{ora}}$, a sufficient condition is $\lVert\widehat{\Sigma}_{[k]}^{-1} - \Sigma_{[k]}^{-1}\rVert_2 = o_{P} (1)$, which holds automatically if $p = o (n)$ by standard matrix concentration bounds. When $p = O (n)$, the bias of $\widehat{\tau}$ is still $o (n^{-1 / 2})$ but the asymptotic variance of $\widehat{\tau}$ is generally not $\sigma^{2}$. In Section (ref), we will discuss the asymptotic variance of $\widehat{\tau}$ when $p \asymp n$. Finally, a couple of remarks are in order before closing this section. \begin{remark} In fact, the bias of $\widehat{\tau}$ is always $o (n^{-1 / 2})$, regardless of how $p$ scales with $n$: when $\widehat{\Sigma}$ is not invertible, plugging-in the pseudo-inverse of $\widehat{\Sigma}_{[k]}$ will not affect the order of the bias. This is a consequence of knowing the probability of treatment assignment by design. For the asymptotic variance of $\widehat{\tau}$ to be no greater than that of $\widehat{\tau}_{\mathrm{unadj}}$, a sufficient condition is to ensure that $\widehat{\tau}$ has the same asymptotic variance as $\sigma^{2}$, which still needs $p = o (n)$. \end{remark} \begin{remark} We use the original $X_i$ rather than the centered $\check X_i = X_i - \bar{X}_{[k]}$ to construct the second-order $U$-statistics. In the superpopulation model under random designs, centered covariates $\check X_i$ inject more complex dependencies that compound the analysis of $U$-statistics. However, we expect that the statistical properties of the resulting $\widehat{\tau}$ continue to hold. A similar observation should hold when using the centered outcomes $\check Y_i = \mathbbm{1}\{A_i=a\}Y_i - \bar{Y}_{[k]a}$. Centering $Y$ is expected to change only the form of the asymptotic variance $\sigma^{2}$ by replacing $Y(a)$ with $Y(a) - E_{[k]}[Y(a)]$, for $a=0,1$. \end{remark} \section{Valid Inference with the Bias-Corrected Estimator} In this section, we demonstrate how to estimate the variance of our new bias-corrected adjusted estimator $\widehat\tau$, so that it can be used in subsequent tasks related to hypothesis testing or statistical inference. Recall from Section (ref) and Theorem (ref) that the asymptotic variance $\sigma^2$ can be decomposed into the following terms: $\sigma^2 = \zeta^2_{H} + \zeta^2_{\mathrm{I},r}(\pi_{[k]}) + \zeta^2_{\mathrm{II}}$. The general strategy that we follow is to estimate each component in $\sigma^{2}$ by a sample average or a $U$-statistic. We now explain how to estimate each of the components of $\widehat{\sigma}^{2}$. First, we estimate $\zeta^2_{H}$ by \[\widehat\zeta^2_{H} = \sum\limits_{k=1}^K p_{n[k]}\Big\{\big(\bar Y_{[k]1} - \sum\limits_{k'=1}^K p_{n[k']}\bar Y_{[k']1}\big) - \big(\bar Y_{[k]0} - \sum\limits_{k'=1}^K p_{n[k']}\bar Y_{[k']0}\big)\Big\}^2,\] because $\zeta_{H}^{2} = E \{E_{[k]} \{r_{i} (1) - r_{i} (0)\} - E \{r_{i} (1) - r_{i} (0)\}\}^{2}$ as stated in Proposition (ref), and we simply estimate each population mean by a corresponding sample average. For the second term $\zeta^2_{\mathrm{I},r}(\pi_{[k]})$, we first recall that \[\zeta^2_{\mathrm{I},r}(\pi_{[k]}) = \sum\limits_{k=1}^K p_{[k]}\Big(\frac{\sigma^2_{[k]r(1)}}{\pi_{[k]}} + \frac{\sigma^2_{[k]r(0)}}{1-\pi_{[k]}}\Big),\] where $r_i(a) = Y_i(a) - X_i^{\top}\beta_{[k]}$, $\beta_{[k]} = (1-\pi_{n[k]})\beta_{[k]}(1) + \pi_{n[k]}\beta_{[k]}(0)$, $\beta_{[k]}(a) = \Sigma_{[k]}^{-1} \eta_{[k]}(a)$, and $\eta_{[k]} (a) = E_{[k]}\{XY(a)\}$. In Appendix (ref), we show that $\zeta^2_{\mathrm{I}, r} (\pi_{[k]})$ can be equivalently represented as \[\zeta^2_{\mathrm{I},r}(\pi_{[k]}) = \sigma^2_Y(\pi_{[k]}) - \sigma^2_{\mathrm{I},\eta}(\pi_{[k]}) - 2 \sigma_{\mathrm{I},\eta(1)\eta(0)},\] where \begin{align*} & \sigma^2_{Y}(\pi_{[k]}) = \sum\limits_{k=1}^K p_{[k]}\Big\{\frac{1}{\pi_{[k]}} \sigma^2_{[k]Y(1)} + \frac{1}{1-\pi_{[k]}} \sigma^2_{[k]Y(0)}\Big\}, \\ & \sigma^2_{\mathrm{I},\eta}(\pi_{[k]}) = \sum\limits_{k=1}^K p_{[k]}\Big\{\frac{1-\pi_{[k]}}{\pi_{[k]}} \eta^{\top}_{[k]}(1)\Sigma_{[k]}^{-1}\eta_{[k]}(1) + \frac{\pi_{[k]}}{1-\pi_{[k]}} \eta^{\top}_{[k]}(0)\Sigma_{[k]}^{-1}\eta_{[k]}(0)\Big\}, \\ & \sigma_{\mathrm{I},\eta(1)\eta(0)} = \sum\limits_{k=1}^K p_{[k]}\Big\{\eta^{\top}_{[k]}(1)\Sigma_{[k]}^{-1}\eta_{[k]}(0)\Big\}. \end{align*} As $\zeta_{H}^{2}$, $\sigma^2_Y(\pi_{[k]})$ can be estimated by the corresponding sample variance. If we directly plug in the OLS coefficients $\widehat\beta_{[k]}(a),\ a\in\{0,1\}$, then the estimators of $\sigma^2_{\mathrm{I},\eta}(\pi_{[k]})$ and $\sigma_{\mathrm{I},\eta(1)\eta(0)}$ are $V$-statistics, which leads to non-negligible bias for the same reason as for $\widehat{\tau}_{\rm OLS}$. Therefore, we can estimate $\zeta^2_{\mathrm{I},r}(\pi_{[k]})$ using $U$-statistics instead, following the same type of construction as our point estimate $\widehat{\tau}$. We denote the resulting estimator of $\zeta^2_{\mathrm{I},r}(\pi_{[k]})$ by $\widehat\zeta^2_{\mathrm{I},r}(\pi_{[k]})$, which is \begin{align*} \widehat\zeta^2_{\mathrm{I},r}(\pi_{n[k]}) =\widehat\sigma^2_Y(\pi_{n[k]}) - \widehat\sigma^2_{\mathrm{I},\eta}(\pi_{n[k]}) - 2 \widehat\sigma_{\mathrm{I},\eta(1)\eta(0)}. \end{align*} The explicit forms of $\widehat\sigma^2_Y(\pi_{n[k]})$, $\widehat\sigma^2_{\mathrm{I},\eta}(\pi_{n[k]})$, and $ \widehat\sigma_{\mathrm{I},\eta(1)\eta(0)}$ are given in Appendix (ref). Finally, in terms of $\zeta^2_{\mathrm{II}} = \zeta^2_{\mathrm{II}, Y (1)} + \zeta^2_{\mathrm{II}, Y (0)} - 2 \zeta^2_{\mathrm{II}, Y [1, 2]} (1, 0, 1)$, we apply the same estimation strategy as for $\zeta^2_{\mathrm{I},r}(\pi_{[k]})$. Specifically, we estimate $\sigma^2_{\mathrm{II},Y(a)[1]}$ by the sample mean, while $\sigma^2_{\mathrm{II},Y(a)[1,2]}$ and $\sigma^2_{\mathrm{II},Y(1,0)}$ are estimated by $U$-statistics. Their explicit forms are also deferred to Appendix (ref). The estimator of $\zeta^2_{\mathrm{II}}$ is simply denoted as $\widehat\zeta^2_{\mathrm{II}}$. In summary, our proposed variance estimator $\widehat{\sigma}^{2}$ of $\sigma^{2}$ reads as follows: $\widehat\sigma^2 = \widehat\zeta^2_{H} + \widehat\zeta^2_{\mathrm{I},r}(\pi_{n[k]}) + \widehat\zeta^2_{\mathrm{II}}$. Our last main theoretical result shows that $\widehat{\sigma}^{2}$ is a consistent estimator of $\sigma^{2}$ under mild assumptions, so we can use $\widehat{\tau}$ and $\widehat{\sigma}^{2}$ to conduct downstream inference. \begin{theorem} Under Assumptions (ref)--(ref), and assuming that $p = o(n)$, the following hold: \begin{enumerate}[label = (\alph*)] • $\widehat\sigma^2$ is a consistent estimator of $\sigma^2$: $\widehat\sigma^2\stackrel{P}\to \sigma^2$. • As a consequence, after scaled by $\widehat{\sigma}$ instead of $\sigma$, $\widehat{\tau}$ is still $\sqrt{n}$-CAN: \begin{align*} \sqrt{n} (\widehat\tau - \tau)/\widehat\sigma \stackrel{d}\to \operatorname{N}(0,1). \end{align*} • Finally, the nominal $(1 - \alpha) \times 100\%$ large-sample Wald confidence interval $$\widehat{\rm CI}_{\alpha} = \left[ \widehat{\tau} - z_{1 - \alpha / 2} \dfrac{\widehat{\sigma}}{\sqrt{n}}, \widehat{\tau} + z_{1 - \alpha / 2} \dfrac{\widehat{\sigma}}{\sqrt{n}} \right]$$ retains the correct coverage probability as $n \rightarrow \infty$. \end{enumerate} \end{theorem} \begin{remark} Similar to Condition (ref), $\widehat\sigma^2$ is still a consistent estimator of $\sigma^2$ if we replace $\widehat{\Sigma}_{[k]}$ by any estimator $\widetilde{\Sigma}_{[k]}$ of $\Sigma_{[k]}$ if: (a) $\widetilde{\Sigma}_{[k]}$ is invertible almost surely and (b) \begin{equation} E_{[k]} \{X_i^{\top} (\widetilde{\Sigma}_{[k]}^{-1} - \Sigma_{[k]}^{-1}) X_i\} = o(p),\ E_{[k]}\{X_i^{\top} (\widetilde{\Sigma}_{[k]}^{-1} - \Sigma_{[k]}^{-1})X_i\}^2 = O\left(\frac{p^3}{n}\right). \end{equation} We verify that $\widehat{\Sigma}_{[k]}$ satisfies Condition (ref) when $p=o(n)$ in Lemma (ref) (see Appendix (ref)). \end{remark} \begin{remark} When $p \asymp n$, $\widehat{\sigma}^2$ is not a consistent estimator of $\sigma^2$ and their difference is of order $p / n$. Based on the empirical results in Section (ref), $\widehat{\sigma}^{2}$ may slightly overestimate $\sigma^{2}$, which makes sure that the resulting Wald CI is conservative. But we can still observe that $\widehat{\sigma}^{2}$ is generally less than the asymptotic variance of $\widehat{\tau}_{\mathrm{unadj}}$, achieving an efficiency gain. \end{remark} \section{Numerical Experiments} \subsection{Simulation studies} In this section, we evaluate the performance of our proposed estimators, including both the oracle estimator $\widehat\tau_{\mathrm{ora}}$ and the feasible estimator $\widehat{\tau}$, and compare them with the unadjusted estimator $\widehat\tau_{\mathrm{unadj}}$ and the OLS estimator $\widehat\tau_{\rm OLS}$. We use $\widehat\sigma^2$ proposed in Theorem (ref) to estimate the variance of $\widehat\tau$. For the oracle estimator $\widehat\tau_{\mathrm{ora}}$, the corresponding variance estimator $\widehat\sigma_{\mathrm{ora}}^{2}$ is the same as $\widehat{\sigma}^{2}$ except that we use $\Sigma_{[k]}$ instead of $\widehat{\Sigma}_{[k]}$. For $\widehat\tau_{\mathrm{unadj}}$ and $\widehat\tau_{\rm OLS}$, we use the variance estimators $\widehat\sigma^2_{\rm unadj}$ and $\widehat\sigma^2_{\rm OLS}$ proposed in Bugni2018 and gu2023regression, respectively. The simulation is carried out under stratified block randomization, using a categorical covariate $X_{1} \in \{1, 2, 3, 4\}$ as the stratification variable, with corresponding probabilities $\{0.2, 0.2, 0.3, 0.3\}$. Specifically, we consider the following data-generating processes. \paragraph{Model 1.} In Model 1, the outcome is continuous, where for $i = 1, \cdots, n$, \begin{align*} Y_i(0) & = X_{1i}+ 2X_{0i}^{\top} \beta_0 - 0.5X_{0i,4}^2 + \epsilon_{0,i}, \text{ and} \\ Y_i(1)&= X_{1i}+0.05 X_{0i}^{\top} \Sigma_0 ^{-1} X_{0i} + \epsilon_{1,i}. \end{align*} Here, $X_{0i}$ is a random vector of dimension $p_{0}$ and $X_{0i}\sim t_5(0,\Sigma)$, with $p_0=30$ and $\Sigma_{i,j} = 0.1^{|i-j|}$. Each component of $\beta_0$ is generated from the uniform distribution on $[-1,1]$ and then rescaled to $\lVert\beta_0\rVert_2=1$. Finally, $\epsilon_{0,i}, \epsilon_{1,i} \sim \operatorname{N}(0,0.1^2)$. \paragraph{Model 2.} In Model 2, the outcome is binary, where for $i = 1, \cdots, n$, \begin{align*} \Pr(Y_i(0)=1) & = \mathrm{expit}\Big(-1+ X_0^{\top}\beta_0 -2X_{0i,1}^2 \Big), \text{ and} \\ \Pr(Y_i(1)=1)&= \mathrm{expit}\Big(-3 +X_0^{\top} \beta_1+2X_{0i,2}^2 + 0.5X_{0i,3}^4 \Big). \end{align*} Here $\mathrm{expit}(z) = 1/(1+\exp(-z))$. $X_0$ is generated in the same fashion as in Model 1. We take $\beta_{1,j} = 0.5,\ j=1\dots p_0$, $\beta_{0,j} = 1.5,\ j=1,2,\dots,p$. We set the sample size to $n = 1000$. The allocation ratio between the treatment and control groups is $1:1$, leading to an expected control (treatment) group sample size of $n_{[k_{\min}]0}=100$ ($n_{[k_{\min}]1} = 100$) in the smallest stratum. We denote the covariates used in the analysis stage as $X_{2i}$ for $i = 1, \cdots, n$. When $p \leq p_0$, $X_{2i}$ consists of the first $p$ components of $X_{0i}$; when $p > p_0$, $X_{2i} = [X_{0i}, Z_i]$, where $Z_i$ is a random vector of dimension $(p - p_0)$ with components $Z_{i, q} \stackrel{i.i.d.}{\sim} t_5,\ q=1,2,\dots,p-p_0$. We set $p = \lceil rn\rceil$, where $r \in \{0.02,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7\}$. We let $r_0 = p_0/n=0.3$ be the ratio between the dimension of $X_{0}$ and the sample size $n$. In the simulation, we draw $R = 2000$ Monte Carlo replicates. The performance of different estimators is evaluated in terms of absolute bias, Monte Carlo standard deviation (SD), the ratio of the Monte Carlo SD to the estimated standard error (SE) ${\rm sd} / {\rm se}$, and the coverage probability (CP) of the associated nominal 95% Wald confidence interval (CI). We also report the Monte Carlo CP, which is computed using the Monte Carlo SD instead of the estimated SE to evaluate the impact of bias on inference. \begin{itemize} • In Figures (ref)(a) and (ref)(a), we report the absolute bias of all estimators included in the comparison. The bias of $\widehat\tau_{\rm OLS}$ increases as the ratio $r = p / n$ increases, until it reaches $r_0 = 0.3$. This is expected because the additional covariates $Z$ in $X_{2}$ are independent of $Y (a)$ and thus will not further increase the bias, as noted in Remark (ref). All other estimators, $\widehat{\tau}_{\mathrm{unadj}}$, $\widehat{\tau}_{\mathrm{ora}}$, and $\widehat{\tau}$, have a bias close to zero when varying $r$. As shown in Figures (ref)(e) and (ref)(e), when evaluating the impact of bias on CP, only the Wald CI centered at $\widehat\tau_{\rm OLS}$ undercovers at around 90% when $r \ge 0.3$. • In Figure (ref)(b), the Monte Carlo SD decreases with increasing $r$ for all three adjusted estimators $\widehat\tau_{\rm OLS}$, $\widehat\tau_{\mathrm{ora}}$, and $\widehat\tau$. All three estimators are more efficient than $\widehat\tau_{\mathrm{unadj}}$. The efficiency gain of $\widehat\tau_{\rm OLS}$ relative to $\widehat\tau$ also increases with $r$, consistent with the claim in Proposition (ref) that the additional term $\zeta^2_{\mathrm{II}}$ is of order $O(p/n)$ and $\zeta^2_{\mathrm{II}} \geq 0$. However, in Figure (ref)(b), $\widehat\tau_{\rm OLS}$ has a larger SD compared to $\widehat\tau_{\rm unadj}$ when $r$ reaches 0.6, which could cause harm. In contrast, $\widehat\tau$ retains efficiency gains compared to $\widehat\tau_{\rm unadj}$. • In Figures (ref)(c) and (ref)(d), we examine the performance of our variance estimators and their corresponding CP. The ratios $\rm sd/se$ for all estimators are within the range $[0.85,1.05]$. As $r$ increases, $\widehat\sigma$ shows a larger decrease in $\rm sd/se$ compared to $\widehat\sigma_{\mathrm{ora}}$, which may be caused by the violation of Condition (ref) by $\widehat{\Sigma}_{[k]}$ when $p = o (n)$ does not hold. In other words, $\widehat{\sigma}$ can overestimate the actual SD when $p$ is near $n$, which is not as concerning as long as $\widehat{\sigma}$ is still smaller than the SD of $\widehat{\tau}_{\mathrm{unadj}}$. The Wald CIs associated with both $\widehat\tau_{\mathrm{ora}}$ and $\widehat\tau_{\mathrm{unadj}}$ have CP around 95%, while $\widehat\tau$ is slightly conservative because its SE overestimates SD. For $\widehat\tau_{\rm OLS}$, the larger bias and the underestimated SE together lead to the undercoverage of its Wald CI for almost all $r$. • In Figures (ref)(a)--(e), we examine the performance of our estimators when the outcome is binary. The results similarly show that our proposed estimator $\widehat\tau$ still has negligible bias, delivers valid inference, and is more efficient than $\widehat\tau_{\mathrm{unadj}}$, further corroborating the practical utility of our proposed estimator. \end{itemize} \begin{figure}[ht] \caption{Simulation results for different $p/n$ ratios of Model 1. Panels (a)–(e) summarize the finite-sample performance of the four estimators: $\widehat\tau_{\rm unadj}$, $\widehat\tau_{\rm OLS}$, $\widehat\tau_{\mathrm{ora}}$, and $\widehat\tau$. (a) Absolute bias; (b) Sampling standard deviation; (c) Ratio of empirical standard deviation to the estimated standard error; (d) Coverage probability; and (e) Monte Carlo coverage probability.} \end{figure} \begin{figure}[ht] \caption{Simulation results for different $p/n$ ratios of Model 2. Panels (a)–(e) summarize the finite-sample performance of the four estimators: $\widehat\tau_{\rm unadj}$, $\widehat\tau_{\rm OLS}$, $\widehat\tau_{\mathrm{ora}}$, and $\widehat\tau$. (a) Absolute bias; (b) Sampling standard deviation; (c) Ratio of empirical standard deviation to the estimated standard error; (d) Coverage probability; and (e) Monte Carlo coverage probability.} \end{figure} \subsection{A semi-synthetic RCT data analysis} In this section, we analyze a semi-synthetic data simulated from an RCT of nefazodone and the cognitive behavioral-analysis system of psychotherapy (CBASP) keller2000comparison. The purpose of this trial is to compare the effects of nefazodone, CBASP, and their combination on chronic depression. We take the combination as the treatment group indexed as $A = 1$ and nefazodone as the control group indexed as $A = 0$. The outcome of interest $Y$ is FinalHAMD, the final score of the 24-item Hamilton rating scale for depression. To generate the semi-synthetic data, we fit the real data with an additive model using the function \texttt{gam} from the \texttt{R} package \texttt{mgcv}. We stratify the data using GENDER and select seven covariates to fit the model, including AGE, HAMA, HAMA_SOMATI, HAMD17, HAMD24, Mstatus2, and TreatPD. We include linear and quadratic transformations of continuous covariates and their interactions, resulting in $p = 30$ covariates in total. We model quadratic and interaction terms with a cubic spline, while the remaining five continuous covariates enter the model linearly. We then sample 600 units with replacement from the real data as our semi-synthetic data. Next, we implement stratified block randomization, using GENDER as the stratified variable with an allocation ratio of 1:1 and a block size of 6. We compare the performance of three estimators based on 2000 Monte Carlos: $\widehat\tau_{\rm unadj}$, $\widehat\tau_{\rm OLS}$ and our new estimator $\widehat\tau$. As shown in Table (ref), the biases of $\widehat\tau_{\mathrm{unadj}}$ and $\widehat\tau$ are close to zero, while $\widehat\tau_{\rm OLS}$ exhibits non-negligible bias. Meanwhile, $\widehat\tau$ is about 30% more efficient than $\widehat\tau_{\mathrm{unadj}}$. The standard error (SE) of $\widehat{\tau}$ is close to its Monte Carlo SD, and is still smaller than that of $\widehat{\tau}_{\mathrm{unadj}}$. However, $\widehat\tau_{\rm OLS}$ suffers from both a larger bias and an underestimated SE, which together lead to an undercoverage (90%) of the nominal 95% Wald CI centered at $\widehat{\tau}_{\rm OLS}$. Taken together, the semi-synthetic data analysis further consolidates our main conclusions. \begin{table}[H] \caption{Results of the semi-synthetic RCT data (all scaled by a factor of $10^2$)} \begin{tabular}{lllll} \Xhline{3\arrayrulewidth} & Bias & SD & SE & CP \\ \Xhline{3\arrayrulewidth} $\widehat\tau_{\rm unadj}$ & 0.4 & 25.2 & 24.7 & 94.8 \\ $\widehat\tau_{\rm OLS}$ & -3.2 & 14.9 & 12.6 & 89.3 \\ $\widehat\tau$ & 0.3 & 19.8 & 20.7 & 96.0\\ \Xhline{3\arrayrulewidth} \end{tabular} \end{table} \section{Discussion} In this paper, we develop a new covariate-adjusted treatment effect estimator $\widehat{\tau}$ based on second-order $U$-statistics for CAR under the superpopulation model, together with its consistent variance estimator. We show that $\widehat{\tau}$ is $\sqrt{n}$-CAN and is more efficient than $\widehat{\tau}_{\mathrm{unadj}}$ in the assumption-lean setting if $p = o (n)$. Empirical results indicate that $\widehat{\tau}$ is competitive in terms of various metrics compared to popular benchmarks. \paragraph{Practical recommendation} In view of the results obtained in our paper, we recommend the following practice, which addresses the quote on page (ref) to some extent. When the number $p$ of adjusted covariates is relatively large compared to $n$, it is recommended to use our proposed point and interval estimates $\widehat{\tau}$ and $\widehat{\sigma}^{2}$ to evaluate the ATE for CAR. The bias of $\widehat{\tau}$ is always negligible compared to its standard deviation and the asymptotic variance of $\widehat{\tau}$ never exceeds that of $\widehat{\tau}_{\mathrm{unadj}}$. This suggestion is especially justified if the covariates being adjusted are prognostic factors of the outcome, as the efficiency gain of $\widehat{\tau}$ is guaranteed and the standard OLS estimator $\widehat{\tau}_{\rm OLS}$ in general has non-negligible bias. When $p$ is small compared to $n$, $\widehat{\tau}$ and $\widehat{\sigma}^{2}$ are still safe to use, having finite sample performance comparable to $\widehat{\tau}_{\rm OLS}$ and $\widehat{\sigma}_{\rm OLS}^{2}$. In summary, we recommend using the proposed estimator $\widehat{\tau}$ as the primary analysis, due to its efficiency gains over $\widehat{\tau}_{\mathrm{unadj}}$ and its reduced bias relative to $\widehat{\tau}_{\rm OLS}$, particularly when the number of covariates $p$ is moderate to large compared to the sample size $n$. At the same time, we suggest reporting $\widehat{\tau}_{\mathrm{unadj}}$ and $\widehat{\tau}_{\rm OLS}$, together with their corresponding CIs, as part of a sensitivity analysis, to provide a more complete picture of their bias–efficiency performance in the given analysis. Nevertheless, any such analysis strategy should be discussed with the relevant regulatory agencies when the results are intended for regulatory submission, in accordance with current regulatory requirements. \paragraph{Extensions} Our theoretical results are stated under the assumption $p = o (n)$, which paves the way to the analysis of the more difficult regime $p \asymp n$. When $p < n / K$, $\widehat{\Sigma}_{[k]}$ is invertible with high probability and $\widehat{\tau}$ remains nearly unbiased. But, as we mentioned in Remark (ref), the asymptotic variance of $\widehat{\tau}$ may exceed that of $\widehat{\tau}_{\mathrm{unadj}}$ or $\widehat{\tau}_{\rm OLS}$. In the superpopulation model, random matrix theory is needed for a more delicate analysis in this regime. When $p > n$, the generalized inverse $\widehat{\Sigma}_{[k]}^{\dag}$ of the sample Gram matrix estimator $\widehat{\Sigma}_{[k]}$ or various shrinkage estimators of $\Sigma_{[k]}^{-1}$ can be used liu2025augmented, ledoit2004well, ding2024eigenvector, and the resulting estimator based on $U$-statistics is still nearly unbiased. The same challenge lies in the analysis of the asymptotic variance, the characterization of asymptotic normality, and the construction of valid CIs zheng2025perturbed, which we leave to future work. Let $b_{[k]} (x; a) = E_{[k]} (Y_{i} (a) | X = x)$ denote the true outcome regression for the treatment group $a$ in the stratum $k$. When $p$ is fixed, a modified version of $\widehat{\tau}$, denoted by $\widehat{\tau}_{\rm eff}$, achieves the semiparametric efficiency bound of $\tau$ under CAR. $\widehat{\tau}_{\rm eff}$ simply replaces $X$ by a set of basis transformations of $X$, say $\{z_{l}\}_{l = 1}^{k}$ with $k = k (n) \rightarrow \infty$ as $n \rightarrow \infty$. Under standard smoothness assumptions on $b_{[k]}$, with $k$ growing with $n$, the asymptotic variance of $\widehat{\tau}_{\rm eff}$ simply replaces $r_{i} (a) = Y_{i} (a) - X_{i}^{\top} \beta_{[k]}$ by $Y_{i} (a) - b_{[k]} (X_{i}; a)$ in $\sigma^{2}_{\rm OLS}$, defined in Proposition (ref), which is the semiparametric efficiency bound of $\tau$ under CAR rafi2023efficient. \paragraph{Other future directions} We conclude our paper by discussing several other future directions. An immediate question is to consider nonlinear outcome working models, probably using generalized linear models with canonical link functions as a starting point guo2023generalized, cohen2024no. Finally, another important problem to investigate urgently is the integration of variable selection methods van2024automated into our procedure without deteriorating the statistical properties established here. \putbib[reference]