Integrating Heterogeneous Information in Randomized Experiments: A Unified Calibration Framework

spacing{1.5} \begin{center} {\bf Calibration for Treatment Effect Estimation in Randomized Controlled Trials: a Unified Approach for Covariates Adjustment and Information Borrowing} \end{center}

bibunit\begin{abstract} In modern randomized experiments, large-scale data collection increasingly yields rich baseline covariates and auxiliary information from multiple sources. Such information offers opportunities for more precise treatment effect estimation, but it also raises the challenge of integrating heterogeneous information coherently without compromising validity. Covariate-adaptive randomization (CAR) is widely used to improve covariate balance at the design stage, but it typically balances only a small set of covariates used to form strata, making covariate adjustment at the analysis stage essential for more efficient estimation of treatment effects. Beyond standard covariate adjustment, it is often desirable to incorporate auxiliary information, including cross-stratum information, predictions from various machine learning models, and external data from historical trials or real-world sources. While this auxiliary information is widely available, existing covariate adjustment methods under CAR primarily exploit within-stratum covariates and do not provide a coherent mechanism for integrating it. We propose a unified calibration framework that integrates such information through an information proxy vector and calibration weights defined by a convex optimization problem. The resulting estimator recovers many recent covariate adjustment procedures as special cases while providing a systematic mechanism for both internal and external information borrowing within a single framework. We establish large-sample validity and a no-harm efficiency guarantee, showing that incorporating additional information sources cannot increase asymptotic variance, and we extend the theory to settings in which both the number of strata and the number of information sources grow with the sample size. Simulation studies and an empirical analysis of a field experiment on savings behavior in Uganda and Malawi demonstrate the strong finite-sample performance and practical utility of our method. \end{abstract} {\it Keywords:} Calibration weights; Covariate-adaptive randomization; Covariate adjustment; Information integration; Machine learning \doublespacing \section{Introduction} In modern randomized experiments, ensuring balance in baseline covariates across treatment groups is critical for reducing bias and improving the credibility of the trial results. Covariate-adaptive randomization (CAR) methods, such as stratified biased coin randomization (efron1971Forcing), stratified block randomization (zelen1974Randomization), and minimization (pocock1975Sequential,taves1974Minimization), are widely used to achieve this balance during the design stage. In practice, however, CAR is typically implemented using only a small set of covariates to form strata, so balance is not directly enforced for many other baseline covariates. Moreover, some important pre-treatment covariates may only be observed after randomization, for instance when they are collected together with the outcome variable (bai2024Covariate). Consequently, covariate adjustment at the statistical analysis stage plays a crucial complementary role. By incorporating baseline covariates into the statistical analysis, covariate adjustment methods correct for residual imbalances, thereby improving the precision and efficiency of treatment effect estimates and strengthening the validity of the trial. The adjustment of baseline covariates in the statistical analysis stage has been studied for a long time (see, e.g., tsiatis2008Covariate,zhang2008Improving,lin2013Agnostic and the references therein). Under the CAR design, ma2022Regression,ye2022Inference,gu2023RegressionBased adjusted for additional covariates using linear regression, while liu2023Lassoadjusted applied Lasso (tibshirani1996Regression). These approaches resulted in average treatment effect (ATE) estimators that are more efficient than the naive difference-in-means estimator and the ordinary least squares estimator (regressing outcomes on strata indicators) proposed by bugni2018Inference and bugni2019Inference. This holds true even in the absence of strong evidence for a linear relationship between covariates and outcomes. To further enhance efficiency, recent work by rafi2023Efficient, tu2024Unified and bannick2025General introduced unified frameworks for nonlinear covariate adjustment. Relying on the augmented inverse probability weighting (AIPW, robins1994Estimation), these frameworks incorporate machine learning techniques, such as random forests and deep neural networks, to adjust for covariates. These machine learning methods can be effective in capturing complex nonlinear relationships between covariates and outcomes. However, much of the existing covariate adjustment literature can be interpreted as focusing on a particular form of internal information borrowing, namely using baseline covariates from the current trial, typically within each stratum, to improve efficiency. This focus leaves relatively little room for other practically important forms of information borrowing, and it limits the extent to which standard adjustment methods can integrate heterogeneous information in a systematic way. Internally, efficiency can often be improved by borrowing information across strata when the outcome--covariate relationship is stable, and by combining predictions from multiple machine learning methods when no single learner is uniformly reliable. Externally, there is increasing interest in leveraging historical trials and real-world data to support analyses of concurrent trials, especially when sample sizes are constrained by cost, ethical considerations, or recruitment challenges FDA2019RareDiseases,gu2024incorporatingexternaldataanalyzing. Existing AIPW-based nonlinear adjustment frameworks (tu2024Unified,rafi2023Efficient,bannick2025General), however, are not designed to integrate these heterogeneous information sources, as they typically rely on a single nuisance estimate and lack a systematic mechanism for combining multiple internal predictors or external information sources. To fill this gap, we propose a unified calibration framework for integrating heterogeneous information that accommodates both internal and external borrowing, thereby enabling estimators that jointly exploit covariates, cross-stratum information, and auxiliary data sources. A central methodological feature of our approach is the use of calibration weights under CAR designs. Although calibration has been studied in survey sampling (deville1992Calibration,kwon2025Debiased), missing data problems (qin2007EmpiricalLikelihoodBased,tan2014Secondorder), and observational studies (chan2016Globally), its analysis under CAR designs raises distinct theoretical issues. Unlike the aforementioned settings where samples are typically independent and identically distributed (i.i.d.), CAR designs induce complex dependence structures among treatment assignments within strata. We address this dependence through a conditional asymptotic argument. Specifically, we condition on the realized stratum indicators and treatment assignments, treat them as fixed, and establish large-sample results using conditional laws of large numbers and conditional central limit theorems. This approach provides a tractable framework for inference under CAR and yields proof techniques that extend to regimes with a growing number of strata and an increasing number of information sources, which may be of independent interest. We summarize the main contributions as follows. \begin{itemize} • A unified calibration framework. We introduce a calibration-based framework for estimation and inference under CAR. This framework is unified in three aspects. First, it provides a common formulation that recovers many recent covariate adjustment procedures as special cases (e.g., bugni2019Inference,cohen2024Noharm,ma2022Regression,tu2024Unified,bannick2025General,ye2022Inference,ye2023Better,liu2023Lassoadjusted). Second, it places internal and external information borrowing within a single architecture, yielding a systematic approach to integrating heterogeneous information sources. Third, it applies broadly across CAR schemes satisfying Assumption (ref), so that the resulting inference procedure is not tied to a particular randomization method. • Flexible and robust information borrowing. We develop practical constructions of the information proxy vector that accommodate a wide range of internal and external information sources. Internally, the framework can borrow across strata and aggregate heterogeneous machine learning predictions. Externally, it can incorporate information from historical trials and real-world data. Importantly, our framework is model-agnostic regarding information sources. We refer to this property as robustness, meaning that the validity of our statistical inference holds even if the utilized information is biased or generated by inaccurate models. • General inference theory under CAR. We provide a rigorous theoretical foundation, proving that our estimator is asymptotically normal with a consistently estimable variance. We establish a no-harm efficiency guarantee, ensuring that utilizing additional information sources improves, or at worst maintains, estimation efficiency. Furthermore, we develop proof techniques tailored to CAR-induced dependence, distinct from existing arguments (e.g., bugni2018Inference,bugni2019Inference,ma2022Regression,liu2023Lassoadjusted). These techniques accommodate a growing number of strata and an increasing number of information sources, offering tools that may be of independent interest in related problems. \end{itemize} The paper is organized as follows. Figure (ref) provides an overview. Section (ref) introduces the setting and our unified calibration framework. Section (ref) discusses practical strategies for constructing the information proxy vector $\bs{\xi}_{n}$, including approaches that borrow auxiliary information from both internal and external sources. Section (ref) establishes the large-sample properties of the proposed estimator, including asymptotic normality, consistent variance estimation, and efficiency comparisons. Section (ref) presents theoretical extensions, including settings with diverging numbers of strata and a growing dimension of $\bs{\xi}_{n}$, as well as general discrepancy measures. Finally, Sections (ref) and (ref) assess finite-sample performance via simulation and illustrate the method using experimental data from dupas2018Bankinga. All proofs are provided in the Appendix (Supplementary Material, ma2026Integrating). \begin{figure}[!tbh] \resizebox{\textwidth}{!}{ \tikzset{ base/.style = { rectangle, rounded corners, draw=black, text centered, font=, blur shadow={shadow blur steps=5} }, topnode/.style = { base, fill=blue!10, draw=blue!80!black, line width=1.5pt, minimum width=10cm, minimum height=1.2cm, font= }, midnode/.style = { base, fill=orange!10, draw=orange!80!black, line width=1.2pt, minimum width=5.5cm, minimum height=2cm, text width=5.2cm }, subnode/.style = { base, fill=gray!5, draw=gray!60!black, dashed, minimum width=3.5cm, minimum height=1.2cm, text width=3.2cm, font= }, connector/.style = { ->, >=stealth, line width=1.2pt, color=black!70 }, subconnector/.style = { -, dashed, line width=0.8pt, color=gray!70 } } \begin{tikzpicture}[node distance=2cm and 2.5cm] \node (sec2) [topnode] {Sec 2: A Unified Calibration Framework}; \node (sec4) [midnode, below=2.5cm of sec2] {Sec 4: Asymptotics\\ Establish normality & efficiency comparison}; \node (sec3) [midnode, left=of sec4] {\textbf{Sec 3: Construction of $\bs{\xi}_{n}$}\\ Leverage internal & external information}; \node (sec5) [midnode, right=of sec4] {\textbf{Sec 5: Extensions}\\ Two extensions of the asymptotic results}; \node (sub3a) [subnode, below=1cm of sec3, xshift=-2cm] {\textbf{Internal:} \\ cross-stratum information & various machine learning & cross-fitting}; \node (sub3b) [subnode, below=1cm of sec3, xshift=2cm] {\textbf{External:} \\ historical trials & real-world data}; \node (sub4a) [subnode, below=1cm of sec4, xshift=-2cm] {asymptotic normality & valid inference}; \node (sub4b) [subnode, below=1cm of sec4, xshift=2cm] {guaranteed \\ efficiency gain}; \node (sub5a) [subnode, below=1cm of sec5, xshift=-2cm] {diverging number of strata and information sources}; \node (sub5b) [subnode, below=1cm of sec5, xshift=2cm] {general discrepancy measure $D(v)$ & second-order bias}; \draw [connector] (sec2.south) -- +(0,-0.8) -| (sec3.north); \draw [connector] (sec2.south) -- (sec4.north); \draw [connector] (sec2.south) -- +(0,-0.8) -| (sec5.north); \draw [subconnector] (sec3.south) -- ++(0,-0.5) -| (sub3a.north); \draw [subconnector] (sec3.south) -- ++(0,-0.5) -| (sub3b.north); \draw [subconnector] (sec4.south) -- ++(0,-0.5) -| (sub4a.north); \draw [subconnector] (sec4.south) -- ++(0,-0.5) -| (sub4b.north); \draw [subconnector] (sec5.south) -- ++(0,-0.5) -| (sub5a.north); \draw [subconnector] (sec5.south) -- ++(0,-0.5) -| (sub5b.north); \end{tikzpicture} } \caption{Roadmap of the proposed unified calibration framework and theoretical analysis. } \end{figure} \textit{Notation.} We maintain the following notation conventions throughout the paper. For any column vector $\bs x=\left(x_{1},x_{2},\ldots,x_{d}\right)^{\top}\in\mathbb{R}^{d}$, where $\mathbb{R}^{d}$ is the $d$-dimensional Euclidean space, $\left\Vert \bs x\right\Vert =\left(\bs x^{\top}\bs x\right)^{1/2}$ denotes its Euclidean norm. For any matrix $A=\left(a_{ij}\right)_{n\times m}$, $\left\Vert A\right\Vert $ denotes its maximum singular value, i.e., the operator norm, $\left\Vert A\right\Vert _{F}=\sqrt{\mathrm{tr}\left(AA^{\top}\right)}$ denotes it Frobenius norm, and $A^{+}$ denotes its Moore-Penrose inverse. For two positive non-random sequence $a_{n},b_{n}$ and random vector sequence $X_{n}$, $X_{n}=o_{P}(a_{n})$ means $P\left(\left\Vert X_{n}\right\Vert >a_{n}\epsilon\right)\to0$ as $n\to\infty$ for any $\epsilon>0$ and $X_{n}=O_{P}(a_{n})$ means for any $\epsilon>0$, there exists a constant $M>0$ such that $\limsup_{n\to\infty}P\left(\left\Vert X_{n}\right\Vert \geq a_{n}M\right)<\epsilon$. The notation $\1(\cdot)$ denotes the indicator function, which takes the value 1 if the condition inside the parentheses is true and 0 otherwise. \section{A unified calibration framework} \subsection{Preliminaries} Let $A_{i}\in\{0,1\}$ ($i=1,\ldots,n$) denote the treatment assignment indicator, where $A_{i}=1$ indicates that the $i$-th unit is assigned to the treatment group and $A_{i}=0$ otherwise. We assume the assignments $\{A_{i}\}_{i=1}^{n}$ are generated via a CAR design satisfying Assumption (ref). Consequently, the variables $A_{i}$ are typically not i.i.d. Adopting the potential outcomes framework (imbens2015Causal), we define $Y_{i}(1)$ and $Y_{i}(0)$ as the potential outcomes under treatment and control, respectively. The observed outcome $Y_{i}$ is determined by $Y_{i}=A_{i}Y_{i}(1)+(1-A_{i})Y_{i}(0)$. The experimental design stratifies units into $K$ strata, with $B_{i}\in\{1,\ldots,K\}$ indicating the stratum of unit $i$. For notational convenience, let $[k]=\{i:B_{i}=k\}$ denote the set of indices for units belonging to stratum $k$. To ensure non-empty strata, we assume positive assignment probabilities: $p_{[k]}=P(B_{i}=k)>0$ for all $k\in\{1,\ldots,K\}$ and $i\in\{1,\ldots,n\}$. The target treatment allocation proportion stratum $k$ is $\pi_{[k]}=P(A_{i}=1\mid B_{i}=k)\in(0,1)$. Each unit has a $p$-dimensional baseline covariate vector $\bs X_{i}=(X_{i1},\ldots,X_{ip})^{\top}\in\mathcal{X}\subset\mathbb{R}^{p}$, which may be either low- or high-dimensional. We assume that the covariate vector $\bs X_{i}$ contains the stratum indicator $B_{i}$, but to highlight the stratum indicator $B_{i}$, we sometimes use the notation $(\bs X_{i},B_{i})$. Let subscripts 1 and 0 denote treatment and control groups, respectively. The treatment group contains $n_{1}=\sum_{i=1}^{n}A_{i}$ units and the control group contains $n_{0}=\sum_{i=1}^{n}(1-A_{i})$ units. For stratum-specific quantities, we use subscript $[k]$: let $n_{[k]}=\sum_{i\in[k]}1$ denote the stratum size, with $n_{1[k]}=\sum_{i\in[k]}A_{i}$ and $n_{0[k]}=\sum_{i\in[k]}(1-A_{i})$ representing treated and control units in stratum $k$, respectively. The stratum proportion and treatment allocation proportion are defined as $p_{n[k]}=n_{[k]}/n$ and $\pi_{n[k]}=n_{1[k]}/n_{[k]}$, correspondingly. Our parameter of interest is the population average treatment effect (ATE): \[ \tau=\mathbb{E}[Y_{i}(1)-Y_{i}(0)]. \] Under CAR, the ATE parameter $\tau$ can be consistently estimated by aggregating the treatment effect estimates from each stratum. A simple estimator for this is the stratified difference-in-means estimator (bugni2019Inference,ma2022Regression): \[ \widehat{\tau}_{\mathrm{sdim}}:=\sum_{k=1}^{K}p_{n[k]}\left(\overline{Y}_{1[k]}-\overline{Y}_{0[k]}\right), \] where $\overline{Y}_{a[k]}:=\frac{1}{n_{a[k]}}\sum_{i\in[k]}\1(A_{i}=a)Y_{i}$ for $a\in\{0,1\}$ and $k=1,\ldots,K$. To improve estimation efficiency, recent literature has proposed nonlinear covariate adjustment methods (tu2024Unified,rafi2023Efficient,bannick2025General). These methods typically use the augmented inverse probability weighting (AIPW, robins1994Estimation) estimator and rely on estimating the conditional mean function $h_{a[k]}^{*}(\bs X_{i})=\mathbb{E}\left[Y_{i}(a)\mid\bs X_{i},B_{i}=k\right]$. However, AIPW-based adjustment methods are limited in that they cannot combine different estimates of $h_{a[k]}^{*}(\bs X_{i})$. For instance, these estimates may leverage internal information via cross-stratum borrowing and heterogeneous machine learning predictions, or incorporate external information from historical trials, real-world data, and functional forms suggested by domain experts (see Section (ref) for a detailed discussion). In the following, we will introduce a unified adjustment framework that can combine these diverse estimates of $h_{a[k]}^{*}(\bs X_{i})$, leading to more efficient and robust ATE estimators. \subsection{Construction of the calibration estimator} Let $D(v):\mathbb{R}\to\mathbb{R}$ be a twice continuously differentiable and strictly convex function that measures the discrepancy from $v$ to 1, e.g., $D(v)=(v-1)^{2}/2$ and $D(v)=v-\log v$. Suppose $\bs{\xi}_{n}:\mathcal{X}\to\mathbb{R}^{d}$ ($n\geq1$) is a sequence of $\mathbb{R}^{d}$-valued (possibly) random functions of $\bs X_{i}$, which we refer to as the \textit{information proxy vector}. Typically, the elements of $\bs{\xi}_{n}$ may be different estimates of the conditional mean function $h_{a[k]}^{*}(\bs X_{i})$. Based on $\bs{\xi}_{n}$, we propose the calibration estimator: \begin{equation} \widehat{\tau}_{\mathrm{cal}}:=\widehat{\tau}_{\mathrm{sdim}}+\frac{1}{n}\sum_{i=1}^{n}\widehat{w}_{i}r_{i}, \end{equation} where $r_{i}:=\sum_{k=1}^{K}\bigl\{\frac{A_{i}}{\pi_{n[k]}}(Y_{i}-\overline{Y}_{1[k]})-\frac{1-A_{i}}{1-\pi_{n[k]}}(Y_{i}-\overline{Y}_{0[k]})\bigr\}\1(B_{i}=k)$ and the calibration weights $\widehat{w}_{i}$'s ($i=1,\ldots,n$) solve the calibration problem: \begin{equation} \begin{cases} (\widehat{w}_{1},\ldots,\widehat{w}_{n})=\min_{w_{i},1\leq i\leq n}\sum_{i=1}^{n}D(w_{i})\ \text{subject to }\\ \frac{1}{n}\sum_{i=1}^{n}w_{i}\left\{ A_{i}-\pi_{n[k]}\right\} \1(B_{i}=k)\left\{ \bs{\xi}_{n}(\bs X_{i})-\overline{\bs{\xi}}_{n[k]}\right\} =0, & \forall k=1,\ldots,K. \end{cases} \end{equation} In ((ref)), $\overline{\bs{\xi}}_{n[k]}:=\frac{1}{n_{[k]}}\sum_{i\in[k]}\bs{\xi}_{n}(\bs X_{i})$ is defined as the stratum-specific sample mean for $\bs{\xi}_{n}(\bs X_{i})$. The calibration problem ((ref)) is a convex optimization problem involving $dK$ linear constraints. This can be efficiently solved using standard convex optimization software via its dual formulation (see boyd2004Convex). The calibration estimator ((ref)) consists of two components: (i) the stratified difference-in-means estimator $\widehat{\tau}_{\mathrm{sdim}}$ and (ii) a correction term constructed from weighted residuals, where the weights are determined by the calibration problem ((ref)). Analogous to linear regression analysis, the residuals $r_{i}=\sum_{k=1}^{K}\bigl\{\frac{A_{i}}{\pi_{n[k]}}(Y_{i}-\overline{Y}_{1[k]})-\frac{1-A_{i}}{1-\pi_{n[k]}}(Y_{i}-\overline{Y}_{0[k]})\bigr\}\1(B_{i}=k)$, $i=1,\ldots,n$, represent the part of $\sum_{k=1}^{K}\bigl\{\frac{A_{i}}{\pi_{n[k]}}Y_{i}-\frac{1-A_{i}}{1-\pi_{n[k]}}Y_{i}\bigr\}\1(B_{i}=k)$ that remains unexplained by the stratum mean. Geometrically, the residuals $r_{i}$ arise from projecting $(\sum_{k=1}^{K}\bigl\{\frac{A_{i}}{\pi_{n[k]}}Y_{i}-\frac{1-A_{i}}{1-\pi_{n[k]}}Y_{i}\bigr\}\1(B_{i}=k):1\leq i\leq n)$ onto the orthogonal complement of the space spanned by $(1,...,1)$, which inherently implies $\sum_{i=1}^{n}r_{i}=0$. The calibration problem ((ref)) enforces the balance of $\bs{\xi}_{n}(\bs X_{i})-\overline{\bs{\xi}}_{n[k]}$ across different treatment groups within strata when the observed samples are weighted by $\widehat{w}_{i}$, $i=1,\ldots,n$, thereby naturally incorporating the information contained in $\bs{\xi}_{n}$ into the calibration weights $\widehat{w}_{i}$'s. This incorporation of information allows the calibration weights $\widehat{w}_{i}$ to further explain the variability in $\sum_{k=1}^{K}\bigl\{\frac{A_{i}}{\pi_{n[k]}}Y_{i}-\frac{1-A_{i}}{1-\pi_{n[k]}}Y_{i}\bigr\}\1(B_{i}=k)$. Consequently, the weighted average of the residuals, $\frac{1}{n}\sum_{i=1}^{n}\widehat{w}_{i}r_{i}$, acts as a correction term representing the part explained by the calibration weights. If $\widehat{w}_{i}=1$ for all $i=1,\ldots,n$, which occurs, for instance, when $\bs{\xi}_{n}(\bs X)-\overline{\bs{\xi}}_{n[k]}$ is already balanced by the randomization procedure, then $\widehat{\tau}_{\mathrm{cal}}$ reduces to $\widehat{\tau}_{\mathrm{sdim}}$. Our calibration estimator ((ref)) provides a unified information integration framework. With an appropriate choice of the information proxy vector $\boldsymbol{\xi}_{n}$, it recovers many existing covariate adjustment methods, which can be viewed as a special case of internal information borrowing (e.g., liu2023Lassoadjusted,cohen2024Noharm,tu2024Unified,bannick2025General,ye2022Inference,ye2023Better,gu2024incorporatingexternaldataanalyzing,bugni2019Inference,ma2022Regression). For instance, when $\boldsymbol{\xi}_{n}$ is specified as a Lasso-based estimate, ((ref)) becomes analogous to the Lasso-adjusted estimator of liu2023Lassoadjusted. Beyond internal adjustment, ((ref)) also accommodates external information proxies, in parallel to external-borrowing approaches such as gu2024incorporatingexternaldataanalyzing, and thus enables new estimators that jointly leverage internal and external information. A key property of our calibration estimator ((ref)) is its invariance under affine transformations of $\bs{\xi}_{n}(\bs X)$. Specifically, replacing $\bs{\xi}_{n}(\bs X)$ with $\mathbf{Q}\bs{\xi}_{n}(\bs X)+\bs q$ in the optimization problem ((ref)) results in the same calibration estimator $\widehat{\tau}_{\mathrm{cal}}$, where $\mathbf{Q}\in\mathbb{R}^{d\times d}$ is any invertible $d\times d$ matrix and $\bs q\in\mathbb{R}^{d}$ is any $d$-dimensional vector. If $\bs{\xi}_{n}(\bs X)$ serves as an estimate of the conditional mean function $h_{a[k]}^{*}(\bs X_{i})$, this property allows for misspecification up to an affine transformation, offering a flexibility that is absent in standard AIPW-based covariate adjustment methods. \section{Strategies for constructing the information proxy $\protect\bs{\xi}_{n}$} The efficacy and flexibility of the proposed unified framework relies on the construction of $\bs{\xi}_{n}$, which serves as a proxy for auxiliary information. Accordingly, we categorize the strategies for determining $\bs{\xi}_{n}$ into internal and external information borrowing. Sections (ref)--(ref) focus on internal strategies, encompassing cross-stratum information borrowing, the aggregation of different machine learning predictions, and cross-fitting techniques, whereas Section (ref) discusses the incorporation of external information such as historical trials and real-world data. \subsection{Cross-stratum information borrowing} Suppose $\widehat{h}_{a[k]}(\cdot)$, where $a\in\{0,1\}$ and $1\leq k\leq K$, are the estimators for $h_{a[k]}^{*}(\cdot)$. According to Theorem (ref) (to be established in Section (ref)), if $\widehat{h}_{a[k]}(\cdot)$ is consistent to $h_{a[k]}^{*}(\cdot)$, then taking $\bs{\xi}_{n}(\bs X)=\left(\sum_{k=1}^{K}\widehat{h}_{1[k]}(\bs X)\1(B=k),\sum_{k=1}^{K}\widehat{h}_{0[k]}(\bs X)\1(B=k)\right)^{\top}$ in ((ref)) will lead to a semiparametrically efficient estimator for $\tau$. However, it is important to note that this approach ensures each stratum performs covariate adjustment using only the information available within that stratum itself.\textcolor{red}{ }However, borrowing information across strata is often advantageous, particularly when the relationship between covariates $\bs X$ and potential outcomes $(Y(1),Y(0))$ is stable across different strata. Our framework accommodates this naturally by taking $\bs{\xi}_{n}(\bs X)=\left(\widehat{h}_{1[k]}(\bs X),\widehat{h}_{0[k]}(\bs X):1\leq k\le K\right)^{\top}$ in ((ref)). This specification allows each stratum to use information from all strata, leading to a more efficient calibration estimator. \subsection{Integration of heterogeneous machine learning predictions} In practice, a variety of machine learning methods can be used to estimate the conditional mean function $h_{a[k]}^{*}(\bs X)$, such as deep neural networks (lecun2015deep,jiao2023Deep,farrell2021Deepa), random forests (breiman2001random,wager2018estimation), Lasso (tibshirani1996Regression), and others. Our framework provides a natural approach to combine these machine learning estimates. Suppose $\widehat{h}_{a[k]}^{\text{rf}}(\bs X)$ and $\widehat{h}_{a[k]}^{\text{nn}}(\bs X)$ are the estimates of $h_{a[k]}^{*}(\bs X)$ based on random forests and deep neural networks, respectively. We can then let $\bs{\xi}_{n}(\bs X)=\left(\widehat{h}_{a[k]}^{\text{rf}}(\bs X),\widehat{h}_{a[k]}^{\text{nn}}(\bs X):a\in\{0,1\},1\leq k\le K\right)^{\top}$ in ((ref)) to derive a calibration estimator for $\tau$. According to Theorem (ref) (to be established in Section (ref)), this estimator will be more efficient than those based solely on random forests or neural networks when a single machine learning method fails to fully capture the conditional mean function $h_{a[k]}^{*}(\bs X)$. This improvement enhances the efficiency of covariate adjustment methods. \subsection{Implementation via cross-fitting and sample splitting} When estimating $\bs{\xi}_{n}$ via machine learning methods, it is often desirable to rely on an independent sample in order to mitigate overfitting and to make Assumption (ref) more plausible. Such independence can be achieved through sample-splitting. However, sample-splitting typically reduces efficiency. To recover efficiency, we adopt the cross-fitting technique (see chernozhukov2018Double,tu2024Unified,bannick2025General,rafi2023Efficient). Specifically, we follow the sample-splitting procedure outlined in rafi2023Efficient, which ensures that each fold contains data from every stratum and treatment arm. Suppose the units $\{1,\ldots,n\}$ are partitioned into two equal folds, $I_{0}$ and $I_{1}$. The cross-fitted calibration estimator is then defined as $\widehat{\tau}_{\mathrm{cal}}^{\mathrm{CF}}:=\frac{1}{2}\widehat{\tau}_{\mathrm{cal}}^{(0)}+\frac{1}{2}\widehat{\tau}_{\mathrm{cal}}^{(1)}$, where \begin{align*} \widehat{\tau}_{\mathrm{cal}}^{(\iota)} & =\sum_{k=1}^{K}p_{n[k]}^{(\iota)}\left(\overline{Y}_{1[k]}^{(\iota)}-\overline{Y}_{0[k]}^{(\iota)}\right)\\ & \quad+\frac{1}{\left|I_{\iota}\right|}\sum_{i\in I_{\iota}}\widehat{w}_{i}\sum_{k=1}^{K}\left\{ \frac{A_{i}}{\pi_{n[k]}^{(\iota)}}\left(Y_{i}-\overline{Y}_{1[k]}^{(\iota)}\right)-\frac{1-A_{i}}{1-\pi_{n[k]}^{(\iota)}}\left(Y_{i}-\overline{Y}_{0[k]}^{(\iota)}\right)\right\} \1(B_{i}=k) \end{align*} with $\left|I_{\iota}\right|=\sum_{i\in I_{\iota}}1$, $p_{n[k]}^{(\iota)}=\sum_{i\in I_{\iota}}\1(B_{i}=k)/\left|I_{\iota}\right|$, $\overline{Y}_{a[k]}^{(\iota)}=\sum_{i\in I_{\iota}}\1(A_{i}=a,B_{i}=k)Y_{i}/\sum_{i\in I_{\iota}}\1(A_{i}=a,B_{i}=k)$, $\pi_{n[k]}^{(\iota)}=\sum_{i\in I_{\iota}}\1(A_{i}=a,B_{i}=k)/\sum_{i\in I_{\iota}}\1(B_{i}=k)$ for $\iota\in\{0,1\}$. The weights $\widehat{w}_{i}$'s ($i\in I_{\iota}$) are obtained by solving the following calibration problem: \[ \begin{cases} (\widehat{w}_{i})_{i\in I_{\iota}}=\min_{w_{i}:i\in I_{\iota}}\sum_{i\in I_{\iota}}D(w_{i})\ \text{subject to }\\ \sum_{i\in I_{\iota}}w_{i}\left\{ A_{i}-\pi_{n[k]}^{(\iota)}\right\} \1(B_{i}=k)\left\{ \bs{\xi}_{n}^{(1-\iota)}(\bs X_{i})-\overline{\bs{\xi}}_{n[k]}^{(1-\iota)}\right\} =0, & \forall k=1,\ldots,K, \end{cases} \] where $\bs{\xi}_{n}^{(1-\iota)}$ is measurable with respect to $\{(Y_{i},\bs X_{i},A_{i}):i\in I_{\iota}\}$ and $\overline{\bs{\xi}}_{n[k]}^{(1-\iota)}:=\sum_{i\in I_{\iota}}\bs{\xi}_{n}(\bs X_{i})\1(B_{i}=k)/\sum_{i\in I_{\iota}}\1(B_{i}=k)$ denotes the stratum-specific sample mean for $\bs{\xi}_{n}^{(1-\iota)}(\bs X)$. The asymptotic properties of $\widehat{\tau}_{\mathrm{cal}}^{\mathrm{CF}}$ can be established by combining the proof of Theorem (ref) with the arguments from the proof of rafi2023Efficient and we omit the details here. We will evaluate the finite-sample performance of this cross-fitted estimator $\widehat{\tau}_{\mathrm{cal}}^{\mathrm{CF}}$ in our simulation studies. \subsection{Leveraging historical and real-world data} Beyond internal strategies, constructing $\bs{\xi}_{n}$ using external information offers a powerful way to enhance efficiency. In this subsection, we focus on two primary sources: historical clinical trials and real-world data. While these sources provide valuable information, they often differ distributionally from the current trial. A key distinction of our framework is its ability to leverage such heterogeneous external data without requiring restrictive similarity assumptions. First, historical clinical trials often provide readily available data from settings similar to the current study. Since treatments often vary between trials, existing literature has predominantly focused on leveraging historical control data to enhance the precision of estimates in the current trial (pocock1975Sequential,callegaro2023Historical). These methods usually rely on the assumption that the historical and current control arms are comparable. In contrast, a distinct advantage of our unified framework is its assumption-lean nature: we impose no assumptions on the validity or direct transferability of the borrowed information. Consequently, we are not restricted to borrowing solely from historical control arms; we can flexibly incorporate information from historical studies involving treatments that may differ from those in the current trial. Second, beyond historical clinical trials, real-world data represents another vast information source. Real-world data encompasses data generated from routine healthcare delivery and patient monitoring, including electronic health records (EHRs), disease registries, and increasingly, data from wearable devices. Unlike tightly controlled clinical trials, observational studies based on real-world data typically involve much larger and more diverse patient populations, which can significantly strengthen the statistical power of the analysis. Crucially, both historical trial data and real-world data may exhibit a “covariate shift”, i.e., the distribution of covariates $\bs X$ in these external sources often differs from that of the current trial. However, a common phenomenon is that while the marginal distribution of $\bs X$ changes, the conditional distribution of the potential outcomes $(Y(1),Y(0))$ given the covariates $\bs X$ remains comparable. This stability is important because the optimal choice of $\bs{\xi}_{n}(\bs X)$, which ensures that the resulting estimator $\widehat{\tau}_{\mathrm{cal}}$ is semiparametrically efficient, depends solely on the conditional distribution of $(Y(1),Y(0))$ given $\bs X$ (see Theorem (ref)). In practice, we can leverage this by estimating the conditional mean functions $h_{1[k]}^{*}(\bs X)$ and $h_{0[k]}^{*}(\bs X)$ ($1\leq k\leq K$) from these external datasets and setting $\bs{\xi}_{n}(\bs X)$ as the vector of these estimates. Our Theorem (ref) guarantees that incorporating these estimates from the external datasets will result in a calibration estimator with no greater asymptotic variance. Consequently, our framework provides a robust strategy for improving estimation efficiency without the risk of “negative transfer”. Unlike existing methods that rely on Bayesian frameworks (hobbs2011Hierarchical,ibrahim2015Power) or Trans-Lasso (gu2024incorporatingexternaldataanalyzing), where the transfer of information depends on the prior distributions and/or model assumptions, our approach is entirely model-free. Notably, it imposes no similarity constraints between the external data and the target data-generating process, thereby enhancing flexibility and robustness against real-world data complexities. \section{Asymptotic properties} We make the following assumptions. \begin{assumption} $\bs W_{i}=(Y_{i}(1),Y_{i}(0),\bs X_{i}^{\top},B_{i})^{\top},$ $i=1,\ldots,n$, are i.i.d. samples from the population distribution of $\bs W=(Y(1),Y(0),\bs X^{\top},B)^{\top}$, and we denote $\bs W^{(n)}=\{\bs W_{1},\ldots,\bs W_{n}\}$. Besides, $\sup_{1\leq k\leq K}\left[\left|Y_{i}(a)\right|^{2+\epsilon}\mid B_{i}=k\right]<\infty$ for both $a=0,1$, where $0<\epsilon\leq1$ is a constant. \end{assumption} \begin{assumption} [Treatment assignment]Let $A^{(n)}=\{A_{1},\ldots,A_{n}\}$ and $B^{(n)}=\{B_{1},\ldots,B_{n}\}$. The treatment assignment mechanism satisfies the following conditions: \begin{enumerate}[label=(\arabic*)] • conditional on the stratum indicators $B^{(n)}$, the treatment assignments $A^{(n)}$ are independent of $\bs W^{(n)}$, i.e., $\bs W^{(n)}\perp A^{(n)}\mid B^{(n)}$; • $\sup_{1\leq k\leq K}\left|\pi_{n[k]}-\pi_{[k]}\right|\overset{p}{\to}0$ as $n\to\infty$, where $\pi_{[k]}$, $k=1,\ldots,K$, satisfy \[ 0<\inf_{n\geq1}\inf_{1\leq k\leq K}\pi_{[k]}\leq\sup_{n\geq1}\sup_{1\leq k\leq K}\pi_{[k]}<1. \] Moreover, $\lim_{n\to\infty}P(\inf_{1\leq k\leq K}n_{[k]}\geq d+2)=1$, where $d$ denotes the dimension of $\bs{\xi}_{n}$. \end{enumerate} \end{assumption} \begin{assumption} $d$ and $K$ are fixed numbers. The sequence of (possibly) random functions $\bs{\xi}_{n}$ satisfies the following conditions: \begin{enumerate}[label=(\arabic*)] • There exists a sequence of non-stochastic functions $\bs{\xi}_{n}^{*}:\mathcal{X}\to\mathbb{R}^{d}$ ($n\geq1$) such that \begin{enumerate} • $\sup_{n\geq1}\mathbb{E}\left[\left\Vert \bs{\xi}_{n}^{*}(\bs X_{i})\right\Vert ^{2+\epsilon}\right]<\infty$, where $\epsilon$ is defined in Assumption (ref); • the minimal non-zero singular value of $\mathbb{E}\left[\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})^{\top}\mid B_{i}=k\right]$ is larger than some constant $c>0$ uniformly over $n\geq1$ and $k=1,\ldots,K$; • for every $k=1,\ldots,K$, \[ \frac{1}{n_{1[k]}}\sum_{i\in[k]}A_{i}\left\{ \bs{\xi}_{n}(\bs X_{i})-\bs{\xi}_{n}^{*}(\bs X_{i})\right\} -\frac{1}{n_{0[k]}}\sum_{i\in[k]}\left(1-A_{i}\right)\left\{ \bs{\xi}_{n}(\bs X_{i})-\bs{\xi}_{n}^{*}(\bs X_{i})\right\} =o_{P}(n^{-1/2}) \] and $\frac{1}{n_{[k]}}\sum_{i\in[k]}\left\Vert \bs{\xi}_{n}(\bs X_{i})-\bs{\xi}_{n}^{*}(\bs X_{i})\right\Vert ^{2}=o_{P}(1)$\textup{;} • $\liminf_{n\to\infty}(\varsigma_{\widetilde{Y}}^{2}-\varsigma_{\widetilde{Y}\mid\widetilde{\bs{\xi}}_{n}^{*}}^{2})>0$, where $\varsigma_{\widetilde{Y}}^{2}$ and $\varsigma_{\widetilde{Y}\mid\widetilde{\bs{\xi}}_{n}^{*}}^{2}$ are defined in Appendix (ref); \end{enumerate} • For every $k=1,\ldots,K$, \[ \left\Vert \left\{ \frac{1}{n}\sum_{i\in[k]}(A_{i}-\pi_{n[k]})^{2}\left\{ \bs{\xi}_{n}(\bs X_{i})-\overline{\bs{\xi}}_{n[k]}\right\} \left\{ \bs{\xi}_{n}(\bs X_{i})-\overline{\bs{\xi}}_{n[k]}\right\} ^{\top}\right\} ^{+}\right\Vert =O_{P}(1). \] \end{enumerate} \end{assumption} In this section, we focus on the case where the number of strata, $K$, is fixed, and leave the case where $K$ is diverging to Section (ref). Assumptions (ref) and (ref) are standard in the literature on statistical inference under CAR (see, e.g., bugni2018Inference,bannick2025General,bugni2019Inference,ma2022Regression,tu2024Unified,jiang2023Regressionadjusted and references therein). Assumption (ref) is satisfied for many randomization methods, such as stratified block randomization zelen1974Randomization, and Pocock and Simon’s minimization pocock1975Sequential. Assumption (ref)(1)(a) is a mild moment assumption on $\bs{\xi}_{n}^{*}$, the probability limit of $\bs{\xi}_{n}$. Assumption (ref)(1)(b) restricts that the minimal non-zero singular value of $\mathbb{E}\left[\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})^{\top}\mid B_{i}=k\right]$ is bounded away from zero. Importantly, this assumption allows the matrix $\mathbb{E}\left[\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})^{\top}\mid B_{i}=k\right]$ to be singular. This is reasonable, as in many applications, the components of $\bs{\xi}_{n}$ may share common information, leading to cases where the matrix $\mathbb{E}\left[\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})^{\top}\mid B_{i}=k\right]$ is indeed singular. Assumption (ref)(1)(c) requires that $\bs{\xi}_{n}$ converges to its probability limit $\bs{\xi}_{n}^{*}$. This assumption is similar to tu2024Unified, bannick2025General and jiang2023Regressionadjusted. When $\bs{\xi}_{n}$ is derived from linear regression, it simplifies to liu2023Lassoadjusted. When $\bs{\xi}_{n}$ is derived from local linear kernel regression, this assumption can be verified by tu2024Unified. If $\bs{\xi}_{n}$ is a general machine learning estimator that depends on the entire observed dataset, then verifying Assumption (ref)(1)(c) usually requires that the function class containing $\bs{\xi}_{n}$ satisfies the Donsker's condition (see bannick2025General), which may be violated if $\bs{\xi}_{n}$ is derived from complex machine learning algorithms. To bypass the Donsker's condition, one can use sample-splitting and/or cross-fitting (rafi2023Efficient,tu2024Unified,chernozhukov2018Double) techniques as discussed in Section (ref). By using cross-fitting, Assumption (ref)(1)(c) can be replaced by a mean square error convergence condition (see rafi2023Efficient and tu2024Unified), which is satisfied for many machine learning methods (see, e.g., farrell2021Deepa,chi2022Asymptotic,wager2018estimation,jiao2023Deep). Assumption (ref)(1)(d) ensures that the outcome $Y$ cannot be fully explained by $\bs{\xi}(\bs X)$. Assumption (ref)(2) is a technical requirement that ensures the estimated Moore--Penrose inverse remains well-behaved and does not diverge. This condition is necessary to address the discontinuity of the Moore--Penrose inverse operation (see e.g., stewart1977Perturbation), as Assumptions (ref)--(ref) and (ref)(1)(a)--(c) alone do not guarantee the convergence of the inverse. However, in the standard case where the limiting matrix $\mathbb{E}\left[\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})^{\top}\mid B_{i}=k\right]$ is non-singular, the Moore--Penrose inverse coincides with the standard matrix inverse. Since standard inversion is a continuous operation, Assumption (ref)(2) is then automatically implied by the preceding assumptions. Recall that $h_{a[k]}^{*}(\bs X)=\mathbb{E}\left[Y(a)\mid\bs X,B=k\right]$ and let $\widetilde{h}_{a[k]}^{*}(\bs X):=h_{a[k]}^{*}(\bs X)-\mathbb{E}\left[h_{a[k]}^{*}(\bs X)\mid B=k\right]$ and $\widetilde{\bs{\xi}}_{n}^{*}(\bs X):=\bs{\xi}_{n}^{*}(\bs X)-\mathbb{E}\left[\bs{\xi}_{n}^{*}(\bs X)\mid B=k\right]$. We have the following theorem. \begin{thm} Suppose that Assumptions (ref)--(ref) hold and $D(v)=(v-1)^{2}/2$. Then \[ \frac{\sqrt{n}\left(\widehat{\tau}_{\mathrm{cal}}-\tau\right)}{\sqrt{\varsigma_{H}^{2}+\varsigma_{\widetilde{Y}}^{2}-\varsigma_{\widetilde{Y}\mid\widetilde{\bs{\xi}}_{n}^{*}}^{2}}}\overset{d}{\to} N(0,1)\text{ and }\widehat{\varsigma}_{H}^{2}+\widehat{\varsigma}_{\widetilde{Y}}^{2}-\widehat{\varsigma}_{\widetilde{Y}\mid\widetilde{\bs{\xi}}_{n}^{*}}^{2}=\varsigma_{H}^{2}+\varsigma_{\widetilde{Y}}^{2}-\varsigma_{\widetilde{Y}\mid\widetilde{\bs{\xi}}_{n}^{*}}^{2}+o_{P}(1), \] where the definitions of the asymptotic variances $\varsigma_{H}^{2}$, $\varsigma_{\widetilde{Y}}^{2}$, $\varsigma_{\widetilde{Y}\mid\widetilde{\bs{\xi}}_{n}^{*}}^{2}$ and their estimates $\widehat{\varsigma}_{H}^{2},\widehat{\varsigma}_{\widetilde{Y}}^{2}$ and $\widehat{\varsigma}_{\widetilde{Y}\mid\widetilde{\bs{\xi}}_{n}^{*}}^{2}$ can be found in Appendix (ref). Moreover, if for every $k=1,\ldots,K$, there exists a non-stochastic vector $\bs{\alpha}_{k}\in\mathbb{R}^{d}$ such that \[ \left\{ \sqrt{\frac{1-\pi_{[k]}}{\pi_{[k]}}}\widetilde{h}_{1[k]}^{*}(\bs X)+\sqrt{\frac{\pi_{[k]}}{1-\pi_{[k]}}}\widetilde{h}_{0[k]}^{*}(\bs X)\right\} \1(B=k)=\bs{\alpha}_{k}^{\top}\widetilde{\bs{\xi}}_{n}^{*}(\bs X)\1(B=k), \] then \[ \frac{\sqrt{n}\left(\widehat{\tau}_{\mathrm{cal}}-\tau\right)}{\sqrt{\varsigma_{H}^{2}+\varsigma_{\widetilde{Y}}^{2}-\varsigma_{\widetilde{Y}\mid\widetilde{h}^{*}}^{2}}}\overset{d}{\to} N(0,1), \] and the asymptotic variance $\varsigma_{H}^{2}+\varsigma_{\widetilde{Y}}^{2}-\varsigma_{\widetilde{Y}\mid\widetilde{h}^{*}}^{2}$ matches the semiparametric efficiency bound developed in rafi2023Efficient. \end{thm} In Theorem (ref), we focus on the case that $D(v)=(v-1)^{2}/2$; the general $D(v)$ will be handled in Section (ref). Theorem (ref) shows that our calibration estimator is asymptotically normally distributed and its asymptotic variance can be consistently estimated. Moreover, as long as there exists a linear combination of $\bs{\xi}_{n}^{*}(\bs X)\1(B=k)$ that equals $\bigl\{\sqrt{\frac{1-\pi_{[k]}}{\pi_{[k]}}}h_{1[k]}^{*}(\bs X)+\sqrt{\frac{\pi_{[k]}}{1-\pi_{[k]}}}h_{0[k]}^{*}(\bs X)\bigr\}\1(B=k)$, our calibration estimator is semiparametric efficient. To the best of our knowledge, this is a new condition for achieving the efficiency bound in the literature. Existing work (e.g., tu2024Unified) requires that both $h_{1[k]}^{*}(\bs X)$ and $h_{0[k]}^{*}(\bs X)$ be consistently estimated. In contrast, we only require that the linear combination of these functions be consistently estimated, which is weaker. This condition aligns with bai2022Optimality, which states that optimal stratification can be achieved by knowing the function $\left\{ \sqrt{\frac{1-\pi_{[k]}}{\pi_{[k]}}}h_{1[k]}^{*}(\bs X)+\sqrt{\frac{\pi_{[k]}}{1-\pi_{[k]}}}h_{0[k]}^{*}(\bs X)\right\} $. In this context, the optimal stratification refers to the stratification method under which the difference-in-means estimator has the smallest mean squared error. Since the components of $\bs{\xi}_{n}(\bs X)$ can be chosen as different estimates of $\bigl\{\sqrt{\frac{1-\pi_{[k]}}{\pi_{[k]}}}h_{1[k]}^{*}(\bs X)+\sqrt{\frac{\pi_{[k]}}{1-\pi_{[k]}}}h_{0[k]}^{*}(\bs X)\bigr\}$ from various pre-specified models, the calibration estimator is \textit{multiply efficient}: it remains semiparametric efficient as long as at least one of these models is correctly specified. In contrast, AIPW-based covariate adjustment methods lack this kind of multiple efficiency. Theorem (ref) implies that the calibration estimator ((ref)) possesses two theoretical advantages. First, Assumption (ref) imposes no restriction on the choice of randomization scheme, which may include simple randomization, stratified block randomization, or minimization. Since the asymptotic distribution of the calibration estimator is invariant to the randomization method, it satisfies the property referred to in the literature as “universal applicability” (bannick2025General,ye2023Better). Second, because $\varsigma_{\widetilde{Y}\mid\widetilde{\bs{\xi}}_{n}^{*}}^{2}\geq0$, Theorem (ref), together with ma2022Regression, implies that $\widehat{\tau}_{\mathrm{cal}}$ is always at least as efficient as the stratified difference-in-means estimator $\widehat{\tau}_{\mathrm{sdim}}$, which corresponds to the special case where $\widehat{w}_{i}=1$ for all $i=1,\ldots,n$ in ((ref)). Hence, the calibration estimator guarantees efficiency gains. More generally, we have the following efficiency comparison result. \begin{thm} Rewrite the estimator $\widehat{\tau}_{\mathrm{cal}}$ in ((ref)) as $\widehat{\tau}_{\mathrm{cal}}(\bs{\xi}_{n})$. Let $\mathbf{\Lambda}$ be a $\widetilde{d}\times d$ non-stochastic matrix, where $\widetilde{d}\geq1$ is fixed. Suppose that Assumptions (ref)--(ref) hold for both $\bs{\xi}_{n}$ and $\mathbf{\Lambda}\bs{\xi}_{n}$. Then $\widehat{\tau}_{\mathrm{cal}}(\bs{\xi}_{n})$ and $\widehat{\tau}_{\mathrm{cal}}(\mathbf{\Lambda}\bs{\xi}_{n})$ satisfy \[ \frac{\sqrt{n}\left(\widehat{\tau}_{\mathrm{cal}}(\bs{\xi}_{n})-\tau\right)}{\sqrt{\varsigma_{H}^{2}+\varsigma_{\widetilde{Y}}^{2}-\varsigma_{\widetilde{Y}\mid\widetilde{\bs{\xi}}_{n}^{*}}^{2}}}\overset{d}{\to} N(0,1)\text{ and }\frac{\sqrt{n}\left(\widehat{\tau}_{\mathrm{cal}}(\mathbf{\Lambda}\bs{\xi}_{n})-\tau\right)}{\sqrt{\varsigma_{H}^{2}+\varsigma_{\widetilde{Y}}^{2}-\varsigma_{\widetilde{Y}\mid\mathbf{\Lambda}\widetilde{\bs{\xi}}_{n}^{*}}^{2}}}\overset{d}{\to} N(0,1), \] respectively. Furthermore, the asymptotic variances satisfy $\varsigma_{H}^{2}+\varsigma_{\widetilde{Y}}^{2}-\varsigma_{\widetilde{Y}\mid\widetilde{\bs{\xi}}_{n}^{*}}^{2}\leq\varsigma_{H}^{2}+\varsigma_{\widetilde{Y}}^{2}-\varsigma_{\widetilde{Y}\mid\mathbf{\Lambda}\widetilde{\bs{\xi}}_{n}^{*}}^{2}.$ \end{thm} Theorem (ref) implies that incorporating additional elements into $\bs{\xi}_{n}$ satisfies a no-harm property; i.e., it is guaranteed to improve or at least maintain efficiency. Consequently, from a theoretical perspective, it appears optimal to include as many elements as possible in $\bs{\xi}_{n}$. However, as the dimension of $\bs{\xi}_{n}$ becomes large, it is no longer appropriate to treat the dimension $d$ as fixed. We will address the scenario where the dimension of $\bs{\xi}_{n}$ diverges in Section (ref). \section{Extensions} \subsection{The dimension of $\protect\bs{\xi}_{n}(\cdot)$ and the number of strata are diverging} In this section, we allow the dimension of $\bs{\xi}_{n}(\cdot)$ to grow with $n$, that is, $d=d_{n}\to\infty$ as $n\to\infty$. Additionally, we also allow the number of strata to increase with $n$, i.e., $K=K_{n}\to\infty$ as $n\to\infty$, which is commonly encountered in many applications. We begin by stating an assumption. \begin{assumption} The sequence of random functions $\bs{\xi}_{n}(\cdot)$ satisfies the following conditions. \begin{enumerate}[label=(\arabic*)] • There exists a sequence of non-stochastic functions $\bs{\xi}_{n}^{*}(\cdot):\mathcal{X}\to\mathbb{R}^{d}$ ($n\geq1$) such that \begin{enumerate} • $K^{2}r_{n}^{2}/n\to0$, $K^{2}\zeta_{n}^{2}r_{n}/n\to0$ and $\max\left\{ \zeta_{n}^{2}\log(2Kd),r_{n}\log^{2}(2Kd)\right\} /(n\inf_{1\leq k\leq K}p_{[k]})\to0$ as $n\to\infty$, where $r_{n}:=\max_{1\leq k\leq K}\mathrm{rank}\left\{ \mathbb{E}\left[\bs{\xi}_{n}^{*}(\bs X_{i})\bs{\xi}_{n}^{*}(\bs X_{i})^{\top}\mid B_{i}=k\right]\right\} $ and $\zeta_{n}:=\sup_{\bs x\in\mathcal{X}}\left\Vert \bs{\xi}_{n}^{*}(\bs x)\right\Vert $; • the maximal singular value of $\mathbb{E}\left[\bs{\xi}_{n}^{*}(\bs X_{i})\bs{\xi}_{n}^{*}(\bs X_{i})^{\top}\mid B_{i}=k\right]$ is smaller than some finite constant $C>0$ uniformly over $n\geq1$ and $k=1,\ldots,K$; the minimal non-zero singular value of $\mathbb{E}\left[\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})^{\top}\mid B_{i}=k\right]$ is larger than some constant $c>0$ uniformly over $n\geq1$ and $k=1,\ldots,K$; • it holds that \begin{align*} & \sup_{1\leq k\leq K}\left\Vert \frac{1}{n_{1[k]}}\sum_{i\in[k]}A_{i}\left\{ \bs{\xi}_{n}(\bs X_{i})-\bs{\xi}_{n}^{*}(\bs X_{i})\right\} -\frac{1}{n_{0[k]}}\sum_{i\in[k]}\left(1-A_{i}\right)\left\{ \bs{\xi}_{n}(\bs X_{i})-\bs{\xi}_{n}^{*}(\bs X_{i})\right\} \right\Vert \\ = & o_{P}(n^{-1/2}) \end{align*} and $\sup_{1\leq k\leq K}\frac{1}{n_{[k]}}\sum_{i\in[k]}\left\Vert \bs{\xi}_{n}(\bs X_{i})-\bs{\xi}_{n}^{*}(\bs X_{i})\right\Vert ^{2}=O_{P}(n^{-1/2})$\textup{;} • $\liminf_{n\to\infty}\left(\varsigma_{\widetilde{Y}}^{2}-\varsigma_{\widetilde{Y}\mid\widetilde{\bs{\xi}}_{n}^{*}}^{2}\right)>0$. \end{enumerate} • $\sup_{1\leq k\leq K}\left\Vert \left\{ \frac{1}{n_{[k]}}\sum_{i\in[k]}(A_{i}-\pi_{n[k]})^{2}\left\{ \bs{\xi}_{n}(\bs X_{i})-\overline{\bs{\xi}}_{n[k]}\right\} \left\{ \bs{\xi}_{n}(\bs X_{i})-\overline{\bs{\xi}}_{n[k]}\right\} ^{\top}\right\} ^{+}\right\Vert =O_{P}(1)$. \end{enumerate} \end{assumption} Assumption (ref) is very similar to Assumption (ref), except that we allow both the dimension of $\bs{\xi}_{n}$ and the number of strata to diverge. Assumption (ref)(1)(a) restricts the growth rate of $d$ and $K$. Suppose that the matrix $\mathbb{E}\left[\bs{\xi}_{n}^{*}(\bs X_{i})\bs{\xi}_{n}^{*}(\bs X_{i})^{\top}\mid B_{i}=k\right]$ is non-singular and the elements in $\bs{\xi}_{n}^{*}(\bs x)$ are uniformly bounded, we have $r_{n}=d$ and $\zeta_{n}=\sqrt{d}$. Consider the case that the strata have equal size, i.e., $p_{[k]}=1/K$. Then a sufficient condition for Assumption (ref)(1)(a) is $Kd=o(\sqrt{n})$, meaning that we allow the product of the dimension of $\bs{\xi}_{n}$ and the number of strata to grow, but at a rate slower than $\sqrt{n}$. When $K$ is fixed, this condition is consistent with a similar result for the regression-adjusted ATE estimator presented by lei2021Regression under complete randomization and finite-population asymptotics. To further relax the growth rate constraint on $d$, lei2021Regression, lu2025Debiased, and gu2025assumptionleancovariateadjustmentcovariate proposed debiased regression-adjusted estimators that allow $d$ to grow faster than $\sqrt{n}$. In contrast, jiang2025Adjustments achieved the same goal by assuming a correctly specified linear relationship between the potential outcomes and covariates. It is also possible to debias the calibration estimator ((ref)) in a way that would relax Assumption (ref)(1)(a). However, addressing this is beyond the scope of the present paper and will be left for future research. Assumption (ref)(1)(b) requires that largest and smallest non-zero singular values of $\mathbb{E}\left[\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})^{\top}\mid B_{i}=k\right]$ are bounded away from zero and infinity, while still allowing this matrix to be singular. Assumption (ref)(1)(c) strengthens Assumption (ref)(1)(c). Under cross-fitting, Assumption (ref)(1)(c) requires that $\bs{\xi}_{n}$ converges to its probability limit $\bs{\xi}_{n}^{*}$ at rate $n^{-1/4}$ uniformly over $k=1,\ldots,K$. Such a rate is achievable for many modern machine learning methods (see, e.g., farrell2021Deepa,jiao2023Deep). Assumption (ref)(2) parallels Assumption (ref)(2), except that it is required to hold uniformly over $k=1,\ldots,K$. We are now ready to state the following theorem. \begin{thm} Suppose that Assumptions (ref)--(ref) and (ref) hold and $D(v)=(v-1)^{2}/2$. Then \[ \frac{\sqrt{n}\left(\widehat{\tau}_{\mathrm{cal}}-\tau\right)}{\sqrt{\varsigma_{H}^{2}+\varsigma_{\widetilde{Y}}^{2}-\varsigma_{\widetilde{Y}\mid\widetilde{\bs{\xi}}_{n}^{*}}^{2}}}\overset{d}{\to} N(0,1)\text{ and }\widehat{\varsigma}_{H}^{2}+\widehat{\varsigma}_{\widetilde{Y}}^{2}-\widehat{\varsigma}_{\widetilde{Y}\mid\widetilde{\bs{\xi}}_{n}^{*}}^{2}=\varsigma_{H}^{2}+\varsigma_{\widetilde{Y}}^{2}-\varsigma_{\widetilde{Y}\mid\widetilde{\bs{\xi}}_{n}^{*}}^{2}+o_{P}(1). \] Moreover, when $K$ is a fixed number, if for every $k=1,\ldots,K$, there exists a non-stochastic vector $\bs{\alpha}_{k}\in\mathbb{R}^{d}$ such that \[ \left\{ \sqrt{\frac{1-\pi_{[k]}}{\pi_{[k]}}}\widetilde{h}_{1[k]}^{*}(\bs X)+\sqrt{\frac{\pi_{[k]}}{1-\pi_{[k]}}}\widetilde{h}_{0[k]}^{*}(\bs X)\right\} \1(B=k)=\bs{\alpha}_{k}^{\top}\widetilde{\bs{\xi}}_{n}^{*}(\bs X)\1(B=k), \] then \[ \frac{\sqrt{n}\left(\widehat{\tau}_{\mathrm{cal}}-\tau\right)}{\sqrt{\varsigma_{H}^{2}+\varsigma_{\widetilde{Y}}^{2}-\varsigma_{\widetilde{Y}\mid\widetilde{h}^{*}}^{2}}}\overset{d}{\to} N(0,1), \] and the asymptotic variance $\varsigma_{H}^{2}+\varsigma_{\widetilde{Y}}^{2}-\varsigma_{\widetilde{Y}\mid\widetilde{h}^{*}}^{2}$ matches the semiparametric efficiency bound developed in rafi2023Efficient. \end{thm} Theorem (ref) shows that, even when both the dimension of $\bs{\xi}_{n}$ and the number of strata diverge, the calibration estimator ((ref)) remains asymptotically normal. It is important to note that, because we allow the number of strata to increase with the sample size, Assumption (ref) becomes non-trivial. xin2024inferencecovariateadaptiverandomizationstrata studied inference under CAR with a growing number of strata, and our Assumption (ref) corresponds directly to their Assumptions (B1)--(B3). Moreover, xin2024inferencecovariateadaptiverandomizationstrata derived more primitive conditions under which Assumption (ref) holds in the cases of simple randomization, stratified adaptive biased-coin randomization and stratified block randomization. For further discussion of Assumption (ref), we refer readers to xin2024inferencecovariateadaptiverandomizationstrata. \subsection{General discrepancy measure $D(v)$} In Theorem (ref), our analysis was limited to the quadratic discrepancy measure $D(v)=(v-1)^{2}/2$. In this section, we extend the analysis to a general class of discrepancy measures $D(v)$, which can be theoretically shown to exhibit smaller second-order bias. To proceed, we impose the following mild regularity condition on $D(v)$. \begin{assumption} Let $D^{\prime}$ be the derivative of $D$ and $(D^{\prime})^{-1}$ is the inverse function of $D^{\prime}$. Define $\rho(v):=D\left\{ (D^{\prime})^{-1}(-v)\right\} +v\cdot(D^{\prime})^{-1}(-v)$. We assume $\rho(v)$ is concave and three times continuously differentiable, $\rho^{\prime}(0)=1$, $-\infty<\rho^{\prime\prime}(0)<0$ and there exist constants $\delta>0$ and $0\leq C_{\rho}<\infty$ such that $\left|\rho^{\prime\prime\prime}(v)-\rho^{\prime\prime\prime}(0)\right|\leq C_{\rho}\left|v\right|$ for all $\left|v\right|\leq\delta$. \end{assumption} Assumption (ref) holds for a wide range of commonly used discrepancy measures, and Table (ref) provides several popular examples. The case that $D(v)=v\log v-v$ is related to the exponential tilting estimator studied by kitamura1997Informationtheoretic, and $D(v)=v-\log v$ is related to empirical likelihood estimator studied by qin1994Empirical. \begin{table}[!tbh] \caption{Different choices of $D(v)$ and their corresponding $\rho(v)$.} \begin{tabular}{ccccc} \toprule $D(v)$ & $\rho(v)$ & $\rho^{\prime}(v)$ & $\rho^{\prime\prime}(0)$ & $\rho^{\prime\prime\prime}(0)$\tabularnewline \midrule \midrule $(v-1)^{2}/2$ & $-v^{2}/2+v$ & $-v+1$ & $-1$ & $0$\tabularnewline \midrule $v\log v-v$ & $-e^{-v}$ & $e^{-v}$ & $-1$ & $1$\tabularnewline \midrule $v-\log v$ & $1+\log(1+v)$ & $\frac{1}{1+v}$ & $-1$ & $2$\tabularnewline \bottomrule \end{tabular} \end{table} Since our goal is to characterize the second-order bias of the calibration estimator, we require a strengthened version of Assumption (ref). \begin{assumption} \begin{enumerate}[label=(\arabic*)] • The sequence of random functions $\bs{\xi}_{n}(\cdot)$ satisfies $\frac{1}{n}\sum_{i=1}^{n}\left\Vert \bs{\xi}_{n}(\bs X_{i})\right\Vert ^{4}=O_{P}(1)$; • The non-stochastic functions $\bs{\xi}_{n}^{*}(\cdot):\mathcal{X}\to\mathbb{R}^{d}$ ($n\geq1$) in Assumption (ref) satisfy the following conditions: $\sup_{n\geq1}\mathbb{E}\left[\left\Vert \bs{\xi}_{n}^{*}(\bs X_{i})\right\Vert ^{4}\right]<\infty$, and for each $k=1,\ldots,K$, \[ \left\Vert \frac{1}{n_{1[k]}}\sum_{i\in[k]}A_{i}\left\{ \bs{\xi}_{n}(\bs X_{i})-\bs{\xi}_{n}^{*}(\bs X_{i})\right\} -\frac{1}{n_{0[k]}}\sum_{i\in[k]}\left(1-A_{i}\right)\left\{ \bs{\xi}_{n}(\bs X_{i})-\bs{\xi}_{n}^{*}(\bs X_{i})\right\} \right\Vert =o_{P}(n^{-1}), \] and $\frac{1}{n_{[k]}}\sum_{i\in[k]}\left\Vert \bs{\xi}_{n}(\bs X_{i})-\bs{\xi}_{n}^{*}(\bs X_{i})\right\Vert ^{2}=O_{P}(\Delta_{n}^{2})$\textup{,} where $\Delta_{n}\to0$ as $n\to\infty$ is a sequence of real numbers. \end{enumerate} \end{assumption} We have the following theorem. \begin{thm} Suppose that Assumptions (ref)--(ref) and (ref) hold. If $\sup_{n\geq1}\sup_{1\leq i\leq n}\mathbb{E}\left[\left\Vert \bs{\xi}_{n}(\bs X_{i})\right\Vert ^{2+\epsilon}\right]<\infty$ for some $\epsilon>0$ and $\mathbb{E}\left[\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})^{\top}\mid B_{i}=k\right]$ is non-singular for every $k=1,\ldots,K$, then the conclusion of Theorem (ref) holds. If, in addition, Assumption (ref) holds and $\mathbb{E}\left[Y_{i}(a)^{4}\right]<\infty$ for all $a=0,1$, then \[ \widehat{\tau}_{\mathrm{cal}}-\tau=\frac{1}{n}\sum_{i=1}^{n}\psi_{1,i}+\frac{1}{n}\psi_{2}+o_{P}(n^{-1})+O_{P}(\Delta_{n}n^{-1/2}), \] where $\psi_{1,i}=\sum_{k=1}^{K}\left\{ \left(\frac{A_{i}}{\pi_{n[k]}}-\frac{1-A_{i}}{1-\pi_{n[k]}}\right)\1(B_{i}=k)\cdot Y_{i}-\tau-\bs{\beta}_{[k],\mathcal{C}_{n}}^{\top}\bs{\Xi}_{i,[k]}^{*}\right\} $ with $\mathbb{E}\left[\frac{1}{n}\sum_{i=1}^{n}\psi_{1,i}\right]=0$ and $\mathbb{E}\left[\psi_{2}\right]=\left\{ \frac{\rho^{\prime\prime\prime}(0)}{2\rho^{\prime\prime}(0)^{2}}-1\right\} \sum_{k=1}^{K}\mathrm{tr}\left\{ \mathbb{E}\left[\left(\mathbf{\Sigma}_{[k]}^{\mathcal{C}_{n}}\right)^{-1}\mathbf{\Sigma}_{\bs{\Xi}\bs{\Xi}\epsilon[k]}^{\mathcal{C}_{n}}\right]\right\} $\textup{.} The definitions of $\bs{\beta}_{[k],\mathcal{C}_{n}}^{\top}\bs{\Xi}_{i,[k]}^{*}$ and $\left(\mathbf{\Sigma}_{[k]}^{\mathcal{C}_{n}}\right)^{-1}\mathbf{\Sigma}_{\bs{\Xi}\bs{\Xi}\epsilon[k]}^{\mathcal{C}_{n}}$ are provided in Appendix (ref). \end{thm} The first conclusion of Theorem (ref) implies that, provided that Assumption (ref) holds, different choices of $D(v)$ lead to calibration estimators with the same asymptotic distribution. We also note that this conclusion requires the non-singularity of $\mathbb{E}\left[\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})^{\top}\mid B_{i}=k\right]$ in order to ensure that ((ref)) admits a unique solution asymptotically. In the special case where $D(v)=(v-1)^{2}/2$, such a condition is not necessary. This is because, under the quadratic discrepancy, the weights $\widehat{w}_{i},i=1,\ldots,n$, admit a closed-form expression. Even if (ref) admits infinitely many solutions, one can always select the solution obtained via the Moore-Penrose inverse (see ((ref)) in the proof of Theorem (ref) in Appendix (ref)). In contrast, when $D(v)$ is a general discrepancy measure, the weights $\widehat{w}_{i}$ no longer have an explicit form. In this case, if ((ref)) admits infinitely many solutions, it is unclear which solution an optimization algorithm would converge to. In practice, a regularization term could be introduced to enforce the uniqueness of the solution, although this is beyond the scope of the present paper. Additionally, if $\mathbb{E}\left[\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})\widetilde{\bs{\xi}}_{n}^{*}(\bs X_{i})^{\top}\mid B_{i}=k\right]$ is singular, one could apply singular value decomposition (SVD) to eliminate the redundant components of $\bs{\xi}_{n}$ and use the reduced version of $\bs{\xi}_{n}$ to compute the calibration estimator ((ref)). We leave the rigorous theoretical development of this strategy for future investigation. The second conclusion of Theorem (ref) provides a characterization of the second-order bias under the CAR design, which, to the best of our knowledge, is a novel contribution to the literature. Unlike existing studies on second-order bias (newey2004Higher,tan2014Secondorder), the observed samples under the CAR design are not i.i.d. The proof of this result relies on a conditional argument: we first examine the asymptotic expansion of the calibration estimator conditional on the treatment and stratum indicators, and then apply the law of iterated expectation. From Theorem (ref), we see that when $\Delta_{n}=o(n^{-1/2})$, taking $D(v)=v-\log v$ gives the calibration estimator with zero second-order bias $\mathbb{E}\left[\psi_{2}\right]$. This aligns with the findings in i.i.d. settings, where empirical likelihood-based estimators exhibit smaller second-order bias (newey2004Higher,tan2014Secondorder). \section{Simulation studies} We evaluate the performance of our calibration estimators using Monte Carlo simulations. For $a\in\{0,1\}$ and $1\leq i\leq n$, the potential outcomes are generated as follows: \[ Y_{i}(a)=g_{a}(\bs X_{i})+\epsilon_{a,i}\ \ i=1,\ldots,n,\ a\in\{0,1\}, \] where $\bs X_{i}$, $\epsilon_{a,i}$ and $g_{a}(\cdot)$ will be specified in each model and the triplet $(\bs X_{i},\epsilon_{0,i},\epsilon_{1,i})$ for $1\leq i\leq n$ is i.i.d. We present simulation results for the estimators under three randomization methods: simple randomization, stratified block randomization, and minimization (pocock1975Sequential). The sample size $n$ varies across 500, 1000, and 2000, and the number of covariates is set to be $p=30$. For stratified block randomization, we use a block size of 6. In the minimization method, a biased-coin probability of 0.75 and equal weights are employed. The cross-fitting technique is applied with two folds as demonstrated in Section (ref). Unless otherwise specified, the calibration estimator is obtained by setting $D(v)=(v-1)^{2}/2$. We compare nine estimators: \textbf{(i)} \texttt{cal_rf}: For each stratum $k$, we use random forest to estimate the conditional mean function $g_{0}(\cdot)$ and $g_{1}(\cdot)$, denoting their estimates as $\widehat{g}_{0k}^{\text{rf}}(\cdot)$ and $\widehat{g}_{1k}^{\text{rf}}(\cdot)$. Then \texttt{cal_rf }is obtained by taking $\bs{\xi}_{n}(\bs X_{i})=(\sum_{k=1}^{K}\widehat{g}_{0k}^{\text{rf}}(\bs X_{i})\1(B_{i}=k),\sum_{k=1}^{K}\widehat{g}_{1k}^{\text{rf}}(\bs X_{i})\1(B_{i}=k))^{\top}$ in our calibration estimator ((ref)); \textbf{(ii)} \texttt{cal_nn}: Defined similarly to \texttt{cal_rf}, but with the random forest estimates replaced by neural network estimates $\widehat{g}_{0k}^{\text{rf}}(\cdot)$ and $\widehat{g}_{1k}^{\text{rf}}(\cdot)$; \textbf{(iii)} \texttt{cal_rfnn}: This estimator is obtained by combining the random forest and neural network estimates. It is obtained by taking $\bs{\xi}_{n}(\bs X_{i})=(\sum_{k=1}^{K}\widehat{g}_{0k}^{\sharp}(\bs X_{i})\1(B_{i}=k),\sum_{k=1}^{K}\widehat{g}_{1k}^{\sharp}(\bs X_{i})\1(B_{i}=k):\sharp\in\{\text{rf},\text{nn}\})^{\top}$ in our calibration estimator ((ref)); \textbf{(iv)} \texttt{cal_rflin}: Defined similarly to \texttt{cal_rfnn}, but obtained by combining the random forest and linear regression estimates. \textbf{(v)}\texttt{ cal_rf_g}: Take $\bs{\xi}_{n}(\bs X_{i})=(\widehat{g}_{0k}^{\text{rf}}(\bs X_{i}),\widehat{g}_{1k}^{\text{rf}}(\bs X_{i}):1\leq k\leq K)^{\top}$ and use it in our calibration estimator ((ref)); \textbf{(vi)} \texttt{cal_nn_g}: Defined similarly to \texttt{cal_rf_g}, but with neural network estimates; \textbf{(vii)} \texttt{cal_lin_EL}: Defined similarly to \texttt{cal_rf}, but replacing the random forest estimates with linear regression estimates and taking $D(v)=v-\log v$;\textbf{ (viii)} \texttt{aipw_rf}: We use AIPW-base method (tu2024Unified) with the random forest estimates $\widehat{g}_{0k}^{\text{rf}}(\bs X_{i})$ and $\widehat{g}_{1k}^{\text{rf}}(\bs X_{i})$ to adjust the covariates; \textbf{(ix)} \texttt{aipw_nn}: Defined similarly to \texttt{aipw_rf}, but replacing the random forest estimates with neural network estimates $\widehat{g}_{0k}^{\text{nn}}(\cdot)$ and $\widehat{g}_{1k}^{\text{nn}}(\cdot)$; \textbf{(x)} \texttt{aipw_lin}: Defined similarly to \texttt{aipw_rf}, but replacing the random forest estimates with linear regression estimates; \textbf{(xi)} \texttt{sdim}: The stratified difference-in-means estimator $\widehat{\tau}_{\mathrm{sdim}}$. For each estimator, we report the following metrics based on 300 replications: absolute bias, empirical standard deviation (SD), average estimated standard error (SE), and the empirical coverage probability (CP) of the 95% confidence intervals. \textbf{Model 1.} Model 1 imposes linear models for the conditional mean functions $g_{0}(\bs X)$ and $g_{1}(\bs X)$. In this model, we set \begin{align*} g_{0}(\bs X_{i}) & =\mu_{0}+\sum_{j=1}^{4}\beta_{0j}X_{ij}\ \text{ and }\ g_{1}(\bs X_{i})=\mu_{1}+\sum_{j=1}^{4}\beta_{1j}X_{ij}, \end{align*} with $\mu_{0}=1$, $\mu_{1}=4$, $(\beta_{01},\ldots,\beta_{04})=(75,35,125,80)$, and $(\beta_{11},\ldots,\beta_{14})=(100,80,60,40)$. Additionally, the variables are specified as follows: $\epsilon_{0,i}\sim N(0,1)$, $\epsilon_{1,i}\sim N(0,9)$, $X_{i1}\sim\text{Beta}(3,4)$, $X_{i2}\sim\text{Uniform}(-2,2)$, $X_{i3}$ takes values in $\{-1,1\}$ with equal probability, $X_{i4}$ takes values in $\{3,5\}$ with probabilities 0.6 and 0.4, respectively. These variables are independent of one another. The remaining variables $X_{i5},\dots,X_{ip}$ are independent of $X_{i1},\dots,X_{i4}$ and follow a multivariate normal distribution with zero mean and a covariance matrix where all off-diagonal elements are 0.2, while the diagonal elements are 1. The randomization variable is an additional variable taking values in $\{1,2,3,4\}$ with probabilities 0.2, 0.3, 0.3, and 0.2, respectively, and is independent of $X_{ij}$ for $j=1,\dots,p$. The simulation results for Model 1 are presented in Table (ref). Since both $g_{0}(\bs X_{i})$ and $g_{1}(\bs X_{i})$ are linear in the covariates $X_{ij}$, the linear regression-adjusted estimator, \texttt{aipw_lin}, performs optimally in large samples ($n=1000,2000$), as expected. The performance of \texttt{cal_rflin} is very close to that of \texttt{aipw_lin}, which aligns with the conclusion of Theorem (ref). However, in smaller samples ($n=500$), the \texttt{aipw_lin} estimator does not perform the best, as linear regression is sensitive to outliers. In contrast, \texttt{cal_rflin} maintains stable performance, suggesting that incorporating different estimates of $g_{0}(\bs X_{i})$ and $g_{1}(\bs X_{i})$ into the calibration estimator enhances its robustness. \begin{table}[!tbh] \caption{The comparison of the performance of different estimators under Model 1.} \resizebox{\textwidth}{!}{ \begin{threeparttable} \begin{centering} \begin{tabular}{clrrrrrrrrrrrr} \toprule \multirow{2}{*}{$n$} & \multirow{2}{*}{Estimator} & \multicolumn{4}{c}{Simple Rand.} & \multicolumn{4}{c}{Stratified Block Rand.} & \multicolumn{4}{c}{Minimization}\tabularnewline \cmidrule{3-14} & & Bias & SD & SE & CP & Bias & SD & SE & CP & Bias & SD & SE & CP\tabularnewline \midrule \multirow{11}{*}{500} & \texttt{cal_rf} & 0.12 & 6.85 & 7.39 & 0.973 & 0.12 & 6.74 & 7.31 & 0.960 & 0.36 & 6.98 & 7.33 & 0.953\tabularnewline & \texttt{cal_nn} & 0.33 & 11.38 & 11.03 & 0.960 & 0.65 & 11.01 & 11.05 & 0.950 & 0.63 & 11.47 & 11.00 & 0.933\tabularnewline & \texttt{cal_rfnn} & 0.47 & 6.90 & 7.33 & 0.960 & 0.18 & 6.90 & 7.24 & 0.960 & 0.12 & 7.03 & 7.27 & 0.950\tabularnewline & \texttt{cal_rflin} & 0.08 & 4.36 & 5.23 & 0.983 & 0.35 & 4.78 & 5.19 & 0.967 & 0.37 & 4.60 & 5.19 & 0.973\tabularnewline & \texttt{cal_rf_g} & 0.13 & 5.83 & 6.08 & 0.967 & 0.48 & 6.01 & 6.07 & 0.960 & 0.01 & 5.68 & 6.09 & 0.960\tabularnewline & \texttt{cal_nn_g} & 0.48 & 9.23 & 9.10 & 0.960 & 0.79 & 9.41 & 9.15 & 0.943 & 0.31 & 9.31 & 9.10 & 0.950\tabularnewline & \texttt{cal_lin_EL} & 0.10 & 4.25 & 5.26 & 0.980 & 0.32 & 4.49 & 5.21 & 0.983 & 0.14 & 4.39 & 5.22 & 0.987\tabularnewline & \texttt{aipw_rf} & 0.31 & 9.03 & 8.91 & 0.957 & 0.07 & 8.74 & 8.88 & 0.940 & 0.65 & 9.13 & 8.88 & 0.943\tabularnewline & \texttt{aipw_nn} & 0.29 & 14.17 & 13.55 & 0.923 & 0.74 & 14.45 & 13.54 & 0.910 & 1.22 & 13.98 & 13.52 & 0.950\tabularnewline & \texttt{aipw_lin} & 3.22 & 50.59 & 14.23 & 0.953 & 0.69 & 25.12 & 9.04 & 0.967 & 1.78 & 19.45 & 7.30 & 0.953\tabularnewline & \texttt{sdim} & 0.48 & 12.35 & 12.34 & 0.953 & 0.06 & 11.71 & 12.31 & 0.943 & 0.93 & 12.61 & 12.28 & 0.927\tabularnewline \midrule \multirow{11}{*}{1000} & \texttt{cal_rf} & 0.21 & 4.00 & 4.13 & 0.957 & 0.06 & 3.70 & 4.12 & 0.960 & 0.09 & 3.85 & 4.12 & 0.963\tabularnewline & \texttt{cal_nn} & 0.10 & 6.39 & 6.48 & 0.963 & 0.15 & 6.42 & 6.48 & 0.940 & 0.38 & 6.14 & 6.44 & 0.960\tabularnewline & \texttt{cal_rfnn} & 0.05 & 3.93 & 4.06 & 0.963 & 0.06 & 3.78 & 4.05 & 0.960 & 0.17 & 3.90 & 4.05 & 0.973\tabularnewline & \texttt{cal_rflin} & 0.17 & 2.94 & 3.28 & 0.973 & 0.24 & 2.84 & 3.27 & 0.973 & 0.05 & 3.14 & 3.27 & 0.957\tabularnewline & \texttt{cal_rf_g} & 0.24 & 3.68 & 3.66 & 0.940 & 0.05 & 3.45 & 3.67 & 0.963 & 0.08 & 3.67 & 3.66 & 0.950\tabularnewline & \texttt{cal_nn_g} & 0.40 & 4.48 & 4.73 & 0.967 & 0.18 & 4.83 & 4.76 & 0.950 & 0.55 & 4.74 & 4.74 & 0.950\tabularnewline & \texttt{cal_lin_EL} & 0.12 & 2.87 & 3.27 & 0.973 & 0.29 & 2.76 & 3.27 & 0.977 & 0.05 & 3.06 & 3.27 & 0.970\tabularnewline & \texttt{aipw_rf} & 0.07 & 5.11 & 5.24 & 0.967 & 0.04 & 5.25 & 5.23 & 0.947 & 0.16 & 5.27 & 5.24 & 0.950\tabularnewline & \texttt{aipw_nn} & 0.17 & 7.96 & 7.51 & 0.943 & 0.26 & 7.52 & 7.53 & 0.953 & 0.36 & 7.38 & 7.44 & 0.973\tabularnewline & \texttt{aipw_lin} & 0.08 & 2.81 & 2.90 & 0.943 & 0.32 & 2.74 & 2.91 & 0.960 & 0.05 & 3.03 & 2.91 & 0.940\tabularnewline & \texttt{sdim} & 0.11 & 8.27 & 8.71 & 0.953 & 0.10 & 8.70 & 8.67 & 0.950 & 0.27 & 8.61 & 8.67 & 0.953\tabularnewline \midrule \multirow{11}{*}{2000} & \texttt{cal_rf} & 0.10 & 2.49 & 2.53 & 0.960 & 0.02 & 2.32 & 2.53 & 0.960 & 0.11 & 2.50 & 2.53 & 0.967\tabularnewline & \texttt{cal_nn} & 0.27 & 2.50 & 2.68 & 0.960 & 0.01 & 2.62 & 2.68 & 0.963 & 0.05 & 2.64 & 2.67 & 0.950\tabularnewline & \texttt{cal_rfnn} & 0.17 & 2.23 & 2.33 & 0.960 & 0.04 & 2.17 & 2.33 & 0.963 & 0.04 & 2.32 & 2.33 & 0.953\tabularnewline & \texttt{cal_rflin} & 0.12 & 2.04 & 2.19 & 0.967 & 0.04 & 2.06 & 2.19 & 0.973 & 0.02 & 2.11 & 2.19 & 0.970\tabularnewline & \texttt{cal_rf_g} & 0.13 & 2.37 & 2.38 & 0.950 & 0.07 & 2.22 & 2.38 & 0.957 & 0.14 & 2.33 & 2.38 & 0.953\tabularnewline & \texttt{cal_nn_g} & 0.28 & 2.09 & 2.23 & 0.957 & 0.15 & 2.19 & 2.24 & 0.957 & 0.05 & 2.18 & 2.23 & 0.967\tabularnewline & \texttt{cal_lin_EL} & 0.16 & 2.01 & 2.19 & 0.970 & 0.00 & 2.03 & 2.19 & 0.970 & 0.02 & 2.08 & 2.19 & 0.973\tabularnewline & \texttt{aipw_rf} & 0.12 & 3.04 & 3.05 & 0.957 & 0.15 & 3.00 & 3.05 & 0.953 & 0.20 & 3.10 & 3.05 & 0.947\tabularnewline & \texttt{aipw_nn} & 0.29 & 2.77 & 2.91 & 0.967 & 0.04 & 2.93 & 2.89 & 0.927 & 0.06 & 2.81 & 2.88 & 0.953\tabularnewline & \texttt{aipw_lin} & 0.16 & 1.99 & 2.06 & 0.967 & 0.00 & 2.02 & 2.06 & 0.957 & 0.03 & 2.06 & 2.05 & 0.963\tabularnewline & \texttt{sdim} & 0.22 & 5.87 & 6.13 & 0.967 & 0.41 & 6.28 & 6.13 & 0.940 & 0.37 & 6.12 & 6.13 & 0.953\tabularnewline \bottomrule \end{tabular} \end{centering} \begin{tablenotes}[flushleft] • \textit{Abbreviations:} Rand., Randomization; SD, standard deviation, SE: standard error; CP, coverage probability. \end{tablenotes} \end{threeparttable} } \end{table} \textbf{Model 2.} Model 2 imposes additive but nonlinear models for the conditional mean functions $g_{0}(\bs X)$ and $g_{1}(\bs X)$. In this model, we set \begin{align*} g_{0}(\bs X_{i}) & =\mu_{0}+\beta_{01}\log(X_{i1}+1)+\beta_{02}X_{i2}^{2}+\beta_{03}\exp(X_{i3})+\beta_{04}/(X_{i4}+3)\\ g_{1}(\bs X_{i}) & =\mu_{1}+\beta_{11}\exp(X_{i1}+2)+\beta_{12}/(X_{i1}+1)+\beta_{13}X_{i2}^{2}, \end{align*} with $\mu_{0}=-3$, $\mu_{1}=0$, $(\beta_{01},\ldots,\beta_{04})=(10,24,15,20)$, and $(\beta_{11},\beta_{12},\beta_{13})=(20,17,10)$. Additionally, the variables are specified as follows: $\epsilon_{0,i}\sim N(0,1)$, $\epsilon_{1,i}\sim N(0,9)$, $X_{i1}\sim\text{Beta}(3,4)$, $X_{i2}\sim\text{Uniform}(-2,2)$ with these two variables being independent of each other. The additional covariates $X_{i3},\dots,X_{ip}$ are first generated as in Model 1. Then we randomly select $\left\lfloor p/3\right\rfloor $ covariates from the additional covariates and multiply them by either $X_{i1}$ or $X_{i2}$ with equal probability to form the final additional covariates. The randomization variable is an additional variable taking values in $\{1,2,3,4\}$ with probabilities 0.2, 0.3, 0.3, and 0.2, respectively, and is independent of $X_{ij}$ for $j=1,\dots,p$. The simulation results for Model 2 are presented in Table (ref). The results indicate that random forests-based calibration estimators (\texttt{cal_rf}, \texttt{cal_rfnn}, \texttt{cal_rflm}, \texttt{cal_rf_g}) consistently outperform other methods across different randomization schemes and sample sizes. The AIPW-based estimators (\texttt{aipw_rf} and \texttt{aipw_nn}) tend to perform worse than the calibration estimators (\texttt{cal_rf} and \texttt{cal_nn}), especially when sample size is small ($n=500$). The empirical coverage probabilities indicate that, in most cases, the estimators yield reliable 95% confidence intervals. \begin{table}[!tbh] \caption{The comparison of the performance of different estimators under Model 2.} \resizebox{\textwidth}{!}{ \begin{threeparttable} \begin{centering} \begin{tabular}{clrrrrrrrrrrrr} \toprule \multirow{2}{*}{$n$} & \multirow{2}{*}{Estimator} & \multicolumn{4}{c}{Simple Rand.} & \multicolumn{4}{c}{Stratified Block Rand.} & \multicolumn{4}{c}{Minimization}\tabularnewline \cmidrule{3-14} & & Bias & SD & SE & CP & Bias & SD & SE & CP & Bias & SD & SE & CP\tabularnewline \midrule \multirow{11}{*}{500} & \texttt{cal_rf} & 0.01 & 2.09 & 2.30 & 0.980 & 0.19 & 2.22 & 2.30 & 0.957 & 0.28 & 2.25 & 2.29 & 0.967\tabularnewline & \texttt{cal_nn} & 0.16 & 3.09 & 3.14 & 0.957 & 0.14 & 3.19 & 3.14 & 0.940 & 0.10 & 3.08 & 3.13 & 0.960\tabularnewline & \texttt{cal_rfnn} & 0.01 & 2.14 & 2.30 & 0.977 & 0.20 & 2.25 & 2.30 & 0.960 & 0.25 & 2.26 & 2.29 & 0.957\tabularnewline & \texttt{cal_rflin} & 0.22 & 2.11 & 2.26 & 0.977 & 0.01 & 2.23 & 2.26 & 0.960 & 0.11 & 2.23 & 2.26 & 0.960\tabularnewline & \texttt{cal_rf_g} & 0.10 & 2.00 & 2.16 & 0.957 & 0.19 & 2.25 & 2.16 & 0.937 & 0.15 & 2.18 & 2.16 & 0.950\tabularnewline & \texttt{cal_nn_g} & 0.12 & 3.23 & 3.09 & 0.943 & 0.19 & 3.21 & 3.09 & 0.937 & 0.05 & 3.13 & 3.08 & 0.957\tabularnewline & \texttt{cal_lin_EL} & 0.33 & 2.68 & 2.78 & 0.950 & 0.19 & 2.73 & 2.77 & 0.960 & 0.10 & 2.69 & 2.76 & 0.963\tabularnewline & \texttt{aipw_rf} & 0.20 & 2.35 & 2.42 & 0.950 & 0.04 & 2.50 & 2.42 & 0.940 & 0.09 & 2.46 & 2.42 & 0.957\tabularnewline & \texttt{aipw_nn} & 0.87 & 5.98 & 5.99 & 0.963 & 0.47 & 6.02 & 5.92 & 0.957 & 0.25 & 5.94 & 5.87 & 0.937\tabularnewline & \texttt{aipw_lin} & 2.91 & 196.74 & 43.13 & 0.953 & 3.67 & 142.28 & 37.24 & 0.963 & 6.10 & 95.39 & 33.45 & 0.937\tabularnewline & \texttt{sdim} & 0.22 & 3.07 & 3.12 & 0.950 & 0.14 & 3.19 & 3.12 & 0.940 & 0.08 & 3.06 & 3.12 & 0.970\tabularnewline \midrule \multirow{11}{*}{1000} & \texttt{cal_rf} & 0.22 & 1.47 & 1.51 & 0.957 & 0.15 & 1.48 & 1.50 & 0.950 & 0.29 & 1.54 & 1.50 & 0.943\tabularnewline & \texttt{cal_nn} & 0.04 & 2.24 & 2.15 & 0.960 & 0.09 & 2.16 & 2.13 & 0.963 & 0.27 & 2.15 & 2.14 & 0.950\tabularnewline & \texttt{cal_rfnn} & 0.22 & 1.49 & 1.50 & 0.953 & 0.15 & 1.50 & 1.50 & 0.947 & 0.25 & 1.55 & 1.50 & 0.940\tabularnewline & \texttt{cal_rflin} & 0.08 & 1.44 & 1.47 & 0.960 & 0.08 & 1.43 & 1.47 & 0.957 & 0.15 & 1.53 & 1.46 & 0.937\tabularnewline & \texttt{cal_rf_g} & 0.17 & 1.46 & 1.45 & 0.940 & 0.12 & 1.46 & 1.45 & 0.947 & 0.25 & 1.50 & 1.45 & 0.937\tabularnewline & \texttt{cal_nn_g} & 0.02 & 2.17 & 2.02 & 0.950 & 0.11 & 2.08 & 2.00 & 0.937 & 0.30 & 2.04 & 2.02 & 0.943\tabularnewline & \texttt{cal_lin_EL} & 0.07 & 1.74 & 1.67 & 0.947 & 0.18 & 1.70 & 1.67 & 0.940 & 0.25 & 1.66 & 1.67 & 0.947\tabularnewline & \texttt{aipw_rf} & 0.07 & 1.62 & 1.59 & 0.947 & 0.06 & 1.60 & 1.59 & 0.943 & 0.26 & 1.63 & 1.59 & 0.943\tabularnewline & \texttt{aipw_nn} & 0.17 & 3.36 & 3.22 & 0.943 & 0.05 & 3.40 & 3.22 & 0.950 & 0.27 & 3.28 & 3.26 & 0.960\tabularnewline & \texttt{aipw_lin} & 0.11 & 1.75 & 1.67 & 0.950 & 0.19 & 1.73 & 1.66 & 0.947 & 0.30 & 1.67 & 1.66 & 0.947\tabularnewline & \texttt{sdim} & 0.01 & 2.25 & 2.20 & 0.950 & 0.07 & 2.27 & 2.20 & 0.933 & 0.31 & 2.20 & 2.20 & 0.960\tabularnewline \midrule \multirow{11}{*}{2000} & \texttt{cal_rf} & 0.09 & 1.00 & 1.01 & 0.943 & 0.05 & 1.04 & 1.01 & 0.937 & 0.09 & 0.97 & 1.01 & 0.953\tabularnewline & \texttt{cal_nn} & 0.03 & 1.26 & 1.27 & 0.950 & 0.13 & 1.31 & 1.27 & 0.950 & 0.05 & 1.25 & 1.27 & 0.950\tabularnewline & \texttt{cal_rfnn} & 0.07 & 1.01 & 1.01 & 0.943 & 0.07 & 1.04 & 1.01 & 0.923 & 0.07 & 0.97 & 1.01 & 0.953\tabularnewline & \texttt{cal_rflin} & 0.01 & 1.00 & 0.99 & 0.943 & 0.14 & 1.02 & 0.99 & 0.923 & 0.00 & 0.96 & 0.99 & 0.957\tabularnewline & \texttt{cal_rf_g} & 0.05 & 1.00 & 0.99 & 0.950 & 0.08 & 1.04 & 0.99 & 0.927 & 0.07 & 0.97 & 0.99 & 0.950\tabularnewline & \texttt{cal_nn_g} & 0.03 & 1.10 & 1.07 & 0.957 & 0.13 & 1.08 & 1.07 & 0.930 & 0.04 & 1.08 & 1.07 & 0.953\tabularnewline & \texttt{cal_lin_EL} & 0.01 & 1.15 & 1.13 & 0.943 & 0.12 & 1.14 & 1.13 & 0.943 & 0.03 & 1.13 & 1.13 & 0.963\tabularnewline & \texttt{aipw_rf} & 0.01 & 1.07 & 1.06 & 0.947 & 0.12 & 1.09 & 1.05 & 0.933 & 0.05 & 1.06 & 1.06 & 0.960\tabularnewline & \texttt{aipw_nn} & 0.06 & 1.38 & 1.44 & 0.957 & 0.13 & 1.57 & 1.43 & 0.923 & 0.12 & 1.49 & 1.43 & 0.940\tabularnewline & \texttt{aipw_lin} & 0.00 & 1.14 & 1.12 & 0.930 & 0.11 & 1.14 & 1.12 & 0.937 & 0.02 & 1.11 & 1.12 & 0.960\tabularnewline & \texttt{sdim} & 0.03 & 1.58 & 1.55 & 0.943 & 0.15 & 1.58 & 1.55 & 0.930 & 0.12 & 1.64 & 1.55 & 0.933\tabularnewline \bottomrule \end{tabular} \end{centering} \begin{tablenotes}[flushleft] • \textit{Abbreviations:} Rand., Randomization; SD, standard deviation, SE: standard error; CP, coverage probability. \end{tablenotes} \end{threeparttable} } \end{table} \textbf{Model 3.} Model 3 imposes non-additive and nonlinear models for the conditional mean functions $g_{0}(\bs X)$ and $g_{1}(\bs X)$. In this model, we set \begin{align*} g_{0}(\bs X_{i}) & =\mu_{0}+\beta_{01}X_{i1}X_{i2}/(X_{i1}+X_{i2}+2)+\beta_{02}X_{i1}^{2}(X_{i2}+X_{i3})\\ g_{1}(\bs X_{i}) & =\mu_{1}+\beta_{11}(X_{i2}+X_{i4})+\beta_{12}X_{i2}^{2}/\exp(X_{i1}+2), \end{align*} with $\mu_{0}=5$, $\mu_{1}=2$, $(\beta_{01},\beta_{02})=(42,83)$, and $(\beta_{11},\beta_{12})=(30,75)$. Additionally, the variables are specified as follows: $\epsilon_{0,i}\sim\text{t}(2)$, $\epsilon_{1,i}\sim3\times\text{t}(2)$, $X_{i1}\sim\text{Beta}(3,4)$, $X_{i2}\sim\text{Uniform}(-2,2)$, $X_{i3}\sim N(0,1)$, $X_{i4}\sim\text{Uniform}(0,2)$, respectively. These variables are independent of one another. The remaining variables $X_{i5},\dots,X_{ip}$ are independent of $X_{i1},\dots,X_{i4}$ and follow a multivariate normal distribution with zero mean and a symmetric Toeplitz covariance matrix where the first row is a geometric sequence with initial value 1 and common ratio 0.5. The randomization variable is an additional variable taking values in $\{1,2\}$ with probabilities 0.4, and 0.6, respectively, and is independent of $X_{ij}$ for $j=1,\dots,p$. The simulation results for Model 3 are presented in Table (ref). The results indicate that the random forests-based calibration estimators (\texttt{cal_rf}, \texttt{cal_rfnn}, \texttt{cal_rflm}, \texttt{cal_rf_g}) perform the best across different randomization methods and sample sizes. As the sample size increases, the calibration estimators converge toward better performance with lower SD, while maintaining correct empirical coverage probabilities. However, the \texttt{sdim} estimator consistently lags behind the calibration estimators in terms of SD, especially in smaller sample sizes. This finding is consistent with Theorem (ref), which demonstrates that the calibration estimator is always more efficient than the \texttt{sdim} estimator. \begin{table}[!tbh] \caption{The comparison of the performance of different estimators under Model 3.} \resizebox{\textwidth}{!}{ \begin{threeparttable} \begin{centering} \begin{tabular}{clrrrrrrrrrrrr} \toprule \multirow{2}{*}{$n$} & \multirow{2}{*}{Estimator} & \multicolumn{4}{c}{Simple Rand.} & \multicolumn{4}{c}{Stratified Block Rand.} & \multicolumn{4}{c}{Minimization}\tabularnewline \cmidrule{3-14} & & Bias & SD & SE & CP & Bias & SD & SE & CP & Bias & SD & SE & CP\tabularnewline \midrule \multirow{11}{*}{500} & \texttt{cal_rf} & 0.13 & 2.80 & 2.68 & 0.933 & 0.01 & 2.40 & 2.61 & 0.980 & 0.03 & 2.41 & 2.65 & 0.967\tabularnewline & \texttt{cal_nn} & 0.17 & 3.62 & 3.43 & 0.943 & 0.08 & 3.20 & 3.35 & 0.960 & 0.08 & 3.29 & 3.38 & 0.960\tabularnewline & \texttt{cal_rfnn} & 0.09 & 2.83 & 2.67 & 0.940 & 0.00 & 2.45 & 2.61 & 0.967 & 0.01 & 2.41 & 2.65 & 0.960\tabularnewline & \texttt{cal_rflin} & 0.08 & 2.80 & 2.63 & 0.930 & 0.01 & 2.36 & 2.56 & 0.967 & 0.02 & 2.42 & 2.61 & 0.957\tabularnewline & \texttt{cal_rf_g} & 0.16 & 2.74 & 2.58 & 0.927 & 0.01 & 2.27 & 2.52 & 0.980 & 0.11 & 2.35 & 2.56 & 0.957\tabularnewline & \texttt{cal_nn_g} & 0.08 & 3.23 & 3.16 & 0.933 & 0.08 & 2.94 & 3.08 & 0.953 & 0.08 & 2.93 & 3.12 & 0.960\tabularnewline & \texttt{cal_lin_EL} & 0.09 & 3.00 & 2.82 & 0.923 & 0.10 & 2.51 & 2.76 & 0.970 & 0.03 & 2.66 & 2.80 & 0.967\tabularnewline & \texttt{aipw_rf} & 0.04 & 3.08 & 2.84 & 0.937 & 0.16 & 2.60 & 2.78 & 0.960 & 0.16 & 2.67 & 2.82 & 0.953\tabularnewline & \texttt{aipw_nn} & 0.11 & 5.09 & 4.88 & 0.943 & 0.12 & 4.72 & 4.62 & 0.947 & 0.23 & 4.93 & 4.79 & 0.947\tabularnewline & \texttt{aipw_lin} & 0.13 & 3.16 & 3.01 & 0.937 & 0.12 & 2.79 & 2.89 & 0.930 & 0.02 & 2.83 & 2.95 & 0.977\tabularnewline & \texttt{sdim} & 0.19 & 4.31 & 4.00 & 0.930 & 0.23 & 3.77 & 3.93 & 0.950 & 0.26 & 3.87 & 3.97 & 0.960\tabularnewline \midrule \multirow{11}{*}{1000} & \texttt{cal_rf} & 0.03 & 1.87 & 1.79 & 0.947 & 0.06 & 1.82 & 1.80 & 0.940 & 0.07 & 1.66 & 1.78 & 0.957\tabularnewline & \texttt{cal_nn} & 0.09 & 2.31 & 2.18 & 0.937 & 0.02 & 2.15 & 2.20 & 0.947 & 0.02 & 2.20 & 2.18 & 0.960\tabularnewline & \texttt{cal_rfnn} & 0.06 & 1.89 & 1.79 & 0.937 & 0.08 & 1.82 & 1.79 & 0.953 & 0.08 & 1.68 & 1.78 & 0.960\tabularnewline & \texttt{cal_rflin} & 0.06 & 1.88 & 1.76 & 0.943 & 0.05 & 1.82 & 1.77 & 0.943 & 0.08 & 1.69 & 1.75 & 0.953\tabularnewline & \texttt{cal_rf_g} & 0.01 & 1.87 & 1.75 & 0.933 & 0.06 & 1.76 & 1.75 & 0.957 & 0.11 & 1.67 & 1.74 & 0.953\tabularnewline & \texttt{cal_nn_g} & 0.04 & 2.12 & 2.01 & 0.950 & 0.04 & 1.94 & 2.02 & 0.963 & 0.15 & 2.03 & 2.01 & 0.963\tabularnewline & \texttt{cal_lin_EL} & 0.02 & 1.96 & 1.82 & 0.943 & 0.05 & 1.87 & 1.82 & 0.953 & 0.04 & 1.76 & 1.81 & 0.967\tabularnewline & \texttt{aipw_rf} & 0.00 & 2.04 & 1.87 & 0.943 & 0.01 & 1.91 & 1.88 & 0.947 & 0.01 & 1.81 & 1.86 & 0.957\tabularnewline & \texttt{aipw_nn} & 0.20 & 2.98 & 2.66 & 0.920 & 0.01 & 2.97 & 2.70 & 0.950 & 0.05 & 2.76 & 2.66 & 0.950\tabularnewline & \texttt{aipw_lin} & 0.00 & 1.98 & 1.83 & 0.940 & 0.08 & 1.89 & 1.83 & 0.947 & 0.04 & 1.79 & 1.82 & 0.953\tabularnewline & \texttt{sdim} & 0.07 & 2.69 & 2.81 & 0.967 & 0.05 & 2.77 & 2.81 & 0.950 & 0.00 & 2.77 & 2.80 & 0.967\tabularnewline \midrule \multirow{11}{*}{2000} & \texttt{cal_rf} & 0.10 & 1.31 & 1.22 & 0.943 & 0.01 & 1.20 & 1.21 & 0.947 & 0.06 & 1.41 & 1.24 & 0.937\tabularnewline & \texttt{cal_nn} & 0.08 & 1.42 & 1.36 & 0.947 & 0.08 & 1.37 & 1.36 & 0.957 & 0.11 & 1.59 & 1.39 & 0.943\tabularnewline & \texttt{cal_rfnn} & 0.10 & 1.29 & 1.21 & 0.943 & 0.03 & 1.21 & 1.21 & 0.947 & 0.09 & 1.43 & 1.24 & 0.940\tabularnewline & \texttt{cal_rflin} & 0.10 & 1.29 & 1.20 & 0.940 & 0.01 & 1.20 & 1.20 & 0.943 & 0.06 & 1.41 & 1.23 & 0.940\tabularnewline & \texttt{cal_rf_g} & 0.12 & 1.28 & 1.19 & 0.950 & 0.03 & 1.16 & 1.19 & 0.950 & 0.08 & 1.40 & 1.22 & 0.943\tabularnewline & \texttt{cal_nn_g} & 0.10 & 1.34 & 1.28 & 0.940 & 0.07 & 1.30 & 1.28 & 0.953 & 0.13 & 1.55 & 1.31 & 0.940\tabularnewline & \texttt{cal_lin_EL} & 0.07 & 1.34 & 1.24 & 0.930 & 0.04 & 1.24 & 1.23 & 0.957 & 0.02 & 1.46 & 1.26 & 0.953\tabularnewline & \texttt{aipw_rf} & 0.05 & 1.36 & 1.25 & 0.930 & 0.03 & 1.26 & 1.25 & 0.950 & 0.01 & 1.46 & 1.27 & 0.940\tabularnewline & \texttt{aipw_nn} & 0.08 & 1.54 & 1.45 & 0.930 & 0.05 & 1.48 & 1.45 & 0.963 & 0.06 & 1.77 & 1.49 & 0.933\tabularnewline & \texttt{aipw_lin} & 0.07 & 1.31 & 1.23 & 0.933 & 0.03 & 1.25 & 1.23 & 0.953 & 0.04 & 1.57 & 1.26 & 0.930\tabularnewline & \texttt{sdim} & 0.07 & 2.14 & 1.98 & 0.927 & 0.05 & 2.06 & 1.98 & 0.943 & 0.01 & 2.17 & 2.00 & 0.947\tabularnewline \bottomrule \end{tabular} \end{centering} \begin{tablenotes}[flushleft] • \textit{Abbreviations:} Rand., Randomization; SD, standard deviation, SE: standard error; CP, coverage probability. \end{tablenotes} \end{threeparttable} } \end{table} Models 1-3 assume the homogeneity of the conditional mean functions $g_{0}(\bs X)$ and $g_{1}(\bs X)$ across different strata. In Appendix (ref), we also examine a model that accounts for the heterogeneity of the conditional mean functions $g_{0}(\bs X)$ and $g_{1}(\bs X)$ across different strata. Our findings show that the results are very similar to those obtained under Models 1--3. \section{Empirical application} In this section, we apply our calibration method to the experimental data from dupas2018Bankinga, which conducted covariate-adaptive randomized experiments to assess the impact of subsidized bank account access ($A_{i}$) on total savings ($Y_{i}$) for individuals across three countries: Uganda, Malawi, and Chile. dupas2018Bankinga estimated both the average treatment effects (ATEs) and the quantile treatment effects (QTEs) of the subsidy. More recently, jiang2023Regressionadjusted applied a regression-adjusted estimator to this dataset to analyze the QTE in Uganda. In this section, we aim to estimate the ATEs of subsidized access to bank accounts in Uganda and Malawi. The Uganda sample includes 2,159 observations, stratified into 41 strata based on gender, occupation, and bank branch, following a stratified block randomization design. The Malawi sample includes 2,108 observations, stratified into 78 strata based on occupation, gender, marital status, literacy, and whether the respondent was from the household or market, also following a stratified block randomization design. To avoid overly small strata, we exclude those containing fewer than six samples. This results in a final Uganda sample of 2,115 observations across 37 strata, and a Malawi sample of 1,987 observations across 67 strata. In both countries, half of the households were randomly assigned to receive the bank account subsidy, while the other half served as the control group. All monetary variables are winsorized at the ninety-ninth percentile to mitigate the influence of outliers. After the randomization and the intervention, dupas2018Bankinga conducted 3 rounds of follow-up surveys in Uganda and Malawi. In line with jiang2023Regressionadjusted, we focus on the first-round follow-up survey to assess the impact of the bank account subsidy on total savings. Following dupas2018Bankinga, we consider one baseline covariate $X$: the baseline value of total savings. When estimating the ATE of subsidized bank account access in Uganda, we leverage data from Malawi to obtain certain components of $\bs{\xi}_{n}(X)$ used in our calibration estimator, and vice versa when estimating the ATE for Malawi. For both countries, we consider the following estimators: (i) \texttt{cal_X}: this estimator is obtained by taking $\bs{\xi}_{n}(X)=X$ in ((ref)); (ii) \texttt{cal_X}$^{\beta}$: this estimator is obtained by taking $\bs{\xi}_{n}(X)=(X+1)^{\beta}$ in ((ref)); (iii) \texttt{cal_X_X}$^{\beta}$: this estimator is obtained by taking $\bs{\xi}_{n}(X)=(X,(X+1)^{\beta})^{\top}$ in ((ref)); (iv) \texttt{cal_info_X}: this estimator is obtained by taking $\bs{\xi}_{n}(X)=(X,\widehat{g}_{\text{rf}}^{\text{info}}(X))^{\top}$ in ((ref)), where $\widehat{g}_{\text{rf}}^{\text{info}}(X)$ is the random forest estimator of $\mathbb{E}\left[Y\mid X\right]$ using the data from the other country; (v) \texttt{cal_info_X}$^{\beta}$: this estimator is obtained by taking $\bs{\xi}_{n}(X)=((X+1)^{\beta},\widehat{g}_{\text{rf}}^{\text{info}}(X))^{\top}$ in ((ref)); (vi) \texttt{cal_info_X_X}$^{\beta}$: this estimator is obtained by taking $\bs{\xi}_{n}(X)=(X,(X+1)^{\beta},\widehat{g}_{\text{rf}}^{\text{info}}(X))^{\top}$ in ((ref)); (vii) \texttt{sdim}: the stratified difference-in-means estimator $\widehat{\tau}_{\mathrm{sdim}}$. We include the term $(X+1)^{\beta}$ because an approximately linear relationship is observed between $\log(Y+1)$ and $\log(X+1)$ in both countries. We fit linear regressions of the form $\log(Y+1)=\alpha+\beta\log(X+1)+\epsilon$, yielding an estimated $\beta$ of 0.481 for Uganda and 0.408 for Malawi. Figure (ref) and Table (ref) present the 95% confidence intervals and point estimates obtained using these methods. \begin{figure}[!tbh] \begin{centering} \end{centering} \begin{centering} \end{centering} \caption{The 95% confidence intervals for the ATE of the bank account subsidy on total household savings in Uganda and Malawi, respectively. Top panel: Uganda; bottom panel: Malawi. The estimators whose names contain “\texttt{X}” (or “\texttt{X}$^{\beta}$”) indicate that we add the component $X$ (or $(X+1)^{\beta}$) to $\protect\bs{\xi}(X)$ to obtain the calibration estimator. The estimators whose names contain “\texttt{info}” indicate that we incorporate information from the other country to form one component of $\protect\bs{\xi}(X)$.} \end{figure} \begin{table}[!tbh] \caption{The ATE estimates of the bank account subsidy effect on total household savings in Uganda and Malawi, with standard errors reported in parentheses.} \resizebox{\textwidth}{!}{ \begin{threeparttable} \begin{centering} \begin{tabular}{cccccccc} \toprule & \texttt{sdim} & \texttt{cal_X} & \texttt{cal_X}$^{\beta}$ & \texttt{cal_X_X}$^{\beta}$ & \texttt{cal_info_X} & \texttt{cal_info_X}$^{\beta}$ & \texttt{cal_info_X_X}$^{\beta}$\tabularnewline \midrule \midrule \multirow{2}{*}{Uganda} & 1.289 & 3.687 & 2.979 & 3.459 & 3.604 & 3.313 & 3.426\tabularnewline & (2.980) & (2.726) & (2.691) & (2.661) & (2.722) & (2.697) & (2.645)\tabularnewline \multirow{2}{*}{Malawi} & 0.452 & 1.431 & 0.870 & 0.787 & 1.484 & 0.582 & 1.031\tabularnewline & (1.766) & (1.679) & (1.702) & (1.674) & (1.690) & (1.708) & (1.654)\tabularnewline \bottomrule \end{tabular} \end{centering} \begin{tablenotes}[flushleft] • \textit{Note:} The estimators whose names contain “\texttt{X}“ (or “\texttt{X}$^{\beta}$”) indicate that we add the component $X$ (or $(X+1)^{\beta}$) to $\bs{\xi}(X)$ to obtain the calibration estimator. The estimators whose names contain “\texttt{info}” indicate that we incorporate information from the other country to form one component of $\bs{\xi}(X)$. \end{tablenotes} \end{threeparttable} } \end{table} The results presented in Figure (ref) and Table (ref) lead to two key observations. First, in line with the theoretical results (Theorem (ref)), the standard errors for the \texttt{cal_info_X_X}$^{\beta}$ estimator are the lowest among all estimators in both Uganda and Malawi. For example, the standard errors of \texttt{cal_info_X_X}$^{\beta}$ ATE estimates are 11.2% and 6.3% smaller than those of the stratified difference-in-means ATE estimates in Uganda and Malawi, respectively. Second, in both countries, all ATE estimates are statistically insignificant, suggesting that expanding access to basic bank accounts does not lead to a significant increase in total savings on average. This finding is consistent with the results of dupas2018Bankinga. \putbib