The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
72,587 characters
On regression-adjusted imputation estimators of the average treatment effect
\setlength{\abovedisplayskip}{5pt}
\setlength{\belowdisplayskip}{5pt}
\setlength{\abovedisplayshortskip}{5pt}
\setlength{\belowdisplayshortskip}{5pt}
\hypersetup{colorlinks,breaklinks,urlcolor=blue,linkcolor=blue}
\title{\LARGE On regression-adjusted imputation estimators of the average treatment effect}
\author{Zhexiao Lin\thanks{Department of Statistics, University of California, Berkeley, CA 94720, USA; e-mail: {\tt [email removed]}}~~~and~
Fang Han\thanks{Department of Statistics, University of Washington, Seattle, WA 98195, USA; e-mail: {\tt [email removed]}}
}
\date{}
\maketitle
\vspace{-1em}
\begin{abstract}
Imputing missing potential outcomes using an estimated regression function is a natural idea for estimating causal effects. In the literature, estimators that combine imputation and regression adjustments are believed to be comparable to augmented inverse probability weighting. Accordingly, people for a long time conjectured that such estimators, while avoiding directly constructing the weights, are also doubly robust \citep{imbens2004nonparametric,stuart2010matching}. Generalizing an earlier result of the authors \citep{lin2021estimation}, this paper formalizes this conjecture, showing that a large class of regression-adjusted imputation methods are indeed doubly robust for estimating the average treatment effect. In addition, they are provably semiparametrically efficient as long as both the density and regression models are correctly specified. Notable examples of imputation methods covered by our theory include kernel matching, (weighted) nearest neighbor matching, local linear matching, and (honest) random forests.
\end{abstract}
{\bf Keywords}: double robustness, kernel matching, nearest neighbor matching, random forests, double machine learning.
\section{Introduction}
The problem of estimating the average effect of a binary treatment on a scalar outcome under unconfoundedness and overlap conditions has had a long and rich history \citep{rosenbaum1983central,imbens2015causal}. While nowadays a large literature focuses on propensity score-based methods, alternatives that are based on regression \citep{heckman1997matching,heckman1998matching,heckman1998characterizing,hahn1998role,athey2016recursive,wager2018estimation} and matching \citep{rubin1973matching,abadie2006large,abadie2011bias} still receive persistent attention.
Regression and matching methods
relate causal inference to the imputation methods prevalent in the statistical missing value literature \citep{rubin2004multiple,tsiatis2006semiparametric,little2019statistical}. Indeed, as Guido Imbens and others (cf. \citet[Section IIIB]{imbens2004nonparametric} and \citet[Page 241]{abadie2006large}) have pointed out, both the regression and matching methods are intrinsically imputing the missing potential outcomes using, e.g., kernel matching, local linear matching, random forests, or the nearest neighbor matching. Accordingly, to be aligned with the missing value terminology, we call both of them the {\it imputation methods}.
Employing imputation methods alone can be either inefficient or lacking precision. This was discussions by \cite{robins1995semiparametric} in the missing value, \citet[Section IIID]{imbens2004nonparametric} and \cite{abadie2006large} in the causal inference, and \cite{cassel1976some} and \cite{sarndal2003model} in the survey literature. It stimulates a surge in combining imputation methods with different types of adjustments --- including the celebrated augmented inverse probability weighted (AIPW) estimators \citep{robins1994estimation,scharfstein1999adjusting} as well as its much more recent cousin, the double machine learning estimators \citep{chernozhukov2018double}--- partly in order to encourage more efficient and robust estimators.
This paper is interested in exploring the {\it double robustness} \citep{robins1994estimation,robins1997toward,scharfstein1999adjusting,bang2005doubly,kang2007demystifying} and {\it semiparametric efficiency} properties of the imputation methods when combined with {\it regression adjustments} for {\it correcting the bias}. While being proposed and studied in prominent works \citep{rubin1973use,abadie2011bias}, unlike its counterpart that integrates imputation with weighting --- e.g., propensity score \citep{robins1994estimation,hirano2003efficient} or covariate balancing \citep{chan2016globally,ben2021balancing} --- theoretical results on regression-adjusted imputation methods are extremely scarce. This may be partly explained by the fact that they are fully outcome model driven, and hence it was unclear which part is playing the role of propensity score weighting.
More specifically, in the literature, people have been long time conjecturing that combining imputation with regression adjustments (for the purpose of bias correction) would yield doubly robust estimators. This was made explicit in, e.g., \citet[Section IIID]{imbens2004nonparametric} that ``the benefit associated with combining methods is made explicit in the notion developed by Robins and Ritov (1997) of double robustness'' as well as \citet[Section 5]{stuart2010matching} that ``[matching and regression] have been shown to work best in combination... [t]his is similar to the idea of double robustness''. However, a mathematical formulation of double robustness for regression-adjusted imputation methods is still absent in the literature.
In addition to double robustness, statistical efficiency is vital for justifying any developed estimator. In a landmark paper, \cite{heckman1998matching} underpinned theoretical studies of (bias-uncorrected) imputation methods and showed that
imputation based on covariate kernel matching yields a semiparametrically efficient estimator. Nevertheless, \cite{heckman1998matching}'s result only focuses on estimating the average treatment effect on the treated (ATT). Later, \cite{abadie2006large,abadie2011bias} studied the limit theorems of NN matching for estimating both the ATT and the average treatment effect (ATE). However, the conveyed message therein is mixed, suggesting that NN matching-based imputation --- no matter bias correction is made or not --- is not semiparametrically efficient in estimating either the ATT or ATE. Except for the aforementioned two special cases, efficiency theory on (regression-adjusted) imputation methods is still largely lacking.
This paper aims to offer a general theory towards demystifying the efficiency and robustness properties of regression-adjusted imputation methods. For imputing the missing potential outcomes, we are concerned with a class of nonparametric regression methods called {\it linear smoothers} \citep{buja1989linear,fan2018local,wasserman2006all}, which include all the aforementioned examples (kernel matching, local linear matching, nearest neighbor matching, and random forests). Building on an earlier result of the authors that focuses on the nearest neighbor matching \citep{lin2021estimation}, the new theory shows:
\begin{itemize}
\item[(P1)] a linear smoother can implicitly give rise to a density ratio estimator;
\item[(P2)] imputation methods with regression adjustments in the form of \cite{rubin1973use} and \cite{abadie2011bias} constitute AIPW estimators;
\item[(P3)] these imputation methods are consistent as long as either the density model or the outcome model is correctly specified, and thus {\it doubly robust};
\item[(P4)] they further constitute asymptotically normal estimators of the ATE with the asymptotic variance attaining the semiparametric efficiency lower bound \citep{hahn1998role} if both the density and outcome models are correctly specified, and are thus {\it semiparametrically efficient};
\item[(P5)] the double machine learning \citep{chernozhukov2018double} versions of regression-adjusted imputations --- those that estimate the imputation function and the corrected bias via sample splitting and cross fitting --- can attain the properties in (P3) and (P4) while weakening some conditions.
\end{itemize}
Our results thus provide necessary theoretical support for using regression-adjusted imputation methods and establish them as useful alternatives to the weighting-based ones.
Notably speaking, the results of this paper are built on an earlier work of the authors \citep{lin2021estimation}, who established the double robustness and semiparametrical efficiency theory for \cite{abadie2011bias}'s NN matching-based ATE estimator by allowing the number of matches to diverge with the sample size. Their Lemma 5.1 reveals that \cite{abadie2011bias}'s bias-corrected NN matching estimator can be formulated as an AIPW one, which stimulates us to explore more cases. This leads to the general theory established in Section \ref{sec:general} and the study of more imputation methods elaborated on in Sections \ref{sec:example} and \ref{sec:RF}. Due to the richness of newly obtained results, we feel compelled to disseminate them to peers by writing a second manuscript.
\vspace{0.2cm}
\noindent {\bf Paper organization.} Section \ref{sec:prelim} introduces necessary notation, the preliminary setup, and those regression-adjusted imputation ATE estimators that will be analyzed in subsequent sections. Section \ref{sec:general} lays out our general theory, with examples provided in Sections \ref{sec:example} and \ref{sec:RF}. Specifically, Section \ref{sec:example} concerns imputation using kernel matching, weighted NN, and local linear matching while Section \ref{sec:RF} is focused on imputing the missing potential outcomes using random forests.
\section{Preliminary}\label{sec:prelim}
In the following, for any integers $n,d\ge 1$, we write $\llbracket n\rrbracket:= \{1,2,\ldots,n\}$, and $\bR^d$ to represent the $d$-dimensional real space. A set consisting of distinct elements $x_1,\dots,x_n$ is written as either $\{x_1,\dots,x_n\}$ or $\{x_i\}_{i=1}^{n}$, and the corresponding sequence is denoted by $[x_1,\dots,x_n]$ or $[x_i]_{i=1}^{n}$.
Consider $n$ observations, categorized to two groups, the treated and control, separately with $D_1,\ldots,D_n$ indexing the treatment statuses. More specifically, for each unit $i \in \llbracket n\rrbracket$, we observe $D_i=1$ if in the treated group and $D_i=0$ if in the control group. Let $n_0:=\sum_{i=1}^n (1-D_i)$ and $n_1:=\sum_{i=1}^n D_i$ be the numbers of control and treated units, respectively. Adopting the Neyman-Rubin potential outcome framework \citep{neyman1923applications,rubin1974estimating}, the unit $i$ has two potential outcomes, $Y_i(1)$ and $Y_i(0)$, but we observe only one of them:
\[
Y_i = \begin{cases}
Y_i(0), & \mbox{ if } D_i=0,\\
Y_i(1), & \mbox{ if } D_i=1.
\end{cases}
\]
Let $X_i$ represent the pretreatment covariates of the $i$-th unit.
The data we observe are $[(X_i,D_i,Y_i)]_{i=1}^n$, which are assumed to be independently drawn from the triple $(X,D,Y)$, where $D\in\{0,1\}$ is a binary variable, $X \in \bR^d$, and $Y \in \bR$. Our goal of interest is to estimate the following population ATE,
\begin{align*}
\tau := {\mathrm E}\Big[Y_i(1)-Y_i(0)\Big],
\end{align*}
based on $[(X_i, D_i, Y_i)]_{i=1}^n$.
As stated in the introduction section, this paper is interested in studying the imputation-based ATE estimators. To this end, we consider imputing the missing potential outcomes by regressing the data points in the opposite group against it:
\begin{align*}
\hat{Y}_i^{\rm imp}(0) := \begin{cases}
Y_i, & \mbox{ if } D_i=0,\\
\displaystyle\sum_{j:D_j=0} w_{i\leftarrow j} Y_j, & \mbox{ if } D_i=1,
\end{cases}
\end{align*}
and
\begin{align*}
\hat{Y}_i^{\rm imp}(1) := \begin{cases}
\displaystyle \sum_{j:D_j=1} w_{i\leftarrow j} Y_j, & \mbox{ if } D_i=0,\\
Y_i, & \mbox{ if } D_i=1.
\end{cases}
\end{align*}
Here the $[w_{i\leftarrow j}]_{i,j}$ constitutes the {\it smoothing matrix}, where each entry $w_{i\leftarrow j}$ --- called the {smoothing parameter} --- is learnt from the covariates $X_i$ and those $X_j$'s in the opposite group, i.e., those with $D_j=1-D_i$.
Nonparametric regressors taking the above form are called the {\it linear smoothers} \citep{buja1989linear}. Note that all imputation methods considered in Sections \ref{sec:example} and \ref{sec:RF}, including the kernel regression and local linear regression estimators \citep{heckman1997matching,heckman1998characterizing,heckman1998matching}, the (weighted) NN regression \citep{abadie2006large,abadie2011bias,lin2021estimation}, and the (honest) random forests \citep{athey2016recursive,wager2018estimation,athey2019estimating}, admit such a form.
Unfortunately, imputing the missing potential outcomes alone is often not sufficient for attaining efficiency or even merely root-$n$ consistency. To remedy it, we are interested in correcting the bias via regression adjustments as proposed in \cite{rubin1973use} and \cite{abadie2011bias}. In detail, let's write
\[
\hat{\mu}_0(x)~~~{\rm and}~~~\hat{\mu}_1(x)
\]
to represent the mappings from $\bR^d$ to $\bR$ that estimate the conditional means of the outcomes
\[
\mu_0(x) := {\mathrm E} [Y \,|\, X=x,D=0]~~ {\rm and}~~ \mu_1(x) := {\mathrm E} [Y \,|\, X=x,D=1],
\]
respectively. Of note, in the literature, $\hat{\mu}_0(x)$ and $\hat{\mu}_1(x)$ may differ from the regression imputation methods used in calculating $\hat{Y}_i^{\rm imp}(0)$'s and $\hat{Y}_i^{\rm imp}(1)$'s. For example, \cite{abadie2011bias} used NN regression to impute the missing potential outcomes, but series regressions to correct the bias.
We are then ready to define the regression-adjusted imputed values as
\begin{align*}
\hat{Y}_i(0) := \begin{cases}
Y_i, & \mbox{ if } D_i=0,\\
\displaystyle \sum_{j:D_j=0} w_{i\leftarrow j} (Y_j + \hat{\mu}_0(X_i) - \hat{\mu}_0(X_j)), & \mbox{ if } D_i=1,
\end{cases}
\end{align*}
and
\begin{align*}
\hat{Y}_i(1) := \begin{cases}
\displaystyle \sum_{j:D_j=1} w_{i\leftarrow j} (Y_j + \hat{\mu}_1(X_i) - \hat{\mu}_1(X_j)), & \mbox{ if } D_i=0,\\
Y_i, & \mbox{ if } D_i=1.
\end{cases}
\end{align*}
The according regression-adjusted imputation-based ATE estimator is
\begin{align*}
\hat\tau_w := \frac{1}{n} \sum_{i=1}^n \Big[\hat{Y}_i(1) -\hat{Y}_i(0)\Big].
\end{align*}
The estimator $\hat\tau_w$ has the appealing property of being fully outcome model driven, i.e., both the imputation and the bias correction steps are regression-based. It is conceptually easy to parse. The first goal of this paper is to show that $\hat\tau_w$, while avoiding directly modeling the propensity score, can be formulated as an AIPW one, and the regression imputation is intrinsically estimating the propensity score. The second goal of this paper is to establish a general theory, formulating conditions under which $\hat\tau_w$ is doubly robust and semiparametrically efficient. Examples covered by our general theory shall occupy the rest two sections of this paper.
\section{The general theory}\label{sec:general}
This section lays out the general theory on the regression-adjusted imputation estimator $\hat\tau_w$. Recall the conditional mean estimators $\hat\mu_0$ and $\hat\mu_1$ introduced in the last section. Let the residuals from fitting the outcome models be
\[
\hat{R}_i := Y_i - \hat{\mu}_{D_i}(X_i), ~~i\in\llbracket n\rrbracket,
\]
and the estimator based on the outcome models be
\[
\hat{\tau}^{\rm reg}:= n^{-1} \sum_{i=1}^n \Big[\hat{\mu}_1(X_i) - \hat{\mu}_0(X_i)\Big].
\]
\subsection{A key lemma}
Results in Section \ref{sec:general} are all built on the following key lemma, which gives an AIPW formulation of the ATE estimator $\hat\tau_w$.
\begin{lemma}\label{lemma:mbc} The regression-adjusted imputation estimator $\hat\tau_w$ can be rewritten as
\begin{align}\label{eq:mbc}
\hat\tau_w = \hat{\tau}^{\rm reg} + \frac{1}{n} \sum_{i=1,D_i=1}^n \Big(1 + \sum_{j:D_j=1-D_i} w_{j\leftarrow i}\Big) \hat{R}_i- \frac{1}{n} \sum_{i=1,D_i=0}^n \Big(1 + \sum_{j:D_j=1-D_i} w_{j\leftarrow i}\Big) \hat{R}_i \notag\\
+ \frac{1}{n} \sum_{i=1}^n (2D_i-1) \Big(1 - \sum_{j:D_j=1-D_i} w_{i\leftarrow j} \Big) \hat{\mu}_{1-D_i}(X_i).
\end{align}
\end{lemma}
The sum of the first three terms in \eqref{eq:mbc} has the same form as an AIPW estimator that was studied in \cite{scharfstein1999adjusting} and \cite{bang2005doubly}, among many others. The last is an additional bias term that was induced by those unnormalized $w_{i\leftarrow j}$'s such that
\[
\sum_{j:D_j=1-D_i} w_{i\leftarrow j} \ne 1.
\]
Accordingly, Equation \eqref{eq:mbc} favors a normalized smoothing matrix such that $\sum_{j:D_j=1-D_i} w_{i\leftarrow j}$ adds up to 1. This is an observation interestingly related to the classic arguments in nonparametric regressions; cf. \citet[Section 3.4]{fan2018local} and \citet[Remark 5.23]{wasserman2006all}.
Note that the relation between regression-adjusted imputation and AIPW estimators was for the first time disclosed in \cite{lin2021estimation}, stated as Lemma 5.1 therein and with a focus on NN regression-based imputation. Lemma \ref{lemma:mbc}, on the other hand, delivers the general form that applies to an arbitrary linear smoother.
\subsection{Double robustness}
For presenting the general theory, let us first introduce some additional notation. In the sequel, for any two real sequences $\{a_n\}$ and $\{b_n\}$, we write $a_n = O(b_n)$ if $\lvert a_n \rvert / \lvert b_n \rvert $ is bounded and $a_n = o(b_n)$ if $\lvert a_n \rvert / \lvert b_n \rvert \to 0$. We use $\stackrel{\sf d}{\longrightarrow}$ and $\stackrel{\sf p}{\longrightarrow}$ to denote convergence in distribution and in probability, respectively. For any sequence of random variables $[X_n]$, write $X_n = o_{\mathrm P}(1)$ if $X_n \stackrel{\sf p}{\longrightarrow} 0$ and $X_n = O_{\mathrm P}(1)$ if $X_n$ is bounded in probability. For any vector $x$, we use $\lVert x \rVert$ to denote its Euclidean norm. For any $0<p\le \infty$ and function $f$, let $\lVert f(Z) \rVert_p$, or simply $\lVert f \rVert_p$ if no confusion is possible, to represent $(\int \lvert f(\omega) \rvert^p {\mathrm d} {\mathrm P}_Z(\omega))^{1/p}$, where ${\mathrm P}_Z$ represents the law of a certain random variable $Z$.
In the following, let $U_\omega := Y(\omega) - \mu_{\omega}(X)$ for $\omega \in \{0,1\}$ be the residuals of $Y(0)$ and $Y(1)$ projected on $X$ and let $\cS$ be the support of $X$.
The first set of assumptions concerns the data generating distribution.
\begin{assumption} \phantomsection \label{asp:dr}
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item\label{asp:dr-1} For almost all $x \in \cS$, $D$ is independent of $(Y(0),Y(1))$ conditional on $X=x$, and there exists some constant $\eta > 0$ such that $\eta < {\mathrm P}(D=1 \,|\, X=x) < 1-\eta$.
\item $[(X_i,D_i,Y_i)]_{i=1}^n$ are independent and identically distributed (i.i.d.) following the joint distribution of $(X,D,Y)$.
\item ${\mathrm E} [U^2_\omega \,|\, X=x] $ is uniformly bounded for almost all $x \in \cS$ and $\omega \in \{0,1\}$.
\item ${\mathrm E} [\mu^2_\omega(X)]$ is bounded for $\omega \in \{0,1\}$.
\end{enumerate}
\end{assumption}
Assumption~\ref{asp:dr}\ref{asp:dr-1} is the unconfoundedness and overlap assumptions commonly assumed in the literature. In particular, $e(x):={\mathrm P}(D=1\,|\, X=x)$ is the propensity score \citep{rosenbaum1983central}. The rest conditions in Assumption \ref{asp:dr} constitute standard i.i.d. assumptions and the moment assumptions on the residuals.
The next set of assumptions concerns the smoothing matrix used in the imputation step.
\begin{assumption} \phantomsection \label{asp:weight}
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item\label{asp:weight-1} Let $\pi: \llbracket n\rrbracket \to \llbracket n\rrbracket$ be any permutation. For samples $[(X_i,D_i,Y_i)]_{i=1}^n$ given, let $[w_{i\leftarrow j}]_{D_i+D_j=1}$ be the weights constructed by $[(X_i,D_i,Y_i)]_{i=1}^n$, and $[w^\pi_{i\leftarrow j}]_{D_i+D_j=1}$ be the weights constructed by $[(X_{\pi(i)},D_{\pi(i)},Y_{\pi(i)})]_{i=1}^n$. Then for any $i,j \in \llbracket n\rrbracket$ such that $D_i+D_j=1$ and any permutation $\pi$, we have $w_{i\leftarrow j} = w^\pi_{\pi(i)\leftarrow \pi(j)}$.
\item \label{asp:weight-2} The weights satisfy
\begin{align*}
\lim_{n \to \infty }{\mathrm E} \Big[
\sum_{j:D_j=1-D_1} w_{1\leftarrow j} - 1
\Big]^2 = 0.
\end{align*}
\end{enumerate}
\end{assumption}
Assumption~\ref{asp:weight} is to our knowledge new and is added for aiding the general theory to be presented later. There Assumption~\ref{asp:weight}\ref{asp:weight-1} ensures that the regression smoothing matrix is invariant to the feeding order of sample points, and Assumption~\ref{asp:weight}\ref{asp:weight-2} ensures that the bias term in Lemma~\ref{lemma:mbc} is asymptotically ignorable, which will be automatically satisfied if the smoother preserves the constant curve \citep[Remark 5.23]{wasserman2006all}.
The next set of assumptions quantifies estimation accuracy of the ``density models''.
\begin{assumption} \phantomsection \label{asp:dr1}
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item\label{asp:dr1,o} For $\omega \in \{0,1\}$, there exists a deterministic function $\bar{\mu}_\omega(\cdot):\bR^d \to \bR$ such that ${\mathrm E} [\bar{\mu}^2_\omega(X)]$ is bounded and the estimator $\hat{\mu}_\omega(x)$ satisfies
\[
\lVert \hat{\mu}_\omega - \bar{\mu}_\omega \rVert_\infty = o_{\mathrm P}(1).
\]
\item\label{asp:dr1,w} The weights satisfy
\begin{align*}
\lim_{n \to \infty }{\mathrm E} \Big[
\sum_{j:D_j=1-D_1} w_{j\leftarrow 1} - \Big(D_1 \frac{1-e(X_1)}{e(X_1)} + (1-D_1) \frac{e(X_1)}{1-e(X_1)} \Big)
\Big]^2 = 0.
\end{align*}
\end{enumerate}
\end{assumption}
Assumption \ref{asp:dr1} allows for outcome model misspecification. Here Assumption \ref{asp:dr1}\ref{asp:dr1,o} is a regression misspecification assumption that is Assumption 5.3 in \cite{lin2021estimation}. Assumption \ref{asp:dr1}\ref{asp:dr1,w} is the key assumption that relates regression imputation/linear smoothers to the estimation of density ratios, in the form of $(1-e(x))/e(x)$ and its inverse; in Sections \ref{sec:example} and \ref{sec:RF} we will verify its validity for a variety of regression imputation methods.
In parallel to Assumption \ref{asp:dr1}, the following conditions quantify estimation accuracy of the ``outcome models''.
\begin{assumption} \phantomsection \label{asp:dr2}
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item\label{asp:dr2,o} For $\omega \in \{0,1\}$, the estimator $\hat{\mu}_\omega(x)$ satisfies
\[
\lVert \hat{\mu}_\omega - \mu_\omega \rVert_\infty = o_{\mathrm P}(1).
\]
\item\label{asp:dr2,w} The weights $[w_{1\leftarrow j}]_{D_j=1-D_1}$ are constructed by $[(X_i,D_i)]_{i=1}^n$ only without using the outcome information $[Y_i]_{i=1}^n$.
\item\label{asp:dr2,w2} The weights satisfy
\begin{align*}
{\mathrm E} \Big[ \Big\lvert \sum_{j:D_j=1-D_1} w_{j\leftarrow 1} \Big\rvert \Big] = O(1).
\end{align*}
\end{enumerate}
\end{assumption}
Assumption \ref{asp:dr2} allows for density model misspecification. Here Assumption \ref{asp:dr2}\ref{asp:dr2,o} is Assumption 5.4 in \cite{lin2021estimation}; \cite{chen2015optimal} and \cite{chen2018optimal} verified such conditions for various nonparametric regressors. Assumption~\ref{asp:dr2}\ref{asp:dr2,w} ensures that the responses are not used in the construction of weights, and is satisfied by all examples to be introduced in Sections \ref{sec:example} and \ref{sec:RF}. This assumption is also related to the sample splitting procedures used in the context of double machine learning \citep{chernozhukov2018double} and honest random forests \citep{wager2018estimation}, shown to help avoid overfitting. Given Assumptions \ref{asp:dr} and \ref{asp:weight}, Assumption~\ref{asp:dr2}\ref{asp:dr2,w2} holds automatically as long as all the weights $w_{i\leftarrow j}$'s are nonnegative, or when Assumption~\ref{asp:dr1}\ref{asp:dr1,w} holds.
We would also like to highlight that Assumption~\ref{asp:dr2}\ref{asp:dr2,w2} is only needed for proving double robustness properties.
With the above assumptions, we are now ready to formalize the double robustness property of the regression-adjusted imputation estimator $\hat\tau_w$.
\begin{theorem}[Double robustness of $\hat\tau_w$] \label{thm:dr}
Suppose Assumptions \ref{asp:dr} and \ref{asp:weight} hold, and either Assumption~\ref{asp:dr1} or Assumption~\ref{asp:dr2} is true. We then have
\begin{align*}
\hat\tau_w - \tau \stackrel{\sf p}{\longrightarrow} 0.
\end{align*}
\end{theorem}
Theorem \ref{thm:dr} unveils an interesting phenomenon that, although regression-adjusted imputation methods are {\it fully outcome model driven}, they are doubly robust and an intrinsic statistic coming from imputation captures the role of the propensity score; cf. Assumption \ref{asp:dr1}\ref{asp:dr1,w}. To the authors' knowledge, both the missing value and causal inference literature is largely silent about this phenomena. The most related result to Theorem \ref{thm:dr} resides in simple parametric models.
In detail, the fact that ordinary least square (OLS) is intrinsically a weighted estimator is very well known; cf. \cite{angrist2009mostly} and \cite{imbens2015matching}. In two very interesting papers, \citet[Section 3]{robins2007comment} and \cite{kline2011oaxaca} showed that OLS is also able to offer double robustness guarantee for estimating either a population mean with incomplete data or the ATT. This was developed more sophistically in a recent work of \cite{chattopadhyay2021implied} and other interesting research along this line includes \cite{guo2021generalized} and \cite{cohen2020no}. In the high level, they all bear a similar flavor to Theorem \ref{thm:dr} that a regression/imputation approach, without designing a set of weights (propensity score-based or not) on purpose, automatically satisfies the double robustness property. The difference with ours, on the other hand, is self-explanatory.
\subsection{Semiparametric efficiency}
This section establishes the semiparametric efficiency theory of $\hat\tau_w$. To this end, it appears that we have to put more assumptions on the moments of $U_\omega$, the regression adjustments $\hat\mu_w(\cdot)$, and the smoothing parameters $w_{i\leftarrow j}$'s.
\begin{assumption} \phantomsection \label{asp:se1}
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item ${\mathrm E} [U^2_\omega \,|\, X=x]$ is uniformly bounded away from zero for almost all $x \in \cS$ and $\omega \in \{0,1\}$.
\item There exists some constant $\kappa>0$ such that ${\mathrm E} [\lvert U_\omega \rvert ^{2+\kappa} \,|\, X=x]$ is uniformly bounded for almost all $x \in \cS$ and $\omega \in \{0,1\}$.
\end{enumerate}
\end{assumption}
\begin{assumption} \phantomsection \label{asp:se2}
There exists a positive integer $k$ such that
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item\label{asp:se2,o1} $\max_{t \in \Lambda_{k}} \lVert \partial^t \mu_{\omega} \rVert_\infty$ is bounded, where for any positive integer $k$, $\Lambda_k$ is the set of all $d$-dimensional vectors of nonnegative integers $t=(t_1,\ldots,t_d)$ such that $\sum_{i=1}^d t_i = k$;
\item\label{asp:se2,o2} For $\omega \in \{0,1\}$, the estimator $\hat{\mu}_\omega(x)$ satisfies
\[
\max_{t \in \Lambda_{k}} \lVert \partial^t \hat{\mu}_{\omega} \rVert_\infty = O_{\mathrm P}(1)~~~{\rm and}~~~
\max_{t \in \Lambda_\ell} \lVert \partial^t \hat{\mu}_{\omega} - \partial^t \mu_{\omega} \rVert_\infty = O_{\mathrm P}(n^{-\gamma_\ell}) ~~\mbox{\rm for all}~~ \ell \in \llbracket k-1\rrbracket,
\]
with some constants $\gamma_\ell$'s for $\ell=1,2,\ldots,k-1$;
\item\label{asp:se2,d} The discrepancy satisfies
\begin{align*}
& {\mathrm E} \Big[ \sum_{j:D_j=1-D_1} \lvert w_{1\leftarrow j} \rvert \cdot \lVert X_j - X_1 \rVert^k \Big] = o(n^{-1/2}),\\
& {\mathrm E} \Big[ \sum_{j:D_j=1-D_1} \lvert w_{1\leftarrow j} \rvert \cdot \lVert X_j - X_1 \rVert^\ell \Big] = o(n^{-1/2+ \gamma_\ell} ) ~\mbox{\rm for all}~ \ell \in \llbracket k-1\rrbracket;
\end{align*}
\item\label{asp:se2,w1} The weights satisfy
\begin{align*}
{\mathrm E} \Big[
\sum_{j:D_j=1-D_1} w_{1\leftarrow j} - 1
\Big]^2 = o(n^{-1}).
\end{align*}
\end{enumerate}
\end{assumption}
Assumptions~\ref{asp:se1} and \ref{asp:se2}\ref{asp:se2,o1}-\ref{asp:se2,o2} are Assumptions 5.6 and 5.7 in \cite{lin2021estimation}; check \cite{abadie2011bias} and \cite{chen2018optimal} for results on verifying these requirements. Assumption~\ref{asp:se2}\ref{asp:se2,d} assumes that the linear smoother used in imputing the missing values is a local method, i.e., it will put larger values on the closer ones and smaller values on the farther ones. Lastly, Assumption \ref{asp:se2}\ref{asp:se2,w1}, as a counterpart of Assumption \ref{asp:weight}\ref{asp:weight-2}, requires the bias term in \eqref{eq:mbc} to be root-$n$ ignorable.
We then introduce the semiparametric efficiency lower bound for estimating the ATE \citep{hahn1998role},
\begin{align}\label{eq:hahn}
\sigma^2:= {\mathrm E} \Big[\mu_1(X) - \mu_0(X) + \frac{D(Y-\mu_1(X))}{e(X)} - \frac{(1-D)(Y-\mu_0(X))}{1-e(X)} - \tau \Big]^2.
\end{align}
The following theorem then shows that the asymptotic variance of $\hat\tau_w$ can attain $\sigma^2$.
\begin{theorem}[Semiparametric efficiency of $\hat\tau_w$]\label{thm:mbc}
Suppose Assumptions~\ref{asp:dr}-\ref{asp:se2} hold. We then have
\begin{align*}
\sqrt{n} (\hat\tau_w - \tau) \stackrel{\sf d}{\longrightarrow} N(0,\sigma^2).
\end{align*}
In addition, the variance estimator
\begin{align*}
\hat{\sigma}^2:= \frac{1}{n} \sum_{i=1}^n \Big[\hat{\mu}_1(X_i) - \hat{\mu}_0(X_i) + (2D_i-1)\Big(1 + \sum_{j:D_j=1-D_i} w_{j\leftarrow i}\Big) \hat{R}_i - \hat\tau_w \Big]^2
\end{align*}
is a consistent estimator of $\sigma^2$ in \eqref{eq:hahn}.
\end{theorem}
\subsection{Double machine learning}
Assumptions \ref{asp:se2}\ref{asp:se2,o1}-\ref{asp:se2,o2} are arguably strong regularity conditions for the outcome model $\mu_\omega(\cdot)$. Partly in order to alleviate such requirements, \cite{chernozhukov2018double} introduced the idea of double machine learning via sample splitting and cross fitting. Similar ideas have also been studied in nonparametric statistics; cf. \citet{bickel1982adaptive}, \cite{efromovich1996nonparametric}, and \cite{zheng2010asymptotic}. In the following, let's introduce $\tilde{\tau}_{w,N}$ as a counterpart of $\hat{\tau}_w$ based on \citet[Definition 3.1]{chernozhukov2018double}.
In detail, let $N \ge 2$ represent a fixed number of partitions. For presentation simplicity and also without much loss of generality, assume $n$ to be divisible by $N$. Let $[I_k]_{k=1}^N$ be an $N$-fold random partition of $\llbracket n\rrbracket$, with each of size equal to $n' = n/N$. For each $k \in \llbracket N\rrbracket$ and $\omega \in \{0,1\}$, construct $\hat{\mu}_{\omega,k}(\cdot)$ using data $[(X_i,D_i,Y_i)]_{i=1,i \notin I_k}^n$. Similarly, for regression imputation, we impute each unit's value by regressing it against all units in the opposite group outside the $k$-th fold. More specifically, we calculate the smoothing matrix entries as follows: for any $i,j \in \llbracket n\rrbracket$ with $D_i+D_j=1,i \in I_k,j \notin I_k$, let $w_{j\leftarrow i,k}$ be the weights constructed using data $(X_i,D_i,Y_i) \cup [(X_j,D_j,Y_j)]_{j=1,j \notin I_k}^n$.
We are then ready to define the double machine learning version of $\hat\tau_w$ as follows:
\begin{align*}
\widecheck{\tau}_{w,k} :=& \frac{1}{n'} \sum_{i=1,i \in I_k}^n \Big[\hat{\mu}_{1,k}(X_i) - \hat{\mu}_{0,k}(X_i)\Big] \\
&+ \frac{1}{n'} \sum_{i=1,i \in I_k}^n (2D_i-1)\Big(1 + \sum_{j:D_j=1-D_i,j \notin I_k} w_{j\leftarrow i,k}\Big) \Big(Y_i - \hat{\mu}_{D_i,k}(X_i)\Big)
\end{align*}
and
\begin{align*}
\tilde{\tau}_{w,N} := \frac{1}{N} \sum_{k=1}^N \widecheck{\tau}_{w,k}.
\end{align*}
For establishing the efficiency theory of $\tilde{\tau}_{w,N} $, the following two sets of assumptions are needed.
\begin{assumption} \phantomsection \label{asp:dml1}
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item ${\mathrm E} [U^2_\omega]$ is bounded away from zero for $\omega \in \{0,1\}$.
\item There exists some constant $\kappa>0$ such that ${\mathrm E} [\lvert Y \rvert^{2+\kappa}]$ is bounded.
\end{enumerate}
\end{assumption}
\begin{assumption}\label{asp:dml2}
There exist two positive integers $1 \le p_1,p_2 \le \infty$ with $p_1^{-1} + p_2^{-1} = 1$, two positive real-valued sequences $[r_1] = [r_1]_n, [r_2] = [r_2]_n$ with $r_1r_2 = o(n^{-1/2})$ such that
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item\label{asp:dml2,o} for $\omega \in \{0,1\}$, the estimator $\hat{\mu}_\omega(x)$ satisfies
\[
\lVert \tilde{\mu}_\omega - \mu_\omega \rVert_{p_1} = O_{\mathrm P}(r_1);
\]
\item\label{asp:dml2,w} the weights $[w_{i\leftarrow j}]_{D_i+D_j=1}$ satisfy
\begin{align*}
& \Big\{{\mathrm E} \Big[ \Big\lvert \sum_{j:D_j=1-D_1} w_{j\leftarrow 1} - \Big(D_1 \frac{1-e(X_1)}{e(X_1)} + (1-D_1) \frac{e(X_1)}{1-e(X_1)} \Big) \Big\rvert^{p_2} \Big] \Big\}^{1/p_2} = O(r_2),\\
{\rm and}\quad & {\mathrm E} \Big[ \sum_{j:D_j=1-D_1} w_{j\leftarrow 1} \Big]^\kappa = O(1), ~~ {\rm for~any~} \kappa>0.
\end{align*}
\end{enumerate}
\end{assumption}
We are now ready to introduce the general theory on the double machine learning-based regression-adjusted imputation estimators.
\begin{theorem} \phantomsection \label{thm:dml}
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item\label{thm:dml1} (Double robustness of $\tilde{\tau}_{w,N}$) Under the same conditions as those in Theorem~\ref{thm:dr}, we have
\begin{align*}
\tilde{\tau}_{w,N} - \tau \stackrel{\sf p}{\longrightarrow} 0.
\end{align*}
\item\label{thm:dml2} (Semiparametric efficiency of $\tilde{\tau}_{w,N}$) Under Assumptions~\ref{asp:dr}-\ref{asp:dr2} and \ref{asp:dml1}-\ref{asp:dml2}, we have
\begin{align*}
\sqrt{n} (\tilde{\tau}_{w,N} - \tau) \stackrel{\sf d}{\longrightarrow} N(0,\sigma^2).
\end{align*}
In addition, the variance estimator
\begin{align*}
\hat{\sigma}^2_N:= \frac{1}{n} \sum_{i=1}^n \Big[\hat{\mu}_1(X_i) - \hat{\mu}_0(X_i) + (2D_i-1)\Big(1 + \sum_{j:D_j=1-D_i} w_{j\leftarrow i}\Big) \hat{R}_i - \tilde{\tau}_{w,N} \Big]^2
\end{align*}
is a consistent estimator for $\sigma^2$ in \eqref{eq:hahn}.
\end{enumerate}
\end{theorem}
\section{Examples}\label{sec:example}
This section aims to provide examples so to put the general theory introduced in Section \ref{sec:general} on a solid ground. In the sequel, write $\ind(\cdot)$ to represent the indicator function and $a_n \asymp b_n$ if both $a_n = O(b_n)$ and $b_n = O(a_n)$ holds. For any matrix $A$, we use $\lvert A \rvert$ and $\lVert A \rVert_2$ to denote its determinant and spectral norm. For any set $\cS$, let ${\rm diam}(\cS):=\sup_{x,y\in \cS}\lVert x-y \rVert$ be its diameter.
\subsection{Kernel matching}
We first consider the kernel matching that has been advocated in various settings \citep{heckman1997matching, heckman1998characterizing, heckman1998matching, frolich2004finite, frolich2005matching, huber2013performance}. It leverages the local constant regression (Nadaraya–Watson estimator) to impute the missing values \citep{nadaraya1964estimating, watson1964smooth}.
More specifically, let $H = H_n \in \bR^{d \times d}$ be the {\it bandwidth matrix} and $K(\cdot): \bR^d \to \bR$ be the {\it multivariate kernel function} on $\bR^d$. For any $x \in \bR^d$, define
\[
K_H(x) := \lvert H \rvert^{-1/2} K(H^{-1/2}x).
\]
For any $i,j \in \llbracket n\rrbracket$ such that $D_i+D_j=1$, one can then verify that the weight $w_{i\leftarrow j}$ corresponding to kernel matching is
\begin{align*}
w_{i\leftarrow j}:= \frac{K_H(X_i-X_j)}{\sum_{k:D_k=1-D_i} K_H(X_i-X_k)}.
\end{align*}
Denote the corresponding kernel matching estimator using the above smoothing matrix as well as the double machine learning version of it by
\[
\hat{\tau}_{\rm K}~~{\rm and}~~ \tilde{\tau}_{{\rm K},N}.
\]
Assumptions in Section \ref{sec:general} can then be shown to hold under the following sufficient conditions.
\begin{assumption}\phantomsection \label{asp:kernel} Assume that (i)
$H$ is symmetric and positive definite, and (ii) $K$ constitutes a multivariate symmetric density function.
\end{assumption}
\begin{assumption} \phantomsection \label{asp:kernel,dr}
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item The density of $X$ is bounded and bounded away from zero. The densities of $X \,|\, D=1$ and $X \,|\, D=0$ are continuous almost everywhere.
\item $K$ is bounded with a compact support such that $\lVert H^{1/2} \rVert_2 \to 0$ and $n \lvert H^{1/2} \rvert \to \infty$.
\end{enumerate}
\end{assumption}
\begin{assumption} \phantomsection \label{asp:kernel,se} There exists a positive integer $k$ such that
\begin{itemize}
\item[(i)] Assumptions~\ref{asp:se2}\ref{asp:se2,o1},\ref{asp:se2,o2} hold;
\item[(ii)] we further have $\lVert H^{1/2} \rVert_2^k = o(n^{-1/2})$ and $\lVert H^{1/2} \rVert_2^\ell = o(n^{-1/2+ \gamma_\ell})$ for all $\ell \in \llbracket k-1\rrbracket$.
\end{itemize}
\end{assumption}
\begin{assumption}[double machine learning] \phantomsection \label{asp:kernel,dml}
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item The densities of $X \,|\, D=1$ and $X \,|\, D=0$ are Lipchitz on $\cS$. The diameter and the surface area (Hausdorff measure, \citet[Section 3.3]{evans2018measure}) of $\cS$ are bounded.
\item There exist two positive real-valued sequences $[r_1] = [r_1]_n, [r_2] = [r_2]_n$ with $r_1r_2 = o(n^{-1/2})$ such that
\[
\lVert \tilde{\mu}_\omega - \mu_\omega \rVert_\infty = O_{\mathrm P}(r_1)~~ {\rm for} ~~\omega \in \{0,1\},
\]
and $(n \lvert H^{1/2} \rvert)^{-1/2} + \lVert H^{1/2} \rVert_2 = O(r_2)$.
\end{enumerate}
\end{assumption}
Assumption~\ref{asp:kernel,dr} is standard for establishing consistency of the Nadaraya-Watson estimator. Assumption~\ref{asp:kernel,se} ensures that the discrepancy level in Assumption~\ref{asp:se2} is small. The regularity condition on the support and the smoothness condition on the density function are standard in nonparametric statistics \cite[Section 2]{MR2724359}.
\begin{remark}
A specific common choice of $H$ is $h_n^2 I_d$, where $I_d$ is the $d$-dimensional identity matrix. The bandwidth selection condition in Assumption~\ref{asp:kernel,dr} then reduces to
\[
h_n \to 0\quad {\rm and}\quad nh_n^d \to \infty,
\]
and Assumption~\ref{asp:kernel,se} reduces to
\[
h_n/n^{-1/(2k)} \to 0\quad {\rm and}\quad h_n/n^{(-1/2+\gamma_\ell)/\ell} \to 0 ~~{\rm for}~~ \ell \in \llbracket k-1\rrbracket,
\]
suggesting that the bandwidth cannot be too large; this echos the NN matching case where the number of NNs incorporated also has to be controlled (Theorem 5.2 in \cite{lin2021estimation}). The convergence rate in Assumption~\ref{asp:kernel,dml} reduces to $(nh^d)^{-1/2}+h$, and is the minimax rate of the density estimation over Lipchitz class $n^{-1/(2+d)}$ \cite[Section 2]{MR2724359} by taking $h_n \asymp n^{-1/(2+d)}$.
\end{remark}
The following theorem then verifies the general conditions presented in Section \ref{sec:general} when kernel matching is used for imputing the missing potential outcomes.
\begin{theorem}\label{thm:kernel}
Assume Assumptions \ref{asp:dr} and \ref{asp:kernel} hold. We then have the following four are true.
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item\label{thm:kernel1} Assumptions \ref{asp:weight}, \ref{asp:dr2}\ref{asp:dr2,w}\ref{asp:dr2,w2}, and \ref{asp:se2}\ref{asp:se2,w1} hold;
\item\label{thm:kernel2} Under Assumption~\ref{asp:kernel,dr}, Assumption~\ref{asp:dr1}\ref{asp:dr1,w} holds;
\item\label{thm:kernel3} Under Assumptions~\ref{asp:kernel,dr} and \ref{asp:kernel,se}, Assumption~\ref{asp:se2} holds;
\item\label{thm:kernel4} Under Assumptions~\ref{asp:kernel,dr} and \ref{asp:kernel,dml}, Assumption~\ref{asp:dml2} holds with $p_1,p_2$ chosen to be $\infty$ and $1$.
\end{enumerate}
\end{theorem}
Theorem \ref{thm:kernel} directly yields the following corollary, which establishes the double robustness and semiparametric efficiency properties of $\hat\tau_{\rm K}$ and $\tilde{\tau}_{{\rm K},N} $.
\begin{corollary} \phantomsection \label{crl:kernel}
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item(Double robustness of $ \hat{\tau}_{{\rm K}}$) Suppose Assumptions~\ref{asp:dr} and \ref{asp:kernel} hold and either Assumptions~\ref{asp:dr1}\ref{asp:dr1,o}, \ref{asp:kernel,dr} or Assumption~\ref{asp:dr2}\ref{asp:dr2,o} is true. We then have
\begin{align*}
\hat{\tau}_{{\rm K}} - \tau \stackrel{\sf p}{\longrightarrow} 0.
\end{align*}
\item(Semiparametric efficiency of $\hat{\tau}_{{\rm K}}$) Under Assumptions~\ref{asp:dr}, \ref{asp:dr1}\ref{asp:dr1,o}, \ref{asp:dr2}\ref{asp:dr2,o}, \ref{asp:se1}, \ref{asp:kernel}-\ref{asp:kernel,se}, we have
\begin{align*}
\sqrt{n} (\hat{\tau}_{{\rm K}} - \tau) \stackrel{\sf d}{\longrightarrow} N(0,\sigma^2).
\end{align*}
\item(Double robustness of $\tilde{\tau}_{{\rm K},N}$) Suppose Assumptions~\ref{asp:dr} and \ref{asp:kernel} hold and either Assumptions~\ref{asp:dr1}\ref{asp:dr1,o}, \ref{asp:kernel,dr} or Assumption~\ref{asp:dr2}\ref{asp:dr2,o} is true. We then have
\begin{align*}
\tilde{\tau}_{{\rm K},N} - \tau \stackrel{\sf p}{\longrightarrow} 0.
\end{align*}
\item(Semiparametric efficiency of $\tilde{\tau}_{{\rm K},N}$) Under Assumptions~\ref{asp:dr}, \ref{asp:dr1}\ref{asp:dr1,o}, \ref{asp:dr2}\ref{asp:dr2,o}, \ref{asp:dml1}, \ref{asp:kernel}, \ref{asp:kernel,dr}, \ref{asp:kernel,dml}, it holds true that
\begin{align*}
\sqrt{n} (\tilde{\tau}_{{\rm K},N} - \tau) \stackrel{\sf d}{\longrightarrow} N(0,\sigma^2).
\end{align*}
\end{enumerate}
\end{corollary}
\subsection{Weighted NNs}
NN matching \citep{rubin1973matching,abadie2006large,stuart2010matching} is a popular imputation method that imputes the missing potential outcomes by a NN regression. In the nonparametric statistics literature, it is well known that NN regression, which assigns equal weights to all NNs, can be less efficient. This motivates the development of weighted NNs as useful alternatives to NN regression for boosting statistical efficiency \citep{royall1966class,samworth2012optimal}. The theoretical properties of WNNs for conducting nonparametric regression have been studied in, among many others, \cite{stone1977consistent}, \cite{samworth2012optimal}, and \citet[Chapter 5]{biau2015lectures}.
Consider the $M$-NN that restricts attention to the first $M$ NNs. The weighted nearest neighbor (WNN) regression imputes the missing potential outcomes using a set of preassigned weights $[\gamma_{M,m}]_{m=1}^M$ satisfying
\[
\gamma_{M,m}\geq 0~~~{\rm and}~~~\sum_{m=1}^M \gamma_{M,m} = 1.
\]
The corresponding imputed outcome values are then
\begin{align*}
\hat{Y}_i^{\rm WNN}(0) := \begin{cases}
Y_i, & \mbox{ if } D_i=0,\\
\sum_{m=1}^M \gamma_{M,m} Y_{j_m(i)}, & \mbox{ if } D_i=1,
\end{cases}
~~{\rm and}~~ \hat{Y}_i^{\rm WNN}(1) := \begin{cases}
\sum_{m=1}^M \gamma_{M,m} Y_{j_m(i)}, & \mbox{ if } D_i=0,\\
Y_i, & \mbox{ if } D_i=1.
\end{cases}
\end{align*}
Here $j_m(i)$ represents the index of $m$-th nearest neighbor (NN) of $X_i$ in $\{X_j:D_j=1-D_i\}_{j=1}^n$, i.e., the index $j \in \llbracket n\rrbracket$ such that $D_j=1-D_i$ and
\[
\sum_{\ell=1, D_\ell=1-D_i}^n \ind\Big(\lVert X_\ell -X_i \rVert \le \lVert X_j - X_i \rVert\Big) = m.
\]
For any $i,j \in \llbracket n\rrbracket$ with $D_i+D_j=1$, the corresponding weight $w_{i\leftarrow j}$ is then defined to be
\begin{align*}
w_{i\leftarrow j}:= \sum_{m=1}^M \gamma_{M,m} \ind(j_m(i)=j)
\end{align*}
and the WNN-based ATE estimator and its double machine learning version are then denoted by
\[
\hat{\tau}_{\rm WNN}~~~{\rm and}~~~\tilde{\tau}_{{\rm WNN},N}.
\]
Notably speaking, when $\gamma_{M,m}=1/M$ for all $m \in \llbracket M\rrbracket$, $\hat\tau_{\rm WNN}$ reduces to the standard bias-corrected NN matching that was studied in \cite{abadie2006large,abadie2011bias} and \cite{lin2021estimation}.
\begin{assumption} \phantomsection \label{asp:wnn,dr}
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item\label{asp:wnn,dr,1} The density of $X$ is bounded and bounded away from zero. The densities of $X \,|\, D=1$ and $X \,|\, D=0$ are continuous almost everywhere. The diameters and the surface area (Hausdorff measure, \citet[Section 3.3]{evans2018measure}) of $\cS$ are bounded. There exists a constant $a \in (0,1)$ such that for any $\delta \in (0,{\rm diam}(\cS)]$ and $z \in \cS$,
\[
\lambda(B_{z,\delta} \cap S) \ge a \lambda(B_{z,\delta}),
\]
where $B_{z,\delta}$ represents the closed ball in $\bR^d$ with center at $z$ and radius $\delta$.
\item\label{asp:wnn,dr,2} Assume $M\log n/n \to 0$, $\sum_{m=1}^M \gamma_{M,m}^2 \to 0$, and
\begin{align*}
\limsup_{n \to \infty} n \int_0^\infty \Big[\sum_{m=1}^M \gamma_{M,m}^2 {\mathrm P}\Big( U_{(m-1)} \le t \le U_{(m)} \Big) \Big]^{1/2} {\mathrm d} t \le 1,
\end{align*}
where $(U_{(1)},\ldots,U_{(M)})$ are the first $M$ order statistics of $n$ i.i.d random variables from the uniform distribution on $[0,1]$.
\end{enumerate}
\end{assumption}
\begin{assumption} \phantomsection \label{asp:wnn,se}
Assume that there exists a positive integer $k$ such that
\begin{itemize}
\item[(i)] Assumptions~\ref{asp:se2}\ref{asp:se2,o1},\ref{asp:se2,o2} hold;
\item[(ii)] we further have
\[
\sum_{m=1}^M \gamma_{M,m} (m/n)^{k/d} = o(n^{-1/2})~~~{\rm and}~~~\sum_{m=1}^M \gamma_{M,m} (m/n)^{\ell/d} = o(n^{-1/2+ \gamma_\ell})
\]
for all $\ell \in \llbracket k-1\rrbracket$.
\end{itemize}
\end{assumption}
\begin{assumption}[double machine learning] \phantomsection \label{asp:wnn,dml}
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item The densities of $X \,|\, D=1$ and $X \,|\, D=0$ are Lipchitz on $\cS$.
\item Assume $M/\log n \to \infty$ and the weights satisfy
\begin{align*}
M\max_{m \in \llbracket M\rrbracket}\gamma_{M,m} = O(1)
\end{align*}
and there exists a positive sequence $[r_3]=[r_3]_n$ such that
\begin{align*}
n \int_0^\infty \Big[\sum_{m=1}^M \Big(\gamma_{M,m} - \frac{1}{M}\Big)^2 {\mathrm P}\Big( U_{(m-1)} \le t \le U_{(m)} \Big) \Big]^{1/2} {\mathrm d} t = O(r_3).
\end{align*}
Further assume that there exist two positive sequences $[r_1] = [r_1]_n, [r_2] = [r_2]_n$ with $r_1r_2 = o(n^{-1/2})$ such that $\lVert \tilde{\mu}_\omega - \mu_\omega \rVert_\infty = O_{\mathrm P}(r_1)$ for $\omega \in \{0,1\}$, and
\[
(M/n)^{1/d} + M^{-1/2} + \Big(\sum_{m=1}^M \gamma_{M,m}^2\Big)^{1/2} + r_3 = O(r_2).
\]
\end{enumerate}
\end{assumption}
\begin{remark}
Assumption \ref{asp:wnn,dr}\ref{asp:wnn,dr,1} is Assumption 4.1 in \cite{lin2021estimation}. When $\gamma_{M,m}=1/M$ for all $m \in \llbracket M\rrbracket$, Assumption~\ref{asp:wnn,dr}\ref{asp:wnn,dr,2} is satisfied as long as $M \to \infty$, and the inequality in Assumption~\ref{asp:wnn,dr}\ref{asp:wnn,dr,2} can be automatically satisfied by using Chernoff's inequality, which recovers Theorems 5.1 and 5.2 in \cite{lin2021estimation}.
Similar discussions also apply to Assumption \ref{asp:wnn,dml}.
\end{remark}
\begin{theorem}\label{thm:wnn}
Assume Assumption \ref{asp:dr} holds. We then have the following four are true.
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item\label{thm:wnn1} Assumptions \ref{asp:weight}, \ref{asp:dr2}\ref{asp:dr2,w}\ref{asp:dr2,w2}, and \ref{asp:se2}\ref{asp:se2,w1} hold;
\item\label{thm:wnn2} Under Assumption~\ref{asp:wnn,dr}, Assumption~\ref{asp:dr1}\ref{asp:dr1,w} holds;
\item\label{thm:wnn3} Under Assumptions~\ref{asp:wnn,dr} and \ref{asp:wnn,se}, Assumption~\ref{asp:se2} holds;
\item\label{thm:wnn4} Under Assumptions~\ref{asp:wnn,dr} and \ref{asp:wnn,dml}, Assumption~\ref{asp:dml2} holds with $p_1,p_2$ chosen to be $\infty$ and $1$.
\end{enumerate}
\end{theorem}
\begin{corollary} \phantomsection \label{crl,wnn}
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item(Double robustness of $\hat\tau_{\rm WNN}$) Suppose Assumption~\ref{asp:dr} holds, and either Assumptions~\ref{asp:dr1}\ref{asp:dr1,o} and \ref{asp:wnn,dr} or Assumption~\ref{asp:dr2}\ref{asp:dr2,o} is true. We then have
\begin{align*}
\hat{\tau}_{{\rm WNN}} - \tau \stackrel{\sf p}{\longrightarrow} 0.
\end{align*}
\item(Semiparametric efficiency of $\hat\tau_{\rm WNN}$) Under Assumptions~\ref{asp:dr}, \ref{asp:dr1}\ref{asp:dr1,o}, \ref{asp:dr2}\ref{asp:dr2,o}, \ref{asp:se1}, \ref{asp:wnn,dr}, \ref{asp:wnn,se}, we have
\begin{align*}
\sqrt{n} (\hat{\tau}_{{\rm WNN}} - \tau) \stackrel{\sf d}{\longrightarrow} N(0,\sigma^2).
\end{align*}
\item(Double robustness of $\tilde{\tau}_{{\rm WNN},N} $) Suppose Assumption~\ref{asp:dr} holds, and either Assumptions~\ref{asp:dr1}\ref{asp:dr1,o} and \ref{asp:wnn,dr} or Assumption~\ref{asp:dr2}\ref{asp:dr2,o} is true. We then have
\begin{align*}
\tilde{\tau}_{{\rm WNN},N} - \tau \stackrel{\sf p}{\longrightarrow} 0.
\end{align*}
\item(Semiparametric efficiency of $\tilde{\tau}_{{\rm WNN},N}$) Under Assumptions~\ref{asp:dr}, \ref{asp:dr1}\ref{asp:dr1,o}, \ref{asp:dr2}\ref{asp:dr2,o}, \ref{asp:dml1}, \ref{asp:wnn,dr}, \ref{asp:wnn,dml}, it holds true that
\begin{align*}
\sqrt{n} (\tilde{\tau}_{{\rm WNN},N} - \tau) \stackrel{\sf d}{\longrightarrow} N(0,\sigma^2).
\end{align*}
\end{enumerate}
\end{corollary}
\subsection{Local linear matching}
In nonparametric statistics, local linear regression has been a prominent alternative to local constant regression, proving to be more efficient than the latter, especially along the boundary \citep{fan1992design,fan1993local}. This approach has also been heavily used in ATE estimation for imputing the missing potential outcomes, which is often called ``local linear matching''; cf. \cite{heckman1997matching}, \cite{heckman1998characterizing}, \cite{heckman1998matching}, and \cite{frolich2005matching}.
In detail, for any unit $i \in \llbracket n\rrbracket$, local linear matching uses the local linear regression \citep{fan2018local} to minimize
\begin{align}\label{eq:llr}
\sum_{j:D_j=1-D_i} \Big[Y_j - \beta_0 - \beta^\top (X_j - X_i) \Big]^2 K_H(X_j-X_i),
\end{align}
and then $Y_i(1-D_i)$ is imputed by the solution to the above objective function.
Let $\mB_i \in \bR^{n_{1-D_i} \times (1+d)}$ be the design matrix with the row corresponding to unit $j$ with $D_i + D_j=1$ to be $(1,(X_j-X_i)^\top):=b_{ij}^\top$. Let $\mW_i \in \bR^{n_{1-D_i} \times n_{1-D_i}}$ be the diagonal matrix with the diagonal element corresponding to unit $j$ with $D_i + D_j=1$ to be $K_H(X_j-X_i)$. It is well known that the solution to the minimization problem \eqref{eq:llr} is:
\[
\hat{Y}_i^{\rm LL}(1-D_i) = e_1^\top (\mB_i^\top \mW_i \mB_i)^{-1} \mB_i^\top \mW_i \mY_{1-D_i},
\]
where $e_1 \in \bR^{1+d}$ is the vector with the first element to be 1 and all the rest 0 and $\mY_{\omega} \in \bR^{n_\omega}$ for $\omega \in \{0,1\}$ represents the vector containing entries $Y_j$'s with $D_j=\omega$. For any $i,j \in \llbracket n\rrbracket$ with $D_i+D_j=1$, one could calculate the corresponding weight $w_{i\leftarrow j}$ as
\begin{align*}
w_{i\leftarrow j}:= e_1^\top (\mB_i^\top \mW_i \mB_i)^{-1} b_{ij} K_H(X_j-X_i).
\end{align*}
Denote the corresponding local linear estimator and its double machine learning version by
\[
\hat{\tau}_{\rm LL}~~ {\rm and}~~ \tilde{\tau}_{{\rm LL,}N}.
\]
For analyzing $\hat\tau_{\rm LL}$ and $\tilde{\tau}_{{\rm LL,}N}$, we need to regulate the kernel function $K(\cdot)$ a little bit more. The following assumption is standard in multivariate local linear regression literature (cf. Assumption A1 in \cite{ruppert1994multivariate}). It can be satisfied by many kernels, e.g., the spherically symmetric kernels and product kernels based on symmetric univariate kernels \citep[Chapter 4]{simonoff2012smoothing}.
\begin{assumption} \phantomsection \label{asp:lp}
Assume $\int z K(z) {\mathrm d} z = 0$ and $\int z z^\top K(z) {\mathrm d} z = \mu_2(K) I_d$ with $\mu_2(K)>0$ as a positive real-valued constant that captures the second-order property of $K$.
\end{assumption}
\begin{theorem}\label{thm:lp}
Assume Assumptions \ref{asp:dr} and \ref{asp:kernel} hold. We then have the following four are true.
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item\label{thm:lp1} Assumptions~\ref{asp:weight}, \ref{asp:dr2}\ref{asp:dr2,w}, \ref{asp:se2}\ref{asp:se2,w1} hold; if $K(\cdot)$ is bounded with a compact support and is bounded away from zero, then Assumption \ref{asp:dr2}\ref{asp:dr2,w2} holds;
\item\label{thm:lp2} Under Assumptions~\ref{asp:kernel,dr}, \ref{asp:lp}, Assumption~\ref{asp:dr1}\ref{asp:dr1,w} holds;
\item\label{thm:lp3} Under Assumptions~\ref{asp:kernel,dr}, \ref{asp:kernel,se}, \ref{asp:lp}, Assumption~\ref{asp:se2} holds;
\item\label{thm:lp4} Under Assumptions~\ref{asp:kernel,dr}, \ref{asp:kernel,dml}, \ref{asp:lp}, Assumption~\ref{asp:dml2} holds with $p_1,p_2$ chosen to be $\infty$ and $1$.
\end{enumerate}
\end{theorem}
\begin{corollary} \phantomsection \label{crl:lp}
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item(Double robustness of $\hat{\tau}_{\rm LL}$) Suppose Assumptions~\ref{asp:dr} and \ref{asp:kernel} hold and either Assumptions~\ref{asp:dr1}\ref{asp:dr1,o}, \ref{asp:kernel,dr}, \ref{asp:lp} hold or Assumption~\ref{asp:dr2}\ref{asp:dr2,o} is true and $K(\cdot)$ is bounded with a compact support and is bounded away from zero. We then have
\begin{align*}
\hat{\tau}_{\rm LL} - \tau \stackrel{\sf p}{\longrightarrow} 0.
\end{align*}
\item(Semiparametric efficiency of $\hat{\tau}_{\rm LL}$) Under Assumptions~\ref{asp:dr}, \ref{asp:dr1}\ref{asp:dr1,o}, \ref{asp:dr2}\ref{asp:dr2,o}, \ref{asp:se1}, \ref{asp:kernel}, \ref{asp:kernel,dr}, \ref{asp:kernel,se}, \ref{asp:lp}, we have
\begin{align*}
\sqrt{n} (\hat{\tau}_{\rm LL} - \tau) \stackrel{\sf d}{\longrightarrow} N(0,\sigma^2).
\end{align*}
\item(Double robustness of $\tilde{\tau}_{{\rm LL,}N}$) Suppose Assumptions~\ref{asp:dr} and \ref{asp:kernel} hold and either Assumptions~\ref{asp:dr1}\ref{asp:dr1,o}, \ref{asp:kernel,dr}, \ref{asp:lp} hold or Assumption~\ref{asp:dr2}\ref{asp:dr2,o} is true and $K(\cdot)$ is bounded with a compact support and is bounded away from zero. We then have
\begin{align*}
\tilde{\tau}_{{\rm LL,}N} - \tau \stackrel{\sf p}{\longrightarrow} 0.
\end{align*}
\item(Semiparametric efficiency of $\tilde{\tau}_{{\rm LL,}N}$) Under Assumptions~\ref{asp:dr}, \ref{asp:dr1}\ref{asp:dr1,o}, \ref{asp:dr2}\ref{asp:dr2,o}, \ref{asp:dml1}, \ref{asp:kernel}, \ref{asp:kernel,dr}, \ref{asp:kernel,dml}, \ref{asp:lp}, it holds true that
\begin{align*}
\sqrt{n} (\tilde{\tau}_{{\rm LL,}N} - \tau) \stackrel{\sf d}{\longrightarrow} N(0,\sigma^2).
\end{align*}
\end{enumerate}
\end{corollary}
\section{Random forests}\label{sec:RF}
This section studies random forests as an imputation method to estimate the ATE. Since being invented by Leo Breiman \citep{breiman2001random}, random forests have proven to be practically powerful in conducting regression and classification tasks; cf. the survey of \cite{biau2016random}. However, it is not until very recent that some major advances were made towards using random forests for inferring causal effect \citep{athey2019machine}; notable works include \cite{hill2011bayesian}, \cite{athey2016recursive}, \cite{athey2019machine}, \cite{athey2019generalized}, among many others. Our results in this section aim to contribute to this growing literature, while being focused on the original regression-adjusted imputation estimator without doing sample splitting and cross fitting.
\subsection{Set up}
In the sequel, for any set $A$ with finite elements, let $\lvert A \rvert$ stand for its cardinality. For introducing the random forests to impute the missing potential outcomes, some additional notation is needed and we also adopt some common terms used in the random forests and regression trees literature \citep{breiman2017classification}.
Let's first introduce the causal tree. Let $T^1$ be a generic {\it tree} built on the treated group $\{(X_i,Y_i)\}_{i=1,D_i=1}^n$ and $T^0$ be another generic tree built on the control group $\{(X_i,Y_i)\}_{i=1,D_i=0}^n$. The two trees $T^1$ and $T^0$ accordingly partition the covariates space $\cS\subset \bR^d$ into a set of leaves $L^1$ and $L^0$, respectively. For any test point $x \in \bR^d$, let $L^1(x)$ and $L^0(x)$ be the {\it leaves} of $T^1$ and $T^0$ containing $x$. One could then impute the missing potential outcomes as follows:
\begin{align*}
\hat{Y}_i^{\rm Tree}(0) := \begin{cases}
Y_i, & \mbox{ if } D_i=0,\\
\displaystyle \Big(\Big\lvert \Big\{j:D_j=0,X_j \in L^0(X_i)\Big\} \Big\rvert\Big)^{-1}\sum_{j:D_j=0,X_j \in L^0(X_i)} Y_j, & \mbox{ if } D_i=1,
\end{cases}
\end{align*}
and
\begin{align*}
\hat{Y}_i^{\rm Tree}(1) := \begin{cases}
\displaystyle \Big(\Big\lvert \Big\{j:D_j=1,X_j \in L^1(X_i)\Big\} \Big\rvert\Big)^{-1}\sum_{j:D_j=1,X_j \in L^1(X_i)} Y_j, & \mbox{ if } D_i=0,\\
Y_i, & \mbox{ if } D_i=1.
\end{cases}
\end{align*}
To aggregate many individual causal trees into a {\it causal forest}, we consider {\it subsampling}. In detail, let $B$ be the number of trees and $s$ be the subsample size, which for presentation simplicity are assumed to be identical for the two groups of samples. In the $b$-th round, for building the tree, we sample without replacement the following two size-$s$ subsets
\[
{\mathcal I}^1_b ~~~{\rm and}~~~ {\mathcal I}^0_b
\]
from $\{i:D_i=1\}$ and $\{i:D_i=0\}$, respectively. Of note, for any $b,b'\in\llbracket B\rrbracket$ and any $\omega\in\{0,1\}$, ${\mathcal I}^\omega_b$ and ${\mathcal I}^\omega_{b'}$ could have a nonempty overlap.
Let $T^1_b$ be the tree built on the data $\{(X_i,Y_i)\}_{i \in {\mathcal I}^1_b}$ and $T^0_b$ be the tree built on the data $\{(X_i,Y_i)\}_{i \in {\mathcal I}^0_b}$. All trees are assumed to be constructed using the same base learner. For any test point $x$, let $L^1_b(x)$ and $L^0_b(x)$ be the leaves of $T^1_b$ and $T^0_b$ that contain $x$. The according random forest then imputes the missing potential outcomes as follows:
\begin{align*}
\hat{Y}_i^{\rm RF}(0) := \begin{cases}
Y_i, & \mbox{ if } D_i=0,\\
\displaystyle \frac1B \sum_{b=1}^B \Big[\Big(\Big\lvert \Big\{j \in {\mathcal I}^0_b:X_j \in L^0_b(X_i)\Big\} \Big\rvert\Big)^{-1}\sum_{j \in {\mathcal I}^0_b:X_j \in L^0_b(X_i)} Y_j\Big], & \mbox{ if } D_i=1,
\end{cases}
\end{align*}
and
\begin{align*}
\hat{Y}_i^{\rm RF}(1) := \begin{cases}
\displaystyle \frac1B \sum_{b=1}^B \Big[\Big(\Big\lvert \Big\{j \in {\mathcal I}^1_b:X_j \in L^1_b(X_i)\Big\} \Big\rvert\Big)^{-1}\sum_{j \in {\mathcal I}^1_b:X_j \in L^1_b(X_i)} Y_j\Big], & \mbox{ if } D_i=0,\\
Y_i, & \mbox{ if } D_i=1.
\end{cases}
\end{align*}
It is well known that random forests, formulable as a special type of weighted NN regressions, constitute linear smoothers \citep{lin2006random,biau2010layered}. In particular, for any $i,j \in \llbracket n\rrbracket$ with $D_i+D_j=1$, one could verify that the weight $w_{i\leftarrow j}$ corresponding to the above random forests imputation method is
\begin{align*}
w_{i\leftarrow j}:= \frac1B \sum_{b=1}^B \frac{\ind\Big(j \in {\mathcal I}^{1-D_i}_b:X_j \in L^{1-D_i}_b(X_i)\Big)}{\Big\lvert \Big\{k \in {\mathcal I}^{1-D_i}_b:X_k \in L^{1-D_i}_b(X_i)\Big\} \Big\rvert}.
\end{align*}
We then denote the corresponding regression-adjusted random forest-based imputation ATE estimator by $\hat{\tau}_{\rm RF}$.
\subsection{Inference theory}
In order to verify the conditions in Section \ref{sec:general}, the following assumptions are needed and were intentionally designed to be general.
\begin{assumption} \phantomsection \label{asp:rf,dr}
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item\label{asp:rf,dr,d} The density of $X$ is bounded and bounded away from zero, the densities of $X \,|\, D=1$ and $X \,|\, D=0$ are continuous almost everywhere, and the support $\cS$ is compact.
\item\label{asp:rf,dr,t} We assume $s=s_n=O(n^{1/2})$ and $n/B = O(1)$. In addition, assume that for the tree $T$ built on $s$ i.i.d. sampled points from $(X,Y)\,|\, D=1$ or $(X,Y)\,|\, D=0$ with leaves $\{L_t\}_{t \ge 1}$, it holds true that
\begin{align}\label{assump:RF1}
\lim_{n \to \infty} {\mathrm E} \Big[\Big(\min_{t\ge1}\Big\lvert L_t \Big\rvert\Big)^{-1}\Big] = 0~~~{\rm and}~~~ \lim_{n \to \infty} \int_S {\mathrm E}\Big[{\rm diam}\Big(L_t(x) \cap \cS\Big)\Big] {\mathrm d} x = 0,
\end{align}
where $\lvert L_t \rvert$ represents the number of samples in the leaf $L_t$ for $t\ge1$ and $L_t(x)$ stands for the leaf that contains $x$.
\end{enumerate}
\end{assumption}
\begin{assumption} \phantomsection \label{asp:rf,honest}
The tree is honest, that is, the tree does not use the responses $Y_i$'s to choose the place to split.
\end{assumption}
\begin{assumption} \phantomsection \label{asp:rf,se}
There exists a positive integer $k$ such that
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item Assumptions~\ref{asp:se2}\ref{asp:se2,o1}, \ref{asp:se2,o2} hold;
\item\label{asp:rf,se,2} for a tree $T$ built on $s$ independent observations from $(X,Y)\,|\, D=1$ or $(X,Y)\,|\, D=0$ with leaves $\{L_t\}_{t \ge 1}$, it holds true that
\begin{align*}
&\int_S {\mathrm E}\Big[{\rm diam}^k\Big(L_t(x) \cap \cS\Big)\Big] {\mathrm d} x = o(n^{-1/2})~~\\
{\rm and}~~&\int_S {\mathrm E}\Big[{\rm diam}^\ell\Big(L_t(x) \cap \cS\Big)\Big] {\mathrm d} x = o(n^{-1/2+ \gamma_\ell}) \text{ for all }\ell \in \llbracket k-1\rrbracket.
\end{align*}
\end{enumerate}
\end{assumption}
\begin{remark}\label{remark:RF1}
Assumption~\ref{asp:rf,dr} requires $s_n=O(n^{1/2})$. In the literature, \cite{mentch2016quantifying} required a similar condition, $s_n = o(n^{1/2})$, for establishing asymptotic normality of random forests.
\cite{wager2018estimation} allowed $s_n \asymp n^\beta$ for some $\beta$ that can be close to 1 (cf. Equation (14) therein); we cannot recover their setting due to the extra difficulty in estimating the ATE compared to estimating the conditional ATE. Assumption~\ref{asp:rf,dr} also requires $n/B=O(1)$, which echoes \cite{wager2014confidence}, where the authors recommended a similar $B \asymp n$ condition. Conditions similar to the two leaf size conditions in \eqref{assump:RF1} have been discussed in multiple places. There the first requirement in \eqref{assump:RF1} regulates the smallest size of the terminal leaves, which echoes the discussions in \citet[Section 3]{lin2006random}. The second requirement in \eqref{assump:RF1} is very related to \citet[Lemma 1]{wager2018estimation}; we defer more discussions on it as well as those on Assumption \ref{asp:rf,se}\ref{asp:rf,se,2} to Lemma \ref{lemma:rf,dist} and Proposition \ref{prop:rf,dist} ahead.
\end{remark}
\begin{remark}
The ``honesty'' condition, Assumption~\ref{asp:rf,honest}, corresponds to Definition 2 in \cite{wager2018estimation}. This condition is usually achieved by implementing sample splitting as was suggested and also analyzed in \cite{wager2018estimation}. It is also satisfied by a variety of alternatives to Breiman's original random forests, including the centered forest \citep{biau2008consistency,scornet2016asymptotics} and the purely uniform random forests \citep{genuer2012variance}. Theoretical analysis of the trees constructed using the responses in the same training data is believed to be much more involved, but was managed in several impressive works including \cite{scornet2015consistency}, \cite{chi2020asymptotic}, and \cite{kulowski2022}. Unfortunately, our analysis hinges on a control of the leaf sizes that is seemingly hard to pursue without Assumption \ref{asp:rf,honest}.
\end{remark}
Under the above assumptions, we are then ready to present our main theory on $\hat\tau_{\rm RF}$. Note that, in the following, Theorem \ref{thm:rf}\ref{thm:rf-dr} also gives rise to a consistent random forests-based density ratio estimator, which can be of independent interest.
\begin{theorem}\label{thm:rf}
Assume Assumption~\ref{asp:dr} holds. We then have the following four are true.
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item\label{thm:rf1} Assumptions~\ref{asp:weight}, \ref{asp:dr2}\ref{asp:dr2,w2}, \ref{asp:se2}\ref{asp:se2,w1} hold.
\item\label{thm:rf-dr} Under Assumptions~\ref{asp:rf,dr} and \ref{asp:rf,honest}, Assumption~\ref{asp:dr1}\ref{asp:dr1,w} holds.
\item\label{thm:rf2} Under Assumption~\ref{asp:rf,honest}, Assumption \ref{asp:dr2}\ref{asp:dr2,w} holds.
\item\label{thm:rf3} Under Assumptions~\ref{asp:rf,dr}-\ref{asp:rf,se}, Assumption~\ref{asp:se2} holds.
\end{enumerate}
\end{theorem}
\begin{corollary} \phantomsection \label{crl,rf}
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item(Double robustness of $\hat\tau_{\rm RF}$) Suppose Assumption~\ref{asp:dr} holds and either Assumptions \ref{asp:dr1}\ref{asp:dr1,o}, \ref{asp:rf,dr}, \ref{asp:rf,honest} or Assumptions~\ref{asp:dr2}\ref{asp:dr2,o}, \ref{asp:rf,honest} hold. We then have
\begin{align*}
\hat{\tau}_{{\rm RF}} - \tau \stackrel{\sf p}{\longrightarrow} 0.
\end{align*}
\item(Semiparametric efficiency of $\hat\tau_{\rm RF}$) Under Assumptions~\ref{asp:dr}, \ref{asp:dr1}\ref{asp:dr1,o}, \ref{asp:dr2}\ref{asp:dr2,o}, \ref{asp:se1}, \ref{asp:rf,dr}-\ref{asp:rf,se}, we have
\begin{align*}
\sqrt{n} (\hat{\tau}_{{\rm RF}} - \tau) \stackrel{\sf d}{\longrightarrow} N(0,\sigma^2).
\end{align*}
\end{enumerate}
\end{corollary}
\subsection{Balanced and regular random forests}
The goal here is to decipher the second part of \eqref{assump:RF1} and Assumption \ref{asp:rf,se}\ref{asp:rf,se,2}; cf. the discussions in Remark \ref{remark:RF1}. To this end, we leverage the technical proofs of \citet[Lemma 1]{wager2018estimation} and \citet[Lemma 2]{meinshausen2006quantile}, and provide the convergence rates of the diameters of leaves for some particular trees.
To this end, we introduce the following regularity conditions on the tree growing patter.
\begin{assumption}\label{def:tree}
We consider the following type of trees:
\begin{enumerate}[itemsep=-.5ex,label=(\roman*)]
\item The tree is $\phi$-balanced, i.e., for each terminal leaf, the proportion of splits along the $j$-th axis for each $j\in\llbracket d\rrbracket$ is lower bounded by $\phi/d$ for some $\phi \in (0,1)$ and the splitting directions (i.e., picking which feature to split) are independent of the data;
\item The tree is $(\alpha,\theta)$-regular for some $\alpha \in (0,0.5]$ and some positive integer $\theta$, i.e., at each step of growing the tree, the split leaves at least $\alpha$ of the samples on each side of the split, and the terminal leaves are all of size in $[\theta,\lfloor \theta/\alpha \rfloor]$, where $\lfloor
\cdot \rfloor$ is the floor function.
\end{enumerate}
\end{assumption}
Notably speaking, Assumption 3 in \cite{meinshausen2006quantile} and Definitions 3 and 4 in \cite{wager2018estimation} considered regular and random-split conditions that are similar to Assumption \ref{def:tree}. In practice, Assumption \ref{def:tree} can always be satisfied by controlling how tree grows in the implementation.
For those trees that satisfy Assumption~\ref{def:tree}, we have the following lemma, which controls arbitrary finite moment of the diameter of the terminal leaves.
\begin{lemma}\label{lemma:rf,dist}
Let $\epsilon \in (0,1)$, $p \in \llbracket d\rrbracket$, $\cS = [0,1]^d$, and $T$ be a tree constructed based on $s$ i.i.d. observations from the uniform distribution on $\cS$. As long as $T$ is $(\alpha,\theta)$-regular and $\phi$-balanced, we have for any $x \in \cS$ and any positive integer $k$,
\begin{align*}
{\mathrm E}\Big[{\rm diam}_p^k\Big(L_t(x) \cap \cS\Big)\Big] \le \Big(\frac{s}{\alpha^{-1}\theta}\Big)^{k \frac{\log(1-(1-\epsilon)\alpha)}{\log(\alpha^{-1})} \frac{\phi}{d} }+ \frac{\log(s/\theta)}{\log(\alpha^{-1})} \exp\Big[-\theta\alpha\Big(\log\Big(\frac{1}{1-\epsilon}\Big)-\epsilon\Big)\Big],
\end{align*}
where ${\rm diam}_p(\cdot)$ stands for the diameter along the $p$-th axis.
\end{lemma}
Lemma~\ref{lemma:rf,dist} then yields sufficient conditions guaranteeing the validity of the second part of \eqref{assump:RF1} and Assumption \ref{asp:rf,se}\ref{asp:rf,se,2}.
\begin{proposition}[Sufficient conditions on the leaf sizes]\label{prop:rf,dist}
Assume $\cS$ to be a compact subset of $\bR^d$ and $T$ to be a tree constructed based on $s$ i.i.d. observations following a distribution with density bounded and bounded away from zero on $\cS$. Assume further that $T$ is both regular and balanced. We then have, if $s/\theta \to \infty$ and $\log \log (s/\theta)/\theta \to 0$,
\begin{align}\label{eq:RF-db}
\lim_{n \to \infty} \int_S {\mathrm E}\Big[{\rm diam}\Big(L_t(x) \cap \cS\Big)\Big] {\mathrm d} x = 0.
\end{align}
If it further holds that $(s/\theta)/n^\epsilon \to \infty$ for some $\epsilon>0$ and $\theta/\log n \to \infty$, we then have, for any sufficiently large $k$,
\begin{align}\label{eq:RF-se}
\int_S {\mathrm E}\Big[{\rm diam}^k\Big(L_t(x) \cap \cS\Big)\Big] {\mathrm d} x = o(n^{-1/2}).
\end{align}
\end{proposition}
Of note, in Proposition~\ref{prop:rf,dist} the requirements about $s$ and $\theta$ are much weaker for double robustness (corresponding to \eqref{eq:RF-db}) than for semiparametric efficiency (corresponding to \eqref{eq:RF-se}).
\begin{remark}
Lemma~\ref{lemma:rf,dist} is key to our analysis and is a stronger version of Lemma 1 in \cite{wager2018estimation}. In detail, Lemma 1 in \cite{wager2018estimation} or the proof of Theorem 3 therein can imply that the $k$-th moment of the diameter will always be dominated by
\[
(s/\theta)^{-0.5[\log((1-\alpha)^{-1})/\log(\alpha^{-1})](\phi/d)},
\]
which, however, can not be faster than $n^{-1/2}$ for any positive integer $k$. In contrast, Lemma~\ref{lemma:rf,dist} establishes that we can reach the order $o(n^{-1/2})$ by taking $k$ large enough. This is viable by replacing the random-split condition in \cite{wager2018estimation} with Assumption \ref{def:tree}.
\end{remark}
\begin{remark}
It is worth noting that Lemma~\ref{lemma:rf,dist} and Proposition~\ref{prop:rf,dist} do not require the tree to be honest. This is in line with Lemma 2 in \cite{meinshausen2006quantile} for quantile regression tree using the responses and Lemma 1 in \cite{wager2018estimation} without assuming honesty. It indicates that the results in Lemma~\ref{lemma:rf,dist} and Proposition~\ref{prop:rf,dist} can be applied to more general random forests, e.g., the tree based on CART criteria \citep{breiman2017classification} with consistency analyzed in \cite{scornet2015consistency}. However, the ``regular'' and ``random-split'' conditions enforced in Assumption \ref{def:tree} seem inevitable to our analysis. Later, we require honesty for the double robustness and semiparametric efficiency of $\hat\tau_{\rm RF}$.
\end{remark}
\section*{Acknowledgement}
We thank helpful discussions with Peng Ding, Kevin Guo, and Elizabeth Stuart on the matching procedure, and Yingying Fan on the random forest.
{
\bibliographystyle{apalike}
\bibliography{AMS}
}
\newpage{}