Covariate Adjustment in Randomized Experiments Motivated by Higher-Order Influence Functions

abstract{ Higher-Order Influence Functions (HOIF), developed in a series of papers over the past twenty years, is a fundamental theoretical device for constructing rate-optimal causal-effect estimators from observational studies. However, the value of HOIF for analyzing well-conducted randomized controlled trials (RCT) has not been explicitly explored. In the recent U.S. Food and Drug Administration and European Medicines Agency guidelines on the practice of covariate adjustment in analyzing RCT, in addition to the simple, unadjusted difference-in-mean estimator, it was also recommended to report the estimator adjusting for baseline covariates via a simple parametric working model, such as a linear model. However, when the number of baseline covariates $p$ is large, the recommendation is somewhat murky. In this paper, we show that HOIF-motivated estimators for the treatment-specific mean have significantly improved statistical properties compared to popular adjusted estimators in practice when $p$ is relatively large relative to the sample size $n$. We also characterize the conditions under which the HOIF-motivated estimator improves upon the unadjusted one. More importantly, we demonstrate that several state-of-the-art adjusted estimators proposed recently can be interpreted as particular HOIF-motivated estimators, thereby placing these estimators in a more unified framework. Numerical and empirical studies are conducted to corroborate our theoretical findings. An accompanying R package can be found on \href{https://cran.r-project.org/web/packages/HOIFCar/index.html}{CRAN}.}

bibunit[plainnat] \section{Introduction} Evidence from randomized clinical trials (RCT) is widely regarded as the gold standard for evaluating treatment effects in comparative effectiveness research. Complete randomization, along with its extensions, such as covariate-adaptive randomization and rerandomization pocock1975sequential, morgan2012rerandomization, ma2024new, relieves analysts' burden of justifying the unconfoundedness assumption. Furthermore, the true propensity score in RCT is known to data analysts aronow2025nonparametric, offering another important advantage over observational studies. The most straightforward RCT designs include Completely Randomized Experiments (CRE) and Bernoulli sampling, in which treatments are randomly assigned without leveraging any covariate information. Under CRE or Bernoulli sampling, the Difference-in-Mean estimator, referred to as the unadjusted estimator henceforth, is unbiased, $\sqrt{n}$-consistent, centered and asymptotically normal (CAN) for the average treatment effect (ATE), where $n$ denotes the sample size. However, modern trials typically collect multiple baseline measurements. It is now widely recognized that, both theoretically and empirically, adjusting for baseline covariates (especially those prognostic factors affecting potential outcomes) in either the design stage or the analysis stage can yield greater efficiency than the unadjusted estimator (see e.g. yang2001efficiency, zhang2008improving, lin2013agnostic, ma2020statistical, zhao2021covariate, ma2022regression, ye2023toward and references therein). In this paper, we focus only on the adjusting for baseline covariates in the analysis stage and primarily on CRE. But we will mention parallel results under Bernoulli sampling when doing so facilitates understanding. Covariate adjustment methods in RCT have a long history in statistics freedman2008regressiona, freedman2008regressionb. Adjustment using linear or simple parametric working models ma2022regression, ye2023toward has been endorsed in the latest statistical analysis guidelines issued by the U.S. Food and Drug Administration (FDA) FDA2023. The rationale for the improved efficiency through covariate adjustment, even with a likely misspecified linear working model, can be clearly articulated within the superpopulation framework. As noted in richardson2014causal, the semiparametric theory developed in robins1994estimation motivated key works on covariate adjustment in RCT tsiatis2008covariate, moore2009covariate. Specifically, the theory indicates that the Augmented Inverse Probability Weighting (AIPW) estimator, rather than the Inverse Probability Weighting (IPW) estimator, achieves the smallest possible asymptotic variance, more precisely the Semiparametric Variance Bound (SVB), when the outcome model is correctly specified. The SVB of the ATE is characterized by its (first-order) influence function, a central concept in semiparametric theory. The variance reduction property of the AIPW estimator, constructed based on the influence function, can be explained geometrically: the augmented term subtracted off the IPW estimator can be viewed as a projection of the IPW estimator onto a specific subspace, thereby reducing variance through the contracting norm property of the projection. This intuition holds even with a misspecified linear model. In CRE, the unadjusted estimator coincides with the IPW estimator, while the adjusted estimator (see (ref) for its specific form) is algebraically equivalent to the AIPW estimator, using an estimated linear working model for the outcome regression. Thus, the previous geometric reasoning shows variance reduction for the adjusted over unadjusted estimator in fixed $p$ (dimension of baseline covariates), large $n$ regimes. Section (ref) elaborates this intuition. With technical advancements, modern RCT routinely collect multidimensional baseline covariates. However, per the recent FDA guideline, when $p$ is large relative to $n$, the best practice to adjust for covariates in the analysis stage becomes murky. The recent literature has therefore gradually turned to the development of methods adjusting for higher-dimensional covariates. Since we primarily focus on the case using linear working models, we highlight some of the most relevant papers. ma2022regression and ye2023toward demonstrated that the estimator that adjusts for baseline covariates $\mathbf{x}$ and the interaction between $\mathbf{x}$ and treatment $t$ in a linear working model using OLS is CAN and more efficient or as efficient as the unadjusted estimator in various designs when $p$ is fixed and $n \rightarrow \infty$, in the superpopulation framework. jiang2025adjustments showed that the same results hold for this OLS-based estimator when $p = o (\sqrt{n})$; but when $\sqrt{n} \lesssim p \lesssim n$, a true linear outcome model needs to be further assumed, an assumption often considered unduly strong in RCTs. lei2021regression also considered covariate adjustment in a linear working model, under CRE and the design-based (or randomization-based) framework, which allows $p = O (n^{2 / 3})$ up to log-factors. Unlike jiang2025adjustments, they do not assume that the linear working model is correctly specified when $p \gtrsim \sqrt{n}$. More recently, lu2025debiased devised a debiased estimator that improves on that of lei2021regression: their estimator is never less efficient than, and can sometimes be more efficient than the unadjusted estimator when $p=o(n)$. lu2025debiased note that chang2024exact constructed a similar but exactly unbiased estimator earlier, but only with theoretical results for fixed $p$. Over the past two decades, higher-order influence functions (HOIF) have been developed as a generalization of classical semiparametric theory to construct rate-optimal estimators for parameters like the ATE from observational data robins2008higher. Given the significant role of classical semiparametric theory in covariate adjustment, a natural inquiry arises: can the HOIF of the ATE also inform the development of covariate adjustment methods in RCT? In this paper, we provide an affirmative answer to this question. For ease of presentation, we focus solely on the treatment-specific mean in the treatment arm (and, by symmetry, the control arm) instead of the ATE. In our accompanying R package available on \href{https://cran.r-project.org/web/packages/HOIFCar/index.html}{CRAN}, the corresponding point and interval estimators of ATE are also available. \subsection*{Main contributions and organization} The main contributions of our paper are summarized below: \begin{enumerate}[leftmargin=0.5cm,topsep=0.25pt] • We demonstrate that a HOIF-motivated adjusted estimator has improved statistical properties over standard adjusted or unadjusted estimators. The theoretical analyses are “conceptually simple”, involving only elementary calculations. We also develop a variance estimator of this estimator, deferred to Appendix (ref) due to space limitation. This new variance estimator is recommended because of the improved coverage when $n$ is small, based on empirical observations from our simulation studies. • Our impression is that the HOIF theory remains elusive even among statisticians, so, in Section (ref), we review it in an accessible manner to garner more interest from statisticians. The more important point we hope to make is not to propose a new estimator. Instead, we show that several state-of-the-art adjusted estimators mentioned above are specific HOIF-motivated estimators, placing them in a more unified framework; see Table (ref). • We develop an accompanying user-friendly R package that is available on \href{https://cran.r-project.org/web/packages/HOIFCar/index.html}{CRAN}, which delivers both point and interval estimators for ATE and treatment/control-specific means. \end{enumerate} The rest of the paper is organized as follows. Section (ref) introduces the basic setup. Section (ref) presents a variety of HOIF-motivated adjusted estimators and their statistical properties (Sections (ref) and (ref)), and also draws connections to other state-of-the-art estimators (Section (ref)). Simulation studies and real data analysis are carried out in Section (ref). We conclude the paper in Section (ref). The Appendix contains supplementary technical and empirical results. \section{Notation and Basic Setup} \subsection*{Notation} In this paper, we denote sample size as $n$ and baseline covariates dimension as $p$. Design-based quantities have superscript “$\mathsf{d}$” to distinguish from superpopulation frameworks; e.g., ${\mathbb{E}}^{\mathsf{d}}$, $\mathrm{bias}^{\mathsf{d}}$, and $\mathrm{var}^{\mathsf{d}}$ for expectation, bias, variance. Superscripts are omitted when unambiguous. We also adopt the common asymptotic and stochastic asymptotic notation, including $O (\cdot)$, $o (\cdot)$, $O_{{\mathbb{P}}} (\cdot)$, $o_{{\mathbb{P}}} (\cdot)$, with ${\mathbb{P}}$ the true distribution. For square matrix $\mathbf{M}$, $\mathsf{tr}(\mathbf{M})$ is the trace, $\Vert \mathbf{M} \Vert_{\mathrm{op}}$ the operator norm, and $\mathbf{M}^{-}$ the inverse or pseudoinverse. Vector norm $\Vert \mathbf{v} \Vert$ is the $\ell_{2}$-norm, and $i \in [n]$ denotes $i=1,\dots,n$. For the $n \times p$ covariate matrix $\mathbf{X} \coloneqq (\mathbf{x}_{1}, \cdots, \mathbf{x}_{n})^{\top}$ with $\mathbf{x}_{i} \in {\mathbb{R}}^{p}$ for $i\in [n]$, let $\bar{\mathbf{x}} \coloneqq n^{-1} \sum_{i = 1}^{n} \mathbf{x}_{i} \in {\mathbb{R}}^{p}$ be the row-wise average vector. Following standard practice in covariate adjustment in randomized experiments, we center the covariate/design matrix by $\bar{\mathbf{x}}$ to obtain $\mathbf{X}_{c} \coloneqq \left( \mathbf{x}_{1} - \bar{\mathbf{x}}, \cdots, \mathbf{x}_{n} - \bar{\mathbf{x}} \right)^{\top}$. We then define the $n \times n$ “hat” projection matrix $\mathbf{H} \coloneqq \mathbf{X}_{c} \widehat{\bm{\Sigma}}^{-} \mathbf{X}_{c}^{\top} \equiv \left( H_{i, j}, 1 \leq i, j \leq n \right), \text{ where } \widehat{\bm{\Sigma}} \coloneqq \mathbf{X}_{c}^{\top} \mathbf{X}_{c}.$ In the random design setting of the superpopulation framework, denote $\bm{\mu} \coloneqq {\mathbb{E}} \mathbf{x}$ and $\bm{\Sigma} \coloneqq n {\mathbb{E}} (\mathbf{x} - \bm{\mu}) (\mathbf{x} - \bm{\mu})^{\top}$. For any vector $\mathbf{v} = (v_{1}, \cdots, v_{n})^{\top}$ of length $n$, define \begin{align*} V_{n} (v) \coloneqq \frac{1}{n} \sum_{i = 1}^{n} v_{i}^{2} - \frac{1}{n (n - 1)} \sum_{1 \leq i \neq j \leq n} v_{i} v_{j} \equiv \frac{1}{n - 1} \sum_{i = 1}^{n} \left( v_{i} - \bar{v} \right)^{2}, \end{align*} as the sample variance of $\mathbf{v}$, where $\bar{v} \coloneqq n^{-1} \sum_{i = 1}^{n} v_{i}$. \subsection*{Basic Setup} Throughout this paper, we observe the data matrix: $\mathbf{O} \in {\mathbb{R}}^{n \times (p + 2)} \coloneqq \left( \mathbf{o}_{1}, \cdots, \mathbf{o}_{n} \right)^{\top}$, where $ \mathbf{o}_{i} \coloneqq \left( \mathbf{x}_{i}^{\top}, t_{i}, y_{i} \right)^{\top} \in {\mathbb{R}}^{p + 2}$ for $i \in [n]$, with $t$ and $y$ denoting the treatment indicator and the outcome. Let $\mathbf{t} \coloneqq (t_{1}, \cdots, t_{n})^{\top}$ and $\mathbf{y} \coloneqq (y_{1}, \cdots, y_{n})^{\top}$. We always assume that $p < n$ and $\lim_{n \rightarrow \infty} p / n = \alpha \in [0, 1)$ without further mentioning this assumption. Our result thus covers the relatively more challenging “proportional asymptotic” regime at least for $p < n$. Without loss of generality, we take $t \in \{0, 1\}$ and $y \in {\mathbb{R}}$. In RCTs, $\mathbf{t}$ is determined by exogenous randomization, hence under the investigator's control. In particular, we mainly consider the CRE, which assign $n_{1}$ out of $n$ subjects uniformly into the treatment group ($t = 1$), and the remaining $n_{0}$ subjects into the control group. Let $\pi_{1} \coloneqq n_{1} / n$ and $\pi_{0} \coloneqq 1 - \pi_{1}$ be the treatment and control proportions respectively. Denoting the potential outcome vector under treatment as $\mathbf{y} (1) = (y_{1} (1), \cdots, y_{n} (1))^{\top}$ and under control as $\mathbf{y} (0) = (y_{1} (0), \cdots, y_{n} (0))^{\top}$. By the standard consistency assumption, the observed outcome and potential outcomes are connected by $y \equiv t y (1) + (1 - t) y (0)$. In RCT, randomization licenses the use of observables $\mathbf{O}$ to identify certain causal quantities defined via potential outcomes. In this paper, we assume : \begin{assumption}[Randomization] $\mathbf{t}$ is assigned via CRE or the Bernoulli sampling, so $\mathbf{t} \protect\mathpalette{\protect\independenT}{\perp} \{\mathbf{x}, \mathbf{y} (0), \mathbf{y} (1)\}$. \end{assumption} For example, if one is interested in the treatment-specific mean $\bar{\tau} \coloneqq {\mathbb{E}}^{\mathsf{d}} y (1) \equiv \frac{1}{n} \sum_{i = 1}^{n} y_{i}(1)$, one can identify $\bar{\tau}_{1}$ via $\widehat{\tau}_{\mathsf{unadj}} \coloneqq \frac{1}{n_{1}} \sum_{i = 1}^{n} t_{i} y_{i} \equiv \frac{1}{n} \sum_{i = 1}^{n} \frac{t_{i}}{\pi_{1}} y_{i}.$ $\widehat{\tau}_{\mathsf{unadj}}$ is often called the unadjusted estimator or \textit{the IPW estimator}. It is easy to see that $\widehat{\tau}_{\mathsf{unadj}}$ is unbiased for $\bar{\tau}_{1}$. One can similarly define the design-based control specific mean and the ATE, together with their corresponding unadjusted or IPW estimators. Without loss of generality, we only consider the treatment specific mean $\bar{\tau}$ for two reasons. First, all the results hold for the control specific mean by symmetry. Second, treatment and control specific means are more primitive parameters than the ATE. Occasionally, we also consider the superpopulation framework, under which the observed data is drawn i.i.d. from a common probability distribution ${\mathbb{P}}$: $\mathbf{o}_{i} \overset{\rm i.i.d.}{\sim} {\mathbb{P}}, \,\,\,\, i \in [n].$ In the superpopulation framework, we only consider the Bernoulli sampling of the treatment assignment vector, i.e. $\{t_{i}\}_{i = 1}^{n} \overset{\rm i.i.d.}{\sim} \mathrm{Bernoulli} (\pi_{1})$. The superpopulation treatment specific mean is denoted as $\tau \coloneqq {\mathbb{E}} y (1)$. \section{HOIF-Motivated Covariate Adjustment: Statistical Intuition} Our main theoretical results (in Section (ref)) are stated under the design-based framework. However, in this section, we first explain the main intuition of using HOIF to construct adjusted estimator of $\tau$ under the superpopulation framework. In our own opinion, for most parts, the \emph{statistical intuition} gathered from the superpopulation framework can be carried over to the design-based framework. While familiar to HOIF experts, this section aims to interest practitioners in RCT analysis. Guarded by randomization, the unadjusted estimator $\widehat{\tau}_{\mathsf{unadj}}$ has already fulfilled the following desiderata: \begin{itemize}[leftmargin=0.5cm,topsep=0.25pt] • It is model-free, unbiased and has variance of order $1 / n$ under certain regularity conditions on $\mathbf{X}$ and $\mathbf{y}$ (see Assumption (ref) later); • It is CAN and the variance is easy to estimate. \end{itemize} However, $\widehat{\tau}_{\mathsf{unadj}}$ fails to leverage the information of $\mathbf{x}$. One convincing argument for using $\mathbf{x}$ comes from \textit{semiparametric theory}. It says that the variance of the following random variable, referred to as the (efficient) first-order influence function of $\tau$, characterizes the SVB of any Regular and Asymptotic Linear estimator of $\tau$: \begin{equation} \dot{\tau}_{1, \tau} \equiv \dot{\tau}_{1, \tau} (\mathbf{o}) \coloneqq \frac{t}{\pi_{1}} y - \left( \frac{t}{\pi_{1}} - 1 \right) \mathrm{b} (\mathbf{x}) - \tau, \end{equation} where $\mathrm{b} (\cdot) \coloneqq {\mathbb{E}} (y | \mathbf{x} = \cdot, t = 1)$. This motivates the AIPW estimator: \begin{equation} \widehat{\tau}_{\mathsf{aipw}} \coloneqq \frac{1}{n} \sum_{i = 1}^{n} \frac{t_{i}}{\pi_{1}} y_{i} - \left( \frac{t_{i}}{\pi_{1}} - 1 \right) \widehat{\mathrm{b}} (\mathbf{x}_{i}), \end{equation} where $\widehat{\mathrm{b}}$ estimates $\mathrm{b}$ using parametric models or machine learning algorithms bannick2025general. If $\widehat{\mathrm{b}}$ is consistent for $\mathrm{b}$ and ${\mathbb{P}}$-Donsker (intuitively speaking, sufficiently stable), $\widehat{\tau}_{\mathsf{aipw}}$ attains the SVB asymptotically. However, since $\mathrm{b}$ is not under the investigator's control, $\widehat{\tau}_{\mathsf{aipw}}$ may inflate the asymptotic variance when $\mathrm{plim}_{n \rightarrow \infty} \widehat{\mathrm{b}} \neq \mathrm{b}$. Fortunately, covariate adjustment via a possibly misspecified linear working model with the least square estimator still guarantees efficiency improvement over $\widehat{\tau}_{\mathsf{unadj}}$ when $p$ is fixed lin2013agnostic, ma2022regression, ye2023toward. An oracle version of this adjusted estimator is: \begin{equation} \widetilde{\tau}_{\mathsf{adj}} \coloneqq \frac{1}{n} \sum_{i = 1}^{n} \frac{t_{i}}{\pi_{1}} y_{i} - \left( \frac{t_{i}}{\pi_{1}} - 1 \right) (\mathbf{x}_{i} - \bm{\mu})^{\top} \bm{\beta}, \end{equation} where $\bm{\beta} \coloneqq n \bm{\Sigma}^{-} \cdot {\mathbb{E}} ((\mathbf{x} - \bm{\mu}) \frac{t}{\pi_{1}} (y - \tau))$ is the population projection of $y (1) - \tau$ onto the linear span of $\mathbf{x} - \bm{\mu}$. $\widetilde{\tau}_{\mathsf{adj}}$ has the same form as the AIPW estimator with linear outcome regression. The following result is immediate and well known robins1994estimation. \begin{lemma} Under Assumption (ref), we have $\mathrm{var} (\widehat{\tau}_{\mathsf{unadj}}) = \dfrac{1}{n} \mathrm{var} \left\{ \dfrac{t}{\pi_{1}} y \right\}$ and \begin{align*} \mathrm{var} (\widetilde{\tau}_{\mathsf{adj}}) = \frac{1}{n} \left( \mathrm{var} \left\{ \frac{t}{\pi_{1}} y \right\} - \mathrm{var} \left\{ \left( \frac{t}{\pi_{1}} - 1 \right) (\mathbf{x} - \bm{\mu})^{\top} \bm{\beta} \right\} \right) \leq \mathrm{var} (\widehat{\tau}_{\mathsf{unadj}}). \end{align*} If $(\mathbf{x}_{i} - \bm{\mu})^{\top} \bm{\beta}$ in $\widetilde{\tau}_{\mathsf{adj}}$ is replaced by the true outcome regression function $\mathrm{b} (\cdot)$, the variance of $\widetilde{\tau}_{\mathsf{adj}}$ attains the \normalfont{SVB}. \end{lemma} See Appendix (ref) for the proof of Lemma (ref). When the dimension $p$ of the baseline covariates is small compared to the sample size $n$, one can compute a feasible estimator $\widehat{\tau}_{\mathsf{adj}, 1}^{\dag}$ by estimating $\bm{\beta}$ with ordinary least squares (OLS) between $t (y - \bar{\tau}) / \pi_{1}$ and $\mathbf{x} - \bm{\mu}$ with $\bm{\mu}$ replaced by $\bar{\mathbf{x}}$ and $\bar{\tau}$ replaced by $\widehat{\tau}_{\mathsf{unadj}}$: \begin{equation} \begin{split} \widehat{\bm{\beta}}_{c} & \coloneqq \left( \sum_{l = 1}^{n} (\mathbf{x}_{l} - \bar{\mathbf{x}}) (\mathbf{x}_{l} - \bar{\mathbf{x}})^{\top} \right)^{-1} \sum_{j = 1}^{n} (\mathbf{x}_{j} - \bar{\mathbf{x}}) \frac{t_{j} (y_{j} - \widehat{\tau}_{\mathsf{unadj}})}{\pi_{1}}, \\ \widehat{\tau}_{\mathsf{adj}, 1}^{\dag} & \coloneqq \widehat{\tau}_{\mathsf{unadj}} - \frac{1}{n} \sum_{i = 1}^{n} \left( \frac{t_{i}}{\pi_{1}} - 1 \right) (\mathbf{x}_{i} - \bar{\mathbf{x}})^{\top} \widehat{\bm{\beta}}_{c} \equiv \widehat{\tau}_{\mathsf{unadj}} - \frac{1}{n} \sum_{i = 1}^{n} \sum_{j = 1}^{n} \left( \frac{t_{i}}{\pi_{1}} - 1 \right) H_{i, j} \frac{t_{j} (y_{j} - \widehat{\tau}_{\mathsf{unadj}})}{\pi_{1}}. \end{split} \end{equation} We also define an alternative adjusted estimator not centering $y$ that will appear later: \begin{equation} \widehat{\tau}_{\mathsf{adj}, 1} \coloneqq \widehat{\tau}_{\mathsf{unadj}} - \frac{1}{n} \sum_{i = 1}^{n} \sum_{j = 1}^{n} \left( \frac{t_{i}}{\pi_{1}} - 1 \right) H_{i, j} \frac{t_{j} y_{j}}{\pi_{1}}. \end{equation} Written in the form of (ref) or (ref), one can view the adjusted estimator by linear working models as augmenting the unadjusted estimator with a second-order $V$-statistic. The corresponding least square regression coefficients $\widehat{\bm{\beta}}$ is defined similarly to $\widehat{\bm{\beta}}_{c}$ except for not centering $y$. To directly see the potential negative impact of the augmented $V$-statistic, its mean is, under the design-based framework and CRE, $- \frac{\pi_{0}}{\pi_{1}} \frac{n - 2}{n (n - 1)} \sum_{i = 1}^{n} H_{i, i} y_{i} (1) = O \left( p / n \right),$ under certain regularity conditions (e.g. Assumption (ref) later) on $\mathbf{X}$ and $\mathbf{y}$. The above derivation explains why the usual adjusted estimator $\widehat{\tau}_{\mathsf{adj}, 1}^{\dag}$ or $\widehat{\tau}_{\mathsf{adj}, 1}$ may hurt statistical inference when $p$ is close to $n$. The bias of $\widehat{\tau}_{\mathsf{adj}, 1}^{\dag}$ can be similarly shown to also be of order $O (p / n)$. \subsection*{The Role of HOIF and a Review} To explain how the theory of HOIF directly leads to an improved estimator, we first notice that $\widehat{\tau}_{\mathsf{aipw}}$ is unbiased and also assume that $\widehat{\mathrm{b}}$ is “sufficiently independent” from the sample $\mathbf{O}$. In the derivation below, we might write 0 redundantly as $\frac{1}{\pi_{1}} - \frac{1}{\pi_{1}}$: \begin{equation*} \mathrm{bias} (\widehat{\tau}_{\mathsf{aipw}}) = {\mathbb{E}} \left[ t \left( \frac{1}{\pi_{1}} - \frac{1}{\pi_{1}} \right) \left( \mathrm{b} (\mathbf{x}) - \widehat{\mathrm{b}} (\mathbf{x}) \right) \right] \equiv {\mathbb{E}} \left[ \pi_{1}^{\frac{1}{2}} \left( \frac{1}{\pi_{1}} - \frac{1}{\pi_{1}} \right) \pi_{1}^{\frac{1}{2}} \left( \mathrm{b} (\mathbf{x}) - \widehat{\mathrm{b}} (\mathbf{x}) \right) \right]. \end{equation*} The HOIF theory, in a nutshell, is to approximate the above bias of $\widehat{\tau}_{\mathsf{aipw}}$ by first choosing a set of $k$-dimensional transformations of $\mathbf{x}$, $\bar{\phi}_{k}(\mathbf{x}) = (\phi_{1} (\mathbf{x}), \cdots, \phi_{k} (\mathbf{x}))^{\top}$. $\bar{\phi}_{k}$ is often chosen by some background knowledge on the space the residual $\mathrm{b} - \widehat{\mathrm{b}}$ may lie in. Here we simply take $\phi$ as $\bar{\phi}_{k}(\mathbf{x}) \equiv \mathbf{x} - \bm{\mu}$. Next, we project the two residuals above, $\mathsf{res}_{1} \coloneqq \pi_{1}^{1 / 2} \left( \frac{1}{\pi_{1}} - \frac{1}{\pi_{1}} \right)$ and $\mathsf{res}_{2} \coloneqq \pi_{1}^{1 / 2} \left( \mathrm{b} (\mathbf{x}) - \widehat{\mathrm{b}} (\mathbf{x}) \right)$ onto the linear space spanned by $\pi_{1}^{1 / 2} (\mathbf{x} - \bm{\mu})$: \begin{align*} \widetilde{\mathsf{res}}_{1} & \coloneqq \pi_{1}^{1 / 2} (\mathbf{x} - \bm{\mu})^{\top} \left\{ {\mathbb{E}} [\pi_{1} (\mathbf{x} - \bm{\mu}) (\mathbf{x} - \bm{\mu})^{\top}] \right\}^{-1} {\mathbb{E}} \left[ \pi_{1} (\mathbf{x} - \bm{\mu}) \left( \frac{1}{\pi_{1}} - \frac{1}{\pi_{1}} \right) \right], \\ \widetilde{\mathsf{res}}_{2} & \coloneqq \pi_{1}^{1 / 2} (\mathbf{x} - \bm{\mu})^{\top} \left\{ {\mathbb{E}} [\pi_{1} (\mathbf{x} - \bm{\mu}) (\mathbf{x} - \bm{\mu})^{\top}] \right\}^{-1} {\mathbb{E}} \left[ \pi_{1} (\mathbf{x} - \bm{\mu}) (\mathrm{b} (\mathbf{x}) - \widehat{\mathrm{b}} (\mathbf{x})) \right]. \end{align*} The weight $\pi_{1}^{1 / 2}$ is chosen to ensure that the unknown $\mathrm{b}$ appeared in the above expectation can be replaced by the observed $y$. Armed with the above projections of the residuals, we can decompose $\mathrm{bias} (\widehat{\tau}_{\mathsf{aipw}})$ into two components by Pythagorean theorem: \begin{equation*} \begin{split} \mathrm{bias} (\widehat{\tau}_{\mathsf{aipw}}) \equiv & \ {\mathbb{E}} \left[ \mathsf{res}_{1} \cdot \mathsf{res}_{2} \right] = \widetilde{\mathrm{bias}} (\widehat{\tau}_{\mathsf{aipw}}) + \widetilde{\mathrm{bias}}^{\perp} (\widehat{\tau}_{\mathsf{aipw}}) \equiv 0, \text{ where } \\ \widetilde{\mathrm{bias}} (\widehat{\tau}_{\mathsf{aipw}}) \coloneqq & \ {\mathbb{E}} \left[ \widetilde{\mathsf{res}}_{1} \cdot \widetilde{\mathsf{res}}_{2} \right] \equiv 0 \\ \equiv & \ {\mathbb{E}} \left[ \pi_{1} \left( \frac{1}{\pi_{1}} - \frac{1}{\pi_{1}} \right) (\mathbf{x} - \bm{\mu})^{\top} \right] \left\{ {\mathbb{E}} [\pi_{1} (\mathbf{x} - \bm{\mu}) (\mathbf{x} - \bm{\mu})^{\top}] \right\}^{-1} {\mathbb{E}} \left[ (\mathbf{x} - \bm{\mu}) \pi_{1} (\mathrm{b} (\mathbf{x}) - \widehat{\mathrm{b}} (\mathbf{x})) \right] \\ = & \ {\mathbb{E}} \left[ \left( \frac{t}{\pi_{1}} - 1 \right) (\mathbf{x} - \bm{\mu})^{\top} \right] \left\{ {\mathbb{E}} \left[ (\mathbf{x} - \bm{\mu}) (\mathbf{x} - \bm{\mu})^{\top} \right] \right\}^{-1} {\mathbb{E}} \left[ (\mathbf{x} - \bm{\mu}) \frac{t (y - \widehat{\mathrm{b}} (\mathbf{x}))}{\pi_{1}} \right]. \end{split} \end{equation*} Setting $\widehat{\mathrm{b}} \equiv 0$, so $\widehat{\tau}_{\mathsf{unadj}} \equiv \widehat{\tau}_{\mathsf{aipw}}$, the gist of the HOIF theory is to estimate the “projected bias” $\widetilde{\mathrm{bias}} (\widehat{\tau}_{\mathsf{aipw}}) \equiv \widetilde{\mathrm{bias}} (\widehat{\tau}_{\mathsf{unadj}})$ by a \emph{second-order $U$-statistic} as follows: \begin{align} \widehat{\mathbb{IF}}_{\mathsf{unadj}, 2, 2} & \coloneqq \frac{1}{n (n - 1)} \sum_{1 \leq i \neq j \leq n} \left( \frac{t_{i}}{\pi_{1}} - 1 \right) (\mathbf{x}_{i} - \bar{\mathbf{x}})^{\top} \left\{ \frac{1}{n - 1} \sum_{l = 1}^{n} (\mathbf{x}_{l} - \bar{\mathbf{x}}) (\mathbf{x}_{l} - \bar{\mathbf{x}})^{\top} \right\}^{-1} (\mathbf{x}_{j} - \bar{\mathbf{x}}) \frac{t_{j} y_{j}}{\pi_{1}} \nonumber \\ & \equiv \frac{1}{n} \sum_{1 \leq i \neq j \leq n} \left( \frac{t_{i}}{\pi_{1}} - 1 \right) H_{i, j} \frac{t_{j} y_{j}}{\pi_{1}}. \end{align} Here we adopt the $\widehat{\mathbb{IF}}$ notation first introduced in robins2008higher because $\widehat{\mathbb{IF}}_{\mathsf{unadj}, 2, 2}$ is an estimator of the \emph{second-order influence function} of $\widetilde{\mathrm{bias}} (\widehat{\tau}_{\mathsf{unadj}})$. With $\widehat{\mathbb{IF}}_{\mathsf{unadj}, 2, 2}$, one can construct the following covariate adjusted estimator: \begin{equation} \widehat{\tau}_{\mathsf{adj}, 2} \coloneqq \widehat{\tau}_{\mathsf{unadj}} - \widehat{\mathbb{IF}}_{\mathsf{unadj}, 2, 2}, \end{equation} which is the main estimator that we study in Section (ref). Though $\widehat{\mathbb{IF}}_{\mathsf{unadj}, 2, 2}$ and $\widehat{\tau}_{\mathsf{adj}, 2}$ are motivated under the superpopulation framework, the way we tacitly estimate the precision matrix $\left\{ {\mathbb{E}} \left[ (\mathbf{x} - \bm{\mu}) (\mathbf{x} - \bm{\mu})^{\top} \right] \right\}^{-1}$ in $\widehat{\mathbb{IF}}_{\mathsf{unadj}, 2, 2}$ happens to be the “correct” choice under the design-based framework, as we condition on $\mathbf{X}$. It is also worth noting that $\widehat{\tau}_{\mathsf{adj}, 2}$ can be viewed as a “diagonal/trace-free” version of $\widehat{\tau}_{\mathsf{adj}, 1}$. As for the variance of $\widehat{\tau}_{\mathsf{adj}, 2}$, following well-established statistical theory of HOIF estimators, we have $\mathrm{var} (\widehat{\tau}_{\mathsf{adj}, 2}) = O \left( \frac{1}{n} + \frac{p}{n^{2}} \right)$, where the first factor $1 / n$ is attributed to $\widehat{\tau}_{\mathsf{unadj}}$ and the second factor $p / n^{2}$ comes from the second-order $U$-statistic $\widehat{\mathbb{IF}}_{\mathsf{unadj}, 2, 2}$. When $p = o (n)$, covariate adjustment by $\widehat{\tau}_{\mathsf{adj}, 2}$ never hurts the asymptotic variance (in fact, in Theorem (ref), $\widehat{\tau}_{\mathsf{adj}, 2}$ may reduce the asymptotic variance if $p = o (n)$). This follows from the “guiding principle” liu2017semiparametric that if either the propensity score or the outcome regression is consistently estimated, correcting bias by adding the second-order $U$-statistic does not inflate the asymptotic variance. Even if $p = O (n)$, the variance maintains the parametric $1 / n$ rate. However, demonstrating that $\widehat{\tau}_{\mathsf{adj}, 2}$ actually improves efficiency for $p = O (n)$ requires more careful analysis, deferred to Section (ref). Intuitively, since $\widehat{\mathbb{IF}}_{\mathsf{unadj}, 2, 2}$ estimates a $L_{2} ({\mathbb{P}})$-projection of $\frac{t}{\pi_{1}} y$, efficiency gains are possible even when $p$ is close to $n$. This fact has been known and explicitly noted in liu2020nearly, liu2023hoif. We hope that this review stimulates more interest in the HOIF theory by those developing statistical methods on randomized experiments. \begin{remark} We briefly compare $\widehat{\mathbb{IF}}_{\mathsf{unadj}, 2, 2}$ with estimators from liu2020nearly, liu2020rejoinder, liu2023hoif. Their original proposal addresses observational studies with unknown propensity scores, using $(\mathbf{x}_{i} - \bar{\mathbf{x}})^{\top} \widehat{\bm{\Sigma}}_{1}^{-1} (\mathbf{x}_{j} - \bar{\mathbf{x}})$ instead of $H_{i, j}$ to protect against model misspecification, where \begin{equation} \widehat{\bm{\Sigma}}_{1} \coloneqq \sum_{l = 1}^{n} \frac{t_{l}}{\pi_{1}} (\mathbf{x}_{l} - \bar{\mathbf{x}}) (\mathbf{x}_{l} - \bar{\mathbf{x}})^{\top} \end{equation} Since we consider CRE, the form in (ref) is preferable. In fact, it is much more difficult to analyze the statistical properties of $\widehat{\tau}_{\mathsf{adj}, 2}$ under observational studies and the superpopulation framework. There one has to control the difference between the sample and population precision matrices. Without imposing strong structural assumptions on the propensity score and outcome regression, higher-order $U$-statistics are needed to remove the bias due to estimating the precision matrix of a very large dimension. Finally, we discuss the choice of the transformation $\bar{\phi}_{k}$, largely ignored above. Choosing $\bar{\phi}_{k}(\mathbf{x}) = \mathbf{x} - \bm{\mu}$ makes $\widehat{\tau}_{\mathsf{adj}, 2}$ correct the \emph{own observation bias} of $\widehat{\tau}_{\mathsf{adj}, 1}$. If the analyst has better knowledge of the $\mathbf{x} \rightarrow y (1)$ mechanism, selecting nonlinear transformations of $\mathbf{x}$ better reflect this mechanism may further improve efficiency. \end{remark} \section{HOIF-Motivated Covariate Adjustment: Theoretical Results} To state our main theoretical results on $\widehat{\tau}_{\mathsf{adj}, 2}$, we need to further impose certain regularity conditions on the observed data $\mathbf{O}$. It is worth noting that to make our exposition more accessible, we choose the following easier-to-interpret regularity conditions on the data instead of the mathematically weaker conditions considered in the mainstream literature on design-based inference lei2021regression, lu2025debiased. \begin{assumption}[Regularity conditions on $\mathbf{O}$] The following regularity condition is occasionally imposed on the observed data $\mathbf{O}$: there exists an $n$-independent universal constant $B > 0$ such that $\max \left\{ \Vert y (1) \Vert_{\infty}, \frac{n}{p} \Vert H \Vert_{\infty} \right\} \leq B,$ and $\widehat{\bm{\Sigma}} \equiv \sum_{i = 1}^{n} (\mathbf{x}_{i} - \bar{\mathbf{x}}) (\mathbf{x}_{i} - \bar{\mathbf{x}})^{\top}$ is invertible. \end{assumption} \begin{remark} As will be seen in Appendix (ref), under Assumption (ref), by directly looking at the decomposition of variance formula into a sum of various “gadgets”, one can immediately tell the bias or the variance order (should be no greater than $O (1 / n)$) after covariate adjustment. In Appendix (ref), we will explain that the statistical orders of the “gadgets” can be easily deduced based on simple statistical intuition. \end{remark} \subsection{Statistical Properties of HOIF-Motivated Estimators in RCT} We state our main theoretical results under CRE in the theorem below. The corresponding results for the Bernoulli sampling will be commented in Remark (ref) that follows. \begin{theorem} Under CRE (Assumption (ref)) and the design-based framework, we have the following theoretical guarantees on the HOIF estimator $\widehat{\tau}_{\mathsf{adj}, 2}$: \begin{enumerate}[leftmargin=1cm,topsep=0.25pt] • The bias of $\widehat{\tau}_{\mathsf{adj}, 2}$ has the following form: \begin{equation} \begin{split} \mathrm{bias}^{\mathsf{d}} \left( \widehat{\tau}_{\mathsf{adj}, 2} \right) & = \frac{\pi_{0}}{\pi_{1}} \frac{1}{n (n - 1)} \sum_{1 \leq i \neq j \leq n} H_{i, j} y_{j} (1) = - \frac{\pi_{0}}{\pi_{1}} \frac{1}{n (n - 1)} \sum_{i = 1}^{n} H_{i, i} y_{i} (1). \end{split} \end{equation} In addition, if Assumption (ref) holds, then $\mathrm{bias}^{\mathsf{d}} \left( \widehat{\tau}_{\mathsf{adj}, 2} \right) = O \left( \frac{\pi_{0}}{\pi_{1}} \frac{\alpha}{n} \right).$ • The variance of $\widehat{\tau}_{\mathsf{adj}, 2}$ has the following form: Under Assumption (ref), \begin{align} & \ \mathrm{var}^{\mathsf{d}} \left( \widehat{\tau}_{\mathsf{adj}, 2} \right) = \nu^{\mathsf{d}} + \mathsf{Rem} = O \left( \frac{1}{n} \left\{ 1 + \frac{p}{n} \right\} \right), \end{align} where the exact form of $\mathsf{Rem} = o (1 / n)$ can be deduced from the proof of this claim in Appendix (ref) and the main term $\nu^{\mathsf{d}}$ reads as follows: \begin{equation} \begin{split} \nu^{\mathsf{d}} \coloneqq & \ \underbrace{\left( \frac{\pi_{0}}{\pi_{1}} \right) \frac{1}{n} V_{n} \left[ y_{i} (1) - \sum_{j \neq i} H_{j, i} y_{j} (1) \right]}_{\eqqcolon \, \nu^{\mathsf{d}}_{1} \, = \, O \left( \frac{1}{n} \right)} \\ & + \underbrace{\left( \frac{\pi_{0}}{\pi_{1}} \right)^{2} \frac{1}{n} \left\{ \frac{1}{n} \sum_{i = 1}^{n} H_{i, i} (1 - H_{i, i}) y_{i} (1)^{2} + \frac{1}{n} \sum_{1 \leq i \neq j \leq n} H_{i, j}^{2} y_{i} (1) y_{j} (1) \right\}}_{\eqqcolon \, \nu^{\mathsf{d}}_{2} \, = \, O \left( \frac{1}{n} \frac{p}{n} \right)}. \end{split} \end{equation} When $p = o (n)$, $\nu^{\mathsf{d}}$ can be further simplified to $\nu^{\mathsf{d}}_{1} = \frac{\pi_{0}}{\pi_{1}} \frac{1}{n} V_{n} \left[ y_{i} (1) - \sum_{j \neq i} H_{j, i} y_{j} (1) \right].$ \end{enumerate} \end{theorem} According to the second assertion of Theorem (ref), if $p = o (n)$, $\widehat{\tau}_{\mathsf{adj}, 2}$ always attains smaller asymptotic variance (after scaled by $n$) than $\widehat{\tau}_{\mathsf{unadj}}$; however, if $p = O (n)$, then the \emph{iff condition} for $\widehat{\tau}_{\mathsf{adj}, 2}$ to enjoy improved asymptotic efficiency (after scaled by $n$) compared to $\widehat{\tau}_{\mathsf{unadj}}$ is simply \begin{equation} \begin{split} \nu^{\mathsf{d}} - \mathrm{var}^{\mathsf{d}} (\widehat{\tau}_{\mathsf{unadj}}) \leq 0 \Longleftrightarrow \mathrm{var}^{\mathsf{d}} (\widehat{\tau}_{\mathsf{unadj}}) - \nu^{\mathsf{d}}_{1} \geq \nu^{\mathsf{d}}_{2}. \end{split} \end{equation} Since the criterion for asymptotic efficiency improvement is fully characterized, one can easily conduct numerical experiments to determine whether $\widehat{\tau}_{\mathsf{adj}, 2}$ has smaller asymptotic variance than $\widehat{\tau}_{\mathsf{unadj}}$ for given baseline covariates $\mathbf{X}$ and potential outcomes $\mathbf{y} (1)$. \begin{remark}[Interpreting $\widehat{\tau}_{\mathsf{adj}, 2}$ as a leave-one-out regression adjustment estimator]\leavevmode The first term $\nu^{\mathsf{d}}_{1}$ of $\nu^{\mathsf{d}}$ in (ref) constitutes the sample variance of $y_{i} (1) - \sum_{j \neq i} H_{j, i} y_{j} (1)$ instead of $y_{i} (1)$, as in the variance of $\widehat{\tau}_{\mathsf{unadj}}$. The term being subtracted off, $\sum_{j \neq i} H_{j, i} y_{j} (1)$, can be represented as $\sum_{j \neq i} H_{j, i} y_{j} (1) \equiv (\mathbf{x}_{i} - \bar{\mathbf{x}})^{\top} \widehat{\bm{\beta}}_{-i},$ where $\widehat{\bm{\beta}}_{-i} \coloneqq \widehat{\bm{\Sigma}}^{-} \sum_{j \neq i} (\mathbf{x}_{j} - \bar{\mathbf{x}}) y_{j} (1)$, which is the leave-one-out coefficient estimator for the linear projection of the potential outcomes $y_{i} (1)$ onto the linear span of $\mathbf{x}_{i} - \bar{\mathbf{x}}$, resembling the construction in wu2018loop. In an ongoing work, we apply this idea of using leave-one-out regression adjustment in broader settings, where the outcome regression is fit by a working generalized linear model. \end{remark} \begin{remark} The proof of Theorem (ref) is in Appendix (ref). Under the Bernoulli sampling with $t_{i} \overset{\rm i.i.d.}{\sim} \mathrm{Bernoulli} (\pi_{1})$ for $i \in [n]$, we immediately have $\mathrm{bias}^{\mathsf{d}} (\widehat{\tau}_{\mathsf{adj}, 2}) = 0$. \end{remark} The next result shows that statistical guarantees of $\widehat{\tau}_{\mathsf{adj}, 2}$ parallel to those under the design-based framework in Theorem (ref) continue to hold under the i.i.d. superpopulation framework. We mention in passing that this result has been obtained in previous works by one of the authors of this paper liu2017semiparametric, liu2020nearly, liu2023hoif, so we omit the proof. \begin{proposition} Under CRE (Assumption (ref)) and the superpopulation framework, the following hold: \begin{enumerate}[leftmargin=1cm,topsep=0.25pt] • The bias of $\widehat{\tau}_{\mathsf{adj}, 2}$ has the following form: $\mathrm{bias} (\widehat{\tau}_{\mathsf{adj}, 2}) = {\mathbb{E}} (\widehat{\tau}_{\mathsf{adj}, 2} - \tau) = - \dfrac{\pi_{0}}{\pi_{1}} \dfrac{1}{n - 1} {\mathbb{E}} \left[ H_{1, 1} \mathrm{b} (\mathbf{x}_{1}) \right]$. In addition, if Assumption (ref) holds, then $\mathrm{bias} (\widehat{\tau}_{\mathsf{adj}, 2}) = O \left( \dfrac{\pi_{0}}{\pi_{1}} \dfrac{\alpha}{n} \right)$. • The variance of $\widehat{\tau}_{\mathsf{adj}, 2}$ has the following order: Under Assumption (ref), further suppose that $\widehat{\bm{\Sigma}}$ has bounded eigenvalues, $\mathrm{var} (\widehat{\tau}_{\mathsf{adj}, 2}) = O \left( \dfrac{1}{n} \left\{ 1 + \dfrac{p}{n} \right\} \right)$. \end{enumerate} \end{proposition} \begin{remark}[On the asymptotic distribution of $\widehat{\tau}_{\mathsf{adj}, 2}$] Under Assumptions (ref)--(ref) and the superpopulation framework, if we additionally suppose that there exists a constant $\sigma^{2} > 0$ such that \begin{equation} \lim_{n \rightarrow \infty} n \cdot \mathrm{var} (\widehat{\tau}_{\mathsf{adj}, 2}) \rightarrow \sigma^{2}, \end{equation} we have $\frac{\widehat{\tau}_{\mathsf{adj}, 2} - \tau}{\sqrt{\mathrm{var} (\widehat{\tau}_{\mathsf{adj}, 2})}} \rightsquigarrow N (0, 1).$ The Gaussian limiting distribution of $\widehat{\tau}_{\mathsf{adj}, 2}$ follows from three main steps: (i) showing that the bias of $\widehat{\tau}_{\mathsf{adj}, 2}$ is $o (1 / \sqrt{n})$; (ii) showing that the variance of $\widehat{\tau}_{\mathsf{adj}, 2}$ is $O (1 / n)$; and (iii) invoking the CLT of second-order $U$-statistics (Corollary 1.2 of bhattacharya1992class). We note that, as mentioned in numerous places in liu2020nearly, the asymptotic distribution of $\widehat{\tau}_{\mathsf{adj}, 2}$ can be obtained by applying Corollary 1.2 of bhattacharya1992class. In the design-based framework, (ref) can be replaced by \begin{equation} \lim_{n \rightarrow \infty} n \cdot \nu^{\mathsf{d}} \rightarrow \sigma^{2}, \end{equation} and one needs to further adapt the proof technique in bhattacharya1992class to L\'{e}vy's martingale CLT by adapting the proof technique in bhattacharya1992class, or using results in koike2023high as in lu2025debiased. Conditions (ref) and (ref) essentially require that the variance scaled by $n$ has a limit as $n \rightarrow \infty$; see Appendix (ref) for details. \end{remark} \begin{remark}[Semiparametric efficiency under the superpopulation framework] Under the superpopulation framework, when $p$ is fixed and one imposes smoothness assumption on $\mathrm{b}$, say H\"{o}lder smooth with smoothness index $s > 0$, then it is easy to see that $\widehat{\tau}_{\mathsf{adj}, 2}$, but with $\mathbf{x} - \bar{\mathbf{x}}$ replaced by $\bar{\phi}_{k}(\mathbf{x}) = (\phi_{1} (\mathbf{x}), \cdots, \phi_{k} (\mathbf{x}))$, where $\bar{\phi}_{k}(\cdot)$ denotes low-degree polynomial transformations up to degree $k \asymp \log \log n$, achieves the SVB. \end{remark} \subsection{A Variety of HOIF-Motivated Estimators} In the previous section, we have demonstrated that the HOIF theory motivates an adjusted estimator $\widehat{\tau}_{\mathsf{adj}, 2}$ that (1) reduces the bias of the adjusted estimator $\widehat{\tau}_{\mathsf{adj}, 1}^{\dag}$ and (2) has asymptotic variance no greater than $\widehat{\tau}_{\mathsf{unadj}}$ whenever $p = o (n)$, sometimes more efficient than $\widehat{\tau}_{\mathsf{unadj}}$ when $p = O (n)$ and $p < n$. We now propose other HOIF-motivated estimators in the vicinity of $\widehat{\tau}_{\mathsf{adj}, 2}$, with slightly different statistical properties. To motivate alternatives, we consider a special scenario where the potential outcomes are constant, i.e. $y_{i} (1) \equiv c, i \in [n]$. Without loss of generality, we take $c \equiv 1$. Here both $\widehat{\tau}_{\mathsf{unadj}} \equiv 1$ and $\widehat{\tau}_{\mathsf{adj}, 1}^{\dag} \equiv 1$ are \emph{error-free}, but $\widehat{\tau}_{\mathsf{adj}, 2}$ fails to be \emph{error-free} in the extreme scenario of homogeneous potential outcomes. To restore the \emph{error-free} property in such a case, one could remove the diagonal/trace part from the error-free adjusted estimator $\widehat{\tau}_{\mathsf{adj}, 1}^{\dag}$ defined in (ref) instead of $\widehat{\tau}_{\mathsf{adj}, 1}$ defined in (ref): \begin{equation} \widehat{\tau}_{\mathsf{adj}, 2}^{\dag} \coloneqq \widehat{\tau}_{\mathsf{unadj}} - \widehat{\mathbb{IF}}_{\mathsf{unadj}, 2, 2}^{\dag}, \text{ where } \widehat{\mathbb{IF}}_{\mathsf{unadj}, 2, 2}^{\dag} \coloneqq \frac{1}{n} \sum_{1 \leq i \neq j \leq n} \left( \frac{t_{i}}{\pi_{1}} - 1 \right) H_{i, j} \frac{t_{j} (y_{j} - \widehat{\tau}_{\mathsf{unadj}})}{\pi_{1}}. \end{equation} It is natural to conjecture that the asymptotic variance of $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$ should be the same as that of $\widehat{\tau}_{\mathsf{adj}, 2}$ in (ref), except that $y_{i} (1)$ is replaced by $y_{i} (1) - \bar{\tau}$ for all $i \in [n]$, i.e. \begin{align*} & \mathrm{var}^{\mathsf{d}} (\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}) = \ \left( \frac{\pi_{0}}{\pi_{1}} \right) \frac{1}{n} V_{n} \left[ (y_{i} (1) - \bar{\tau}) - \sum_{j \neq i} H_{j, i} (y_{j} (1) - \bar{\tau}) \right] \\ & + \left( \frac{\pi_{0}}{\pi_{1}} \right)^{2} \frac{1}{n} \left\{ \frac{1}{n} \sum_{i = 1}^{n} H_{i, i} (1 - H_{i, i}) (y_{i} (1) - \bar{\tau})^{2} + \frac{1}{n} \sum_{1 \leq i \neq j \leq n} H_{i, j}^{2} (y_{i} (1) - \bar{\tau}) (y_{j} (1) - \bar{\tau}) \right\} + o (n^{-1}). \end{align*} Despite being \emph{error-free} in the case of homogeneous potential outcomes, similar to $\widehat{\tau}_{\mathsf{adj}, 2}$, $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$ is not unbiased under CRE. If one insists on constructing a \emph{bias-free} adjusted estimator under CRE, the following estimator can be constructed, again building on $\widehat{\tau}_{\mathsf{adj}, 2}$: \begin{equation} \widehat{\tau}_{\mathsf{adj}, 3} \coloneqq \widehat{\tau}_{\mathsf{adj}, 2} + \frac{\pi_{0}}{\pi_{1}} \frac{1}{n (n - 1)} \sum_{i = 1}^{n} H_{i, i} \frac{t_{i} y_{i}}{\pi_{1}}. \end{equation} In Appendix (ref), we show that $\widehat{\tau}_{\mathsf{adj}, 3}$ is unbiased under CRE, but biased under the Bernoulli sampling. The following proposition characterizes the asymptotic variance of $\widehat{\tau}_{\mathsf{adj}, 3}$. Of course, a similar strategy can be employed to also completely remove the bias of $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$, which we do not further pursue. \begin{proposition} Under Assumptions (ref) -- (ref), the following holds: $\mathrm{var}^{\mathsf{d}} (\widehat{\tau}_{\mathsf{adj}, 3}) = \mathrm{var}^{\mathsf{d}} (\widehat{\tau}_{\mathsf{adj}, 2}) + o \left( \frac{1}{n} \right).$ \end{proposition} In words, one could completely remove the bias of $\widehat{\tau}_{\mathsf{adj}, 2}$ due to CRE \emph{for free asymptotically}. The proof can be found in Appendix (ref). The conclusion of Proposition (ref) directly implies the asymptotic normality of $\widehat{\tau}_{\mathsf{adj}, 3}$ under the same conditions as in Remark (ref). Using a similar strategy, one can also construct an exactly unbiased version of $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$ under CRE, denoted as $\widehat{\tau}_{\mathsf{adj}, 3}^{\dag}$: \begin{align*} \widehat{\tau}_{\mathsf{adj}, 3}^{\dag} \coloneqq \widehat{\tau}_{\mathsf{adj}, 2}^{\dag} - 2 \frac{\pi_{0}}{\pi_{1}} \left( 1 - \frac{\pi_{0}}{\pi_{1}} \frac{1}{n - 1} \right) \frac{1}{n - 2} \left\{ \frac{p}{n} \widehat{\tau}_{\mathsf{unadj}} - \frac{1}{n} \sum_{i = 1}^{n} H_{i, i} \frac{t_{i} y_{i}}{\pi_{1}} \right\}. \end{align*} The reason for $\widehat{\tau}_{\mathsf{adj}, 3}^{\dag}$ will become immediately clear once we reveal the bias of $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$ under CRE in Proposition (ref) in the next subsection. \subsection{HOIF-Motivated Estimators: A Unifying Theme of Recent Proposals} As mentioned in the Introduction, HOIF-motivated estimators unify several recently proposed adjusted estimators. Unlike the OLS-based adjusted estimator $\widehat{\tau}_{\mathsf{adj}, 1}^{\dag}$ or $\widehat{\tau}_{\mathsf{adj}, 1}$, these estimators are CAN and has guarantee efficiency gains over $\widehat{\tau}_{\mathsf{unadj}}$ even when $p \gtrsim \sqrt{n}$. Here, we will illustrate that HOIF-motivated estimators unify those from: lei2021regression, lu2025debiased, chang2024exact, and also jiang2025adjustments. We start with lu2025debiased, in which the following debiased adjusted estimator of $\bar{\tau}$ was proposed: \begin{equation} \widehat{\tau}_{\mathsf{db}} \coloneqq \widehat{\tau}_{\mathsf{adj}, 1}^{\dag} + \frac{\pi_{0}}{\pi_{1}} \frac{1}{n} \sum_{i = 1}^{n} \frac{t_{i}}{\pi_{1}} H_{i,i} \left( y_{i} - \widehat{\tau}_{\mathsf{unadj}} \right). \end{equation} As pointed out in lu2025debiased, the estimator proposed in chang2024exact is similar to $\widehat{\tau}_{\mathsf{db}}$ but is exactly unbiased under CRE, and thus we denote the estimator in chang2024exact as $\widehat{\tau}_{\mathsf{db}}^{u}$. The following lemma immediately classifies $\widehat{\tau}_{\mathsf{db}}$ and $\widehat{\tau}_{\mathsf{db}}^{u}$ as HOIF-motivated estimators. The proof is deferred to Appendix (ref). \begin{lemma} The following algebraic equivalences hold: $\widehat{\tau}_{\mathsf{db}} \equiv \widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$ and $\widehat{\tau}_{\mathsf{db}}^{u} \equiv \widehat{\tau}_{\mathsf{adj}, 3}^{\dag}$. \end{lemma} \begin{remark} The sole difference between $\widehat{\tau}_{\mathsf{adj}, 2}$ and $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$, and now also $\widehat{\tau}_{\mathsf{db}}$, is that the former does not center $y_{i}$ by $\widehat{\tau}_{\mathsf{unadj}}$. In a slightly different vein, one can also consider $\widehat{\tau}_{\mathsf{unadj}}$ as an “adjusted” estimator using the linear working model with \emph{only the intercept term but no baseline covariates}: $\widehat{\tau}_{\mathsf{unadj}} \equiv \frac{1}{n} \sum_{i = 1}^{n} \frac{t_{i}}{\pi_{1}} \left( y_{i} - \widehat{\tau}_{\mathsf{unadj}} \right) + \widehat{\tau}_{\mathsf{unadj}}.$ Then by a similar line of reasoning to that of Section (ref), we should use a second-order $U$-statistic to estimate the following “projected bias”: \begin{align*} {\mathbb{E}} \left[ \left( \frac{t}{\pi_{1}} - 1 \right) (\mathbf{x} - \bm{\mu})^{\top} \right] (n^{-1} \bm{\Sigma})^{-} {\mathbb{E}} \left[ (\mathbf{x} - \bm{\mu}) \frac{t}{\pi_{1}} (y - \widehat{\tau}_{\mathsf{unadj}}) \right], \end{align*} leading to the statistic $\widehat{\mathbb{IF}}_{\mathsf{unadj}, 2, 2}^{\dag}$ defined in (ref). We hope that this algebraic equivalence sheds some new light on the bias & variance reduction “mechanism” of $\widehat{\tau}_{\mathsf{db}}$ to readers. \end{remark} Next, we present the statistical properties of $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$, and equivalently, $\widehat{\tau}_{\mathsf{db}}$. \begin{proposition} Under Assumptions (ref)--(ref), we have the following theoretical guarantees on $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$ and equivalently $\widehat{\tau}_{\mathsf{db}}$: \begin{enumerate}[leftmargin=1cm,topsep=0.25pt] • The design-based bias has the following form: \begin{equation} \begin{split} \mathrm{bias}^{\mathsf{d}} \left( \widehat{\tau}_{\mathsf{adj}, 2}^{\dag} \right) & \equiv \mathrm{bias}^{\mathsf{d}} \left( \widehat{\tau}_{\mathsf{db}} \right) = 2 \frac{\pi_{0}}{\pi_{1}^{2}} \frac{n_{1} - 1}{n - 1} \frac{1}{n - 2} \left\{ \frac{p}{n} \bar{\tau} + \frac{1}{n} \sum_{1 \leq i \neq j \leq n} H_{i, j} y_{j} (1) \right\} \\ & = 2 \frac{\pi_{0}}{\pi_{1}} \left( 1 - \frac{\pi_{0}}{\pi_{1}} \frac{1}{n - 1} \right) \frac{1}{n - 2} \left\{ \frac{p}{n} \bar{\tau} - \frac{1}{n} \sum_{i = 1}^{n} H_{i, i} y_{i} (1) \right\}. \end{split} \end{equation} In addition, if Assumption (ref) holds, then $\mathrm{bias}^{\mathsf{d}} \left( \widehat{\tau}_{\mathsf{adj}, 2}^{\dag} \right) \equiv \mathrm{bias}^{\mathsf{d}} \left( \widehat{\tau}_{\mathsf{db}} \right) = O \left( \frac{\pi_{0}}{\pi_{1}} \frac{\alpha}{n} \right).$ • The design-based variance has the following approximation under Assumption (ref): \begin{equation} \mathrm{var}^{\mathsf{d}} \left( \widehat{\tau}_{\mathsf{adj}, 2}^{\dag} \right) \equiv \mathrm{var}^{\mathsf{d}} \left( \widehat{\tau}_{\mathsf{db}} \right) = \nu^{\mathsf{d}\dag} + o (n^{-1}), \end{equation} \begin{align*} \text{where } & \nu^{\mathsf{d}\dag} \coloneqq \ \underbrace{\left( \frac{\pi_{0}}{\pi_{1}} \right) \frac{1}{n} V_{n} \left[ (y_{i} (1) - \bar{\tau}) - \sum_{j \neq i} H_{j, i} (y_{j} (1) - \bar{\tau}) \right]}_{\eqqcolon \, \nu^{\mathsf{d}}_{\mathsf{db}, 1}} \\ & + \underbrace{\left( \frac{\pi_{0}}{\pi_{1}} \right)^{2} \frac{1}{n} \left\{ \frac{1}{n} \sum_{i = 1}^{n} H_{i, i} (1 - H_{i, i}) (y_{i} (1) - \bar{\tau})^{2} + \frac{1}{n} \sum_{1 \leq i \neq j \leq n} H_{i, j}^{2} (y_{i} (1) - \bar{\tau}) (y_{j} (1) - \bar{\tau}) \right\}}_{\eqqcolon \, \nu^{\mathsf{d}}_{\mathsf{db}, 2}}. \end{align*} \end{enumerate} \end{proposition} The proof of Proposition (ref) can be found in Appendix (ref). The asymptotic normality of $\widehat{\tau}_{\mathsf{db}}$ can be obtained in a fashion similar to that of $\widehat{\tau}_{\mathsf{adj}, 2}$, and therefore omitted. It is now also clear why $\widehat{\tau}_{\mathsf{adj}, 3}^{\dag}$ is exactly unbiased under CRE. We also refer to lu2025debiased for analogous results on ATE. Unlike lu2025debiased, we, in fact, characterize the exact variance of $\widehat{\tau}_{\mathsf{db}}$ or $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$, and the asymptotic variance is simply a corollary. We do not further consider the statistical properties of $\widehat{\tau}_{\mathsf{db}}$ under the superpopulation framework. \begin{remark} We briefly compare the variances of $\widehat{\tau}_{\text{adj}, 2}$ and $\widehat{\tau}_{\text{adj}, 2}^{\dag}$ (equivalent to $\widehat{\tau}_{\text{db}}$). Equations (ref) and (ref) show that their asymptotic variances differ only in whether potential outcomes are centered by $\bar{\tau}$. The advantage of $\widehat{\tau}_{\text{adj}, 2}^{\dag}$ emerges when potential outcomes are homogeneous with a small coefficient of variation. Although centering is common in practice, Appendix (ref) provides computer-assisted examples where centering increases asymptotic variance. \end{remark} We are left to discuss the connection of the estimators proposed in lei2021regression and jiang2025adjustments to HOIF-motivated estimators. The estimator proposed in lei2021regression, denoted by $\widehat{\tau}_{\rm ld}$, differs from $\widehat{\tau}_{\mathsf{db}}$ (ref) only in the way $\widehat{\bm{\beta}}_{c}$ in $\widehat{\tau}_{\mathsf{adj}, 1}^{\dag}$ (see (ref)) is computed. Instead of using $\widehat{\bm{\Sigma}}$ in $\widehat{\bm{\beta}}_{c}$, $\widehat{\tau}_{\rm ld}$ estimates $\widehat{\Sigma}$ using covariates only in the treated group, although the covariate distributions in the two groups should be the same by design. $\widehat{\tau}_{\rm ld}$ remains CAN if $p = O (n^{2 / 3})$ up to a log-factor lei2021regression, and the stronger dependence on the dimension $p$ is a result of the mismatch between using an estimated $\widehat{\bm{\Sigma}}_{1}$ (ref) in $\widehat{\tau}_{\mathsf{adj}, 1}^{\dag}$ and using the true $\widehat{\bm{\Sigma}}$ when removing the “diagonal”. The extra term of $\widehat{\tau}_{\rm ld}$ over $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$ contains an error of order $\Vert \widehat{\bm{\Sigma}}^{-1} - \widehat{\bm{\Sigma}}_{1}^{-1} \Vert_{\mathrm{op}} = O_{{\mathbb{P}}} (p^{1 / 2} / n^{3 / 2})$ up to log-factors by matrix Bernstein inequality tropp2015introduction, resulting in $p = O (n^{2 / 3})$; see the end of Appendix (ref) for derivations. Finally, jiang2025adjustments directly uses the $V$-statistic $\widehat{\tau}_{\mathsf{adj}, 1}^{\dag}$, because they either assume $p = o (\sqrt{n})$ or assume that the linear working model for the outcome regression is correctly specified. We now summarize our findings in this section in Table (ref) below, which, in our opinion, is the most important message in our paper. \begin{table}[htbp] \caption{The connection between HOIF-motivated estimators and other recent proposals. Here $\widetilde{o}(\cdot)\coloneqq o(\cdot/(\log n)^{1/3})$. In the first row, the extra term between $\widehat{\tau}_{\rm ld}$ and $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$ can be found in the end of Appendix (ref).} \begin{tabular}{cccc} \Xhline{3\arrayrulewidth} \textbf{Estimator} & \textbf{Correspondence} & \textbf{Linear Model} & \textbf{Dimension} \\ \Xhline{3\arrayrulewidth} \addlinespace \multirowcell{2}{$\widehat{\tau}_{\mathrm{ld}}$ \\ {\scriptsizelei2021regression}} & \multirowcell{2}{$\equiv \widehat{\tau}_{\mathrm{adj}, 2}^{\dag} + \mathrm{rem}$} & \multirowcell{2}{\ding{55}} & \multirowcell{2}{$p = \widetilde{o} \left( n^{2/3} \right)$} \\ [20pt] \hline \addlinespace \addlinespace \multirowcell{2}{$\widehat{\tau}_{\mathrm{db}}$ \\ {\scriptsizelu2025debiased}} & \multirowcell{2}{$\equiv \widehat{\tau}_{\mathrm{adj}, 2}^{\dag}$} & \multirowcell{2}{\ding{55}} & \multirowcell{2}{$p = o(n)$} \\left[20pt] \hline \addlinespace \multirowcell{2}{$\widehat{\tau}_{\mathrm{db}}^{u}$ \\ {\scriptsizechang2024exact}} & \multirowcell{2}{$\equiv \widehat{\tau}_{\mathrm{adj}, 3}^{\dag}$} & \multirowcell{2}{\ding{55}} & \multirowcell{2}{$p = o(n)$} \\left[20pt] \hline \addlinespace \multirowcell{3}{$\widehat{\tau}_{\mathrm{adj}, 1}^{\dag}$ \\ {\scriptsizema2022regression, ye2023toward} \\ {\scriptsizejiang2025adjustments, and etc.}} & \multirowcell{3}{\shortstack{ V-statistic \\ version of $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$}} & \multirowcell{3}{\shortstack{\ding{51} \\ \\ \ding{55}}} & \multirowcell{3}{\shortstack{$p = o(n)$ \\ $p = o(\sqrt{n})$}} \\left[30pt] \addlinespace \Xhline{3\arrayrulewidth} \end{tabular} \end{table} \begin{remark} We finally comment on the connection of our work with two recent preprints: abadie2025unbiased and song2025neumann. abadie2025unbiased use the ridge inverse $\widehat{\bm{\Sigma}}_{\lambda}^{-1} = (\mathbf{X}^{\top} \mathbf{X} + \lambda \mathbf{I})^{-1}$ instead of the non-regularized inverse $\widehat{\bm{\Sigma}}^{-1}$ to improve the numerical stability of the estimator. In the implementation of our proposed estimator, we in fact use the generalized inverse, which is also numerically stable and can be viewed as the limit of $\widehat{\bm{\Sigma}}_{\lambda}^{-1}$ as $\lambda \rightarrow 0$. From our analysis, the adjusted estimator will always have negligible bias compared to the sampling variability if a $U$-statistic is used instead of a $V$-statistic and $\mathbf{t}$ is not used to compute the inverse Gram matrix. Therefore, we conjecture that as long as $\lambda$ is appropriately chosen such that $\widehat{\bm{\Sigma}}_{\lambda}^{-1} - \widehat{\bm{\Sigma}}^{-1}$ is sufficiently close, we can conclude that the ridge-regularized adjusted estimator of abadie2025unbiased has an asymptotic variance sometimes smaller than and never greater than that of $\widehat{\tau}_{\mathsf{unadj}}$. Regarding the estimator proposed in song2025neumann, since they estimate $\widehat{\bm{\Sigma}}^{-1}$ by $\widehat{\bm{\Sigma}}_{t}^{-1}$ as in lei2021regression, as already suggested in liu2020nearly, additional corrections of the estimation bias is needed to relax the requirement on $p$, as done in liu2023hoif or initially in liu2020nearly. It should be noted that song2025neumann focused on the design-based framework, so the true Gram matrix $\widehat{\bm{\Sigma}}$ is known and bias correction using Neumann series is simplified compared to liu2023hoif or liu2020nearly, because they assume the existence of a (possibly fictitious) superpopulation so the true $\bm{\Sigma}$ is unknown. \end{remark} \section{Simulation Studies and Real Data Application} In this section, we conduct numerical experiments and real data analysis to explore our theoretical findings regarding the design-based bias and variance formulas for $\widehat{\tau}_{\mathsf{unadj}}$, $\widehat{\tau}_{\mathsf{adj}, 2}$, $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$ and $\widehat{\tau}_{\mathsf{adj}, 3}$. The accompanying \texttt{R} code for replicating our simulations is available at \href{https://github.com/Cinbo-Wang/HOIF-Car}{the GitHub repository}, which includes both the exact bias and variance formulas, as well as the asymptotic versions presented in our paper. \subsection{Simulation studies} We first generate a very large data matrix $\mathcal{X} \in {\mathbb{R}}^{N\times N}$ with $N = 5000$, with each row $\mathcal{X}_i \sim t_3(0,\Sigma)$, where $\Sigma_{k,l}=0.1^{|k-l|},k,l\in[N]$. Then we generated two types of exogenous noise $\mathcal{\epsilon}_0\in {\mathbb{R}}^{N}$: univariate $t_3$ distribution and “worst-case residual” lei2021regression): $\mathcal{\epsilon}_0=\text{Scale}\left((\mathbb{I}_N - \mathbf{H})(H_{1,1},\cdots,H_{N,N})^\top\right),$ where $\text{Scale}(a_i)\coloneqq \left(a_i-\bar{a}\right)/\left(\sum_{i=1}^{n}(a_i-\bar{a})^2\right)^{1/2}$. Also, we let $\bm{\beta}_j = (-1)^j/\sqrt{j},j\in[N]$. We vary the following in our experiments: \begin{itemize}[leftmargin=0.5cm,topsep=0.25pt] • \textbf{covariates} $\mathbf{X}$: sample size $n \in \{ 50,100,500,1000\}$, and covariate dimension $p = \lceil n \cdot \alpha \rceil$ where $\alpha \in \{ 0.05,\dots,0.7\}$. Then we choose $\mathbf{X} \coloneqq \mathcal{X}_{1:n,1:p} \in {\mathbb{R}}^{n \times p}$. • \textbf{potential outcome model}: linear $f(\mathbf{x})=\mathbf{x}^{\top} \bm{\beta}_{[p]}$, and nonlinear $f(\mathbf{x})=\text{sign}\left(\mathbf{x}^{\top} \bm{\beta}_{[p]}\right)|\mathbf{x}^{\top} \bm{\beta}_{[p]}|^{\frac{1}{2}} + \sin(\mathbf{x}^{\top} \bm{\beta}_{[p]})$. Here, $\bm{\beta}_{[p]}$ represents the first $p$ elements of $\bm{\beta}$. • \textbf{exogenous error distribution} $\epsilon$: for both $t_3$ and “worst-case residual” types, we choose the first $n$ elements in $\epsilon_0$. • \textbf{signal size} $ \gamma\in \{1,2\}$: the potential outcomes are generated according to $y_i (1) = 1 + f(\mathbf{x}_i) + \epsilon_i\cdot \sqrt{V_n\left(f(\mathbf{X})\right)/V_n(\epsilon)} / \sqrt{\gamma}$ for $i\in [n]$. • \textbf{treatment assignment ratio} $\pi_1 \in \{\frac{1}{2},\frac{2}{3}\}$: once the pre-treatment variables $\{\left(\mathbf{x}_i,y_i(1)\right)\}_{i=1}^{n}$ are generated, we fix them and randomly assign $\lceil n\pi_1\rceil$ samples to the treatment arm. (When $\pi_1=\frac{2}{3}$, the “true” $\pi_1$ we used for estimation is $\lceil n\pi_1\rceil / n $.) \end{itemize} The relative efficiencies of different estimators, based on their exact and approximate variances, benchmarked by $\mathrm{var}^{\mathsf{d}} (\widehat{\tau}_{\mathsf{unadj}})$, are presented in Figure (ref). Overall, $\widehat{\tau}_{\mathsf{adj}, 2}$ and $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$($\widehat{\tau}_{\mathsf{db}}$) exhibit very similar relative efficiencies; however, in the nonlinear outcome model setting, there are cases where $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$ demonstrates greater efficiency. As indicated in Remark (ref), $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$ is expected to be more efficient when CoV$^{2}$ is small, which is confirmed by the results shown in Figure (ref). \begin{figure}[H] \caption{Relative efficiencies of $\widehat{\tau}_{\mathsf{adj}, 2}^{\dag}$ (or equivalently $\widehat{\tau}_{\mathsf{db}}$), $\widehat{\tau}_{\mathsf{adj}, 2}$, $\widehat{\tau}_{\mathsf{adj}, 3}$ based on exact and approximate formula.} \end{figure} \begin{figure}[H] \caption{CoV$^{2}$ vs. $\mathrm{var}^{\mathsf{d}} (\widehat{\tau}_{\mathsf{adj}, 2}) / \mathrm{var}^{\mathsf{d}} (\widehat{\tau}_{\mathsf{adj}, 2}^{\dag})$. Each point represents a particular simulation setting in Figure (ref); only the settings with CoV$^{2} \leq 100$ are shown here.} \end{figure} Furthermore, we conduct realistic simulations by actually drawing treatment assignments. Once $\{\mathbf{x}_i, y_i(1)\}_{i = 1}^{n}$ are generated, they are fixed, and random treatment assignments are drawn repeatedly from CRE, with the Monte Carlo repetition size $K = 20000$. The resulting biases and variances are shown in the Appendix (ref) and the overall message is similar. \subsection{Real Data Application} We next apply our proposed estimators to data from the NIDA-CTN-0030 trial, which tested whether adding individual drug counseling to buprenorphine/naloxone treatment and standard medical management (SMM) improved outcomes for prescription opioid dependence. We focused on the first phase, where patients were randomized to either SMM or enhanced medical management (EMM), stratified by heroin use and chronic pain history. For simplicity, we treated the design as CRE by including the randomization stratum as baseline covariates. The outcome of interest is the proportion of positive urine laboratory results among all tests. We include the following baseline covariates: randomization stratum, age, sex, and baseline urine test results. Missing data were replaced by medians. Following wang2023model, outcomes were considered missing after two consecutive missed tests. We constructed a design matrix $\mathbf{X} \in {\mathbb{R}}^{n \times p}$ with $n=587, p=6$, an outcome vector $\mathbf{y} \in {\mathbb{R}}^{n}$, and a treatment assignment vector (296 controls, 291 treated). Our target is the treatment-specific mean $\bar{\tau}$. \begin{table} \begin{center} \caption{Treatment effect estimates of the proportion of positive urine laboratory results among all tests, using data from NIDA-CTN-0030. The four approaches considered are (a) $\widehat{\tau}_{\mathsf{unadj}}$, (b) $\widehat{\tau}_{\mathsf{adj},2}$, (c) $\widehat{\tau}_{\mathsf{adj},2}^{\dag}$ and (d) $\widehat{\tau}_{\mathsf{adj},3}$. All values are multiplied by 10. The superscript “c" refers to the conservative standard error estimates.} \begin{tabular}{cccccc} \Xhline{3\arrayrulewidth} \textbf{Method} & \textbf{Estimate} & \textbf{SE} & \textbf{$95\%$-CI} & \textbf{SE$^{c}$} & \textbf{$95\%$-CI$^{c}$} \\ \Xhline{3\arrayrulewidth} $\widehat{\tau}_{\mathsf{unadj}}$ & $1.1836$ & $0.0372$ & $(1.1108, 1.2564)$ & NA & NA \\ \hline $\widehat{\tau}_{\mathsf{adj},2}$ & $1.1837$ & $0.0327$ & (1.1196, 1.2489) & 0.0328 & (1.1195, 1.2479) \\ \hline $\widehat{\tau}_{\mathsf{adj},2}^{\dag}$ & $1.1755$ & $0.0321$ & (1.1126, 1.2383) & 0.0318 & (1.1131, 1.2379) \\ \hline $\widehat{\tau}_{\mathsf{adj},3}$ & $1.1838$ & $0.0327$ & (1.1196, 1.2489) & 0.0328 & (1.1196, 1.2480) \\ \Xhline{3\arrayrulewidth} \end{tabular} \end{center} \end{table} In Table (ref), we present the point estimate, along with the non-conservative and conservative estimates of standard deviation, and the $95\%$ confidence interval for each of the four different estimators, excluding $\widehat{\tau}_{\mathsf{unadj}}$. Compared to other adjusted estimators, the point estimate of $\widehat{\tau}_{\mathsf{adj},2}$ and the bias-free $\widehat{\tau}_{\mathsf{adj},3}$ that we proposed are closer to the unadjusted estimator $\widehat{\tau}_{\mathsf{unadj}}$. Furthermore, their standard errors are also lower than $\widehat{\tau}_{\mathsf{unadj}}$. \section{Concluding Remarks} In this paper, we demonstrate that HOIF-motivated estimators are natural candidate treatment effect estimators adjusting for high-dimensional baseline covariates in RCT. We prove statistical properties of the HOIF-motivated estimator under CRE or Bernoulli sampling and the design-based framework, which, to the best of our knowledge, is new in the HOIF literature. More importantly, we show that HOIF-motivated estimators place several recently proposed treatment effect estimators adjusting for high-dimensional baseline covariates in a unified framework. This result further consolidates the role of HOIF as a useful template to construct “good” estimators from first principles, without resorting to clever tricks. Finally, there are two main future directions that worth pursuing. First, it is of theoretical interest to use the design-based Riesz representation theory of harshaw2022design to justify that the HOIF-motivated estimators are indeed HOIF in the design-based framework. Second, it is of practical interest to extend our estimators to designs beyond CRE. \section*{Acknowledgments} Section (ref) of this paper is partly motivated by a conversation with Oliver Dukes, who kindly suggested to LL that a somewhat pedagogical/elementary review of HOIF could be useful to practitioners. The authors would also like to express their sincere gratitude to Yujia Gu and \href{https://maweiruc.github.io/}{Wei Ma} for enlightening discussions and Muluneh Alene Addis and Kelly Van Lancker for spotting an issue in the variance estimator in the first version of the R package \texttt{HOIFCar}. \putbib[Master]