EconBase
← Back to paper

On regression-adjusted imputation estimators of the average treatment effect

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

72,594 characters · 16 sections · 126 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

On regression-adjusted imputation estimators of the average treatment effect

{5pt} {5pt} {5pt} {5pt} \hypersetup{colorlinks,breaklinks,urlcolor=blue,linkcolor=blue}

abstractImputing 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 imbens2004nonparametric,stuart2010matching. Generalizing an earlier result of the authors 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.

{\bf Keywords}: double robustness, kernel matching, nearest neighbor matching, random forests, double machine learning.

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 rosenbaum1983central,imbens2015causal. While nowadays a large literature focuses on propensity score-based methods, alternatives that are based on regression heckman1997matching,heckman1998matching,heckman1998characterizing,hahn1998role,athey2016recursive,wager2018estimation and matching 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 rubin2004multiple,tsiatis2006semiparametric,little2019statistical. Indeed, as Guido Imbens and others (cf. imbens2004nonparametric and 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 robins1995semiparametric in the missing value, imbens2004nonparametric and abadie2006large in the causal inference, and cassel1976some and 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 robins1994estimation,scharfstein1999adjusting as well as its much more recent cousin, the double machine learning estimators chernozhukov2018double--- partly in order to encourage more efficient and robust estimators.

This paper is interested in exploring the {\it double robustness} 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 rubin1973use,abadie2011bias, unlike its counterpart that integrates imputation with weighting --- e.g., propensity score robins1994estimation,hirano2003efficient or covariate balancing 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., 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 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, heckman1998matching underpinned theoretical studies of (bias-uncorrected) imputation methods and showed that imputation based on covariate kernel matching yields a semiparametrically efficient estimator. Nevertheless, heckman1998matching's result only focuses on estimating the average treatment effect on the treated (ATT). Later, 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} 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 lin2021estimation, the new theory shows:

itemize• a linear smoother can implicitly give rise to a density ratio estimator; • imputation methods with regression adjustments in the form of rubin1973use and abadie2011bias constitute AIPW estimators; • these imputation methods are consistent as long as either the density model or the outcome model is correctly specified, and thus {\it doubly robust}; • they further constitute asymptotically normal estimators of the ATE with the asymptotic variance attaining the semiparametric efficiency lower bound hahn1998role if both the density and outcome models are correctly specified, and are thus {\it semiparametrically efficient}; • the double machine learning 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.

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 lin2021estimation, who established the double robustness and semiparametrical efficiency theory for 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 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) and the study of more imputation methods elaborated on in Sections (ref) and (ref). Due to the richness of newly obtained results, we feel compelled to disseminate them to peers by writing a second manuscript.

{\bf Paper organization.} Section (ref) introduces necessary notation, the preliminary setup, and those regression-adjusted imputation ATE estimators that will be analyzed in subsequent sections. Section (ref) lays out our general theory, with examples provided in Sections (ref) and (ref). Specifically, Section (ref) concerns imputation using kernel matching, weighted NN, and local linear matching while Section (ref) is focused on imputing the missing potential outcomes using random forests.

Preliminary

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 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 =

casesY_i(0), & if D_i=0,\\ Y_i(1), & if D_i=1.

\] 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,

align*[align* omitted — 61 chars of source]

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:

align*[align* omitted — 189 chars of source]

and

align*[align* omitted — 185 chars of source]

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} buja1989linear. Note that all imputation methods considered in Sections (ref) and (ref), including the kernel regression and local linear regression estimators heckman1997matching,heckman1998characterizing,heckman1998matching, the (weighted) NN regression abadie2006large,abadie2011bias,lin2021estimation, and the (honest) random forests 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 rubin1973use and 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, 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

align*[align* omitted — 220 chars of source]

and

align*[align* omitted — 215 chars of source]

The according regression-adjusted imputation-based ATE estimator is

align*[align* omitted — 94 chars of source]

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.

The general theory

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]. \]

A key lemma

Results in Section (ref) are all built on the following key lemma, which gives an AIPW formulation of the ATE estimator $\hat\tau_w$.

lemmaThe regression-adjusted imputation estimator $\hat\tau_w$ can be rewritten as \begin{align} \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}

The sum of the first three terms in (ref) has the same form as an AIPW estimator that was studied in scharfstein1999adjusting and 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 (ref) 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. fan2018local and wasserman2006all.

Note that the relation between regression-adjusted imputation and AIPW estimators was for the first time disclosed in lin2021estimation, stated as Lemma 5.1 therein and with a focus on NN regression-based imputation. Lemma (ref), on the other hand, delivers the general form that applies to an arbitrary linear smoother.

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.

assumption\phantomsection \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • 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$. • $[(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)$. • ${\mathrm E} [U^2_\omega \,|\, X=x] $ is uniformly bounded for almost all $x \in \cS$ and $\omega \in \{0,1\}$. • ${\mathrm E} [\mu^2_\omega(X)]$ is bounded for $\omega \in \{0,1\}$. \end{enumerate}

Assumption (ref)(ref) 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 rosenbaum1983central. The rest conditions in Assumption (ref) 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.

assumption\phantomsection \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • 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)}$. • 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}

Assumption (ref) is to our knowledge new and is added for aiding the general theory to be presented later. There Assumption (ref)(ref) ensures that the regression smoothing matrix is invariant to the feeding order of sample points, and Assumption (ref)(ref) ensures that the bias term in Lemma (ref) is asymptotically ignorable, which will be automatically satisfied if the smoother preserves the constant curve wasserman2006all.

The next set of assumptions quantifies estimation accuracy of the “density models”.

assumption\phantomsection \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • 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). \] • 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}

Assumption (ref) allows for outcome model misspecification. Here Assumption (ref)(ref) is a regression misspecification assumption that is Assumption 5.3 in lin2021estimation. Assumption (ref)(ref) 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) and (ref) we will verify its validity for a variety of regression imputation methods.

In parallel to Assumption (ref), the following conditions quantify estimation accuracy of the “outcome models”.

assumption\phantomsection \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • For $\omega \in \{0,1\}$, the estimator $\hat{\mu}_\omega(x)$ satisfies \[ \lVert \hat{\mu}_\omega - \mu_\omega \rVert_\infty = o_{\mathrm P}(1). \] • 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$. • 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}

Assumption (ref) allows for density model misspecification. Here Assumption (ref)(ref) is Assumption 5.4 in lin2021estimation; chen2015optimal and chen2018optimal verified such conditions for various nonparametric regressors. Assumption (ref)(ref) ensures that the responses are not used in the construction of weights, and is satisfied by all examples to be introduced in Sections (ref) and (ref). This assumption is also related to the sample splitting procedures used in the context of double machine learning chernozhukov2018double and honest random forests wager2018estimation, shown to help avoid overfitting. Given Assumptions (ref) and (ref), Assumption (ref)(ref) holds automatically as long as all the weights $w_{i\leftarrow j}$'s are nonnegative, or when Assumption (ref)(ref) holds. We would also like to highlight that Assumption (ref)(ref) 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$.

theorem[Double robustness of $\hat\tau_w$] Suppose Assumptions (ref) and (ref) hold, and either Assumption (ref) or Assumption (ref) is true. We then have \begin{align*} \hat\tau_w - \tau \stackrel{\sf p}{\longrightarrow} 0. \end{align*}

Theorem (ref) 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)(ref). 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) resides in simple parametric models.

In detail, the fact that ordinary least square (OLS) is intrinsically a weighted estimator is very well known; cf. angrist2009mostly and imbens2015matching. In two very interesting papers, robins2007comment and 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 chattopadhyay2021implied and other interesting research along this line includes guo2021generalized and cohen2020no. In the high level, they all bear a similar flavor to Theorem (ref) 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.

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.

assumption\phantomsection \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • ${\mathrm E} [U^2_\omega \,|\, X=x]$ is uniformly bounded away from zero for almost all $x \in \cS$ and $\omega \in \{0,1\}$. • 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}
assumption\phantomsection There exists a positive integer $k$ such that \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • $\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$; • 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$; • 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} ) \rm for all \ell \in \llbracket k-1\rrbracket; \end{align*} • 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}

Assumptions (ref) and (ref)(ref)-(ref) are Assumptions 5.6 and 5.7 in lin2021estimation; check abadie2011bias and chen2018optimal for results on verifying these requirements. Assumption (ref)(ref) 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)(ref), as a counterpart of Assumption (ref)(ref), requires the bias term in (ref) to be root-$n$ ignorable.

We then introduce the semiparametric efficiency lower bound for estimating the ATE hahn1998role,

align[align omitted — 157 chars of source]

The following theorem then shows that the asymptotic variance of $\hat\tau_w$ can attain $\sigma^2$.

theorem[Semiparametric efficiency of $\hat\tau_w$] Suppose Assumptions (ref)-(ref) 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 (ref).

Double machine learning

Assumptions (ref)(ref)-(ref) are arguably strong regularity conditions for the outcome model $\mu_\omega(\cdot)$. Partly in order to alleviate such requirements, chernozhukov2018double introduced the idea of double machine learning via sample splitting and cross fitting. Similar ideas have also been studied in nonparametric statistics; cf. bickel1982adaptive, efromovich1996nonparametric, and zheng2010asymptotic. In the following, let's introduce $\tilde{\tau}_{w,N}$ as a counterpart of $\hat{\tau}_w$ based on 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:

align*[align* omitted — 290 chars of source]

and

align*[align* omitted — 88 chars of source]

For establishing the efficiency theory of $\tilde{\tau}_{w,N} $, the following two sets of assumptions are needed.

assumption\phantomsection \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • ${\mathrm E} [U^2_\omega]$ is bounded away from zero for $\omega \in \{0,1\}$. • There exists some constant $\kappa>0$ such that ${\mathrm E} [\lvert Y \rvert^{2+\kappa}]$ is bounded. \end{enumerate}
assumptionThere 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*)] • 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); \] • 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}

We are now ready to introduce the general theory on the double machine learning-based regression-adjusted imputation estimators.

theorem\phantomsection \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • (Double robustness of $\tilde{\tau}_{w,N}$) Under the same conditions as those in Theorem (ref), we have \begin{align*} \tilde{\tau}_{w,N} - \tau \stackrel{\sf p}{\longrightarrow} 0. \end{align*} • (Semiparametric efficiency of $\tilde{\tau}_{w,N}$) Under Assumptions (ref)-(ref) and (ref)-(ref), 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 (ref). \end{enumerate}

Examples

This section aims to provide examples so to put the general theory introduced in Section (ref) 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.

Kernel matching

We first consider the kernel matching that has been advocated in various settings heckman1997matching, heckman1998characterizing, heckman1998matching, frolich2004finite, frolich2005matching, huber2013performance. It leverages the local constant regression (Nadaraya–Watson estimator) to impute the missing values 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

align*[align* omitted — 91 chars of source]

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) can then be shown to hold under the following sufficient conditions.

assumption\phantomsection Assume that (i) $H$ is symmetric and positive definite, and (ii) $K$ constitutes a multivariate symmetric density function.
assumption\phantomsection \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • 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. • $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}
assumption\phantomsection There exists a positive integer $k$ such that \begin{itemize} • Assumptions (ref)(ref),(ref) hold; • 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}
assumption[double machine learning] \phantomsection \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • The densities of $X \,|\, D=1$ and $X \,|\, D=0$ are Lipchitz on $\cS$. The diameter and the surface area (Hausdorff measure, evans2018measure) of $\cS$ are bounded. • 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}

Assumption (ref) is standard for establishing consistency of the Nadaraya-Watson estimator. Assumption (ref) ensures that the discrepancy level in Assumption (ref) is small. The regularity condition on the support and the smoothness condition on the density function are standard in nonparametric statistics MR2724359.

remarkA 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) then reduces to \[ h_n \to 0\quad {\rm and}\quad nh_n^d \to \infty, \] and Assumption (ref) 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 lin2021estimation). The convergence rate in Assumption (ref) reduces to $(nh^d)^{-1/2}+h$, and is the minimax rate of the density estimation over Lipchitz class $n^{-1/(2+d)}$ MR2724359 by taking $h_n \asymp n^{-1/(2+d)}$.

The following theorem then verifies the general conditions presented in Section (ref) when kernel matching is used for imputing the missing potential outcomes.

theoremAssume Assumptions (ref) and (ref) hold. We then have the following four are true. \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • Assumptions (ref), (ref)(ref)(ref), and (ref)(ref) hold; • Under Assumption (ref), Assumption (ref)(ref) holds; • Under Assumptions (ref) and (ref), Assumption (ref) holds; • Under Assumptions (ref) and (ref), Assumption (ref) holds with $p_1,p_2$ chosen to be $\infty$ and $1$. \end{enumerate}

Theorem (ref) 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} $.

corollary\phantomsection \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • (Double robustness of $ \hat{\tau}_{{\rm K}}$) Suppose Assumptions (ref) and (ref) hold and either Assumptions (ref)(ref), (ref) or Assumption (ref)(ref) is true. We then have \begin{align*} \hat{\tau}_{{\rm K}} - \tau \stackrel{\sf p}{\longrightarrow} 0. \end{align*} • (Semiparametric efficiency of $\hat{\tau}_{{\rm K}}$) Under Assumptions (ref), (ref)(ref), (ref)(ref), (ref), (ref)-(ref), we have \begin{align*} \sqrt{n} (\hat{\tau}_{{\rm K}} - \tau) \stackrel{\sf d}{\longrightarrow} N(0,\sigma^2). \end{align*} • (Double robustness of $\tilde{\tau}_{{\rm K},N}$) Suppose Assumptions (ref) and (ref) hold and either Assumptions (ref)(ref), (ref) or Assumption (ref)(ref) is true. We then have \begin{align*} \tilde{\tau}_{{\rm K},N} - \tau \stackrel{\sf p}{\longrightarrow} 0. \end{align*} • (Semiparametric efficiency of $\tilde{\tau}_{{\rm K},N}$) Under Assumptions (ref), (ref)(ref), (ref)(ref), (ref), (ref), (ref), (ref), 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}

Weighted NNs

NN matching 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 royall1966class,samworth2012optimal. The theoretical properties of WNNs for conducting nonparametric regression have been studied in, among many others, stone1977consistent, samworth2012optimal, and 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

align*[align* omitted — 352 chars of source]

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

align*[align* omitted — 78 chars of source]

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 abadie2006large,abadie2011bias and lin2021estimation.

assumption\phantomsection \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • 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, 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$. • 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}
assumption\phantomsection Assume that there exists a positive integer $k$ such that \begin{itemize} • Assumptions (ref)(ref),(ref) hold; • 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}
assumption[double machine learning] \phantomsection \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • The densities of $X \,|\, D=1$ and $X \,|\, D=0$ are Lipchitz on $\cS$. • 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}
remarkAssumption (ref)(ref) is Assumption 4.1 in lin2021estimation. When $\gamma_{M,m}=1/M$ for all $m \in \llbracket M\rrbracket$, Assumption (ref)(ref) is satisfied as long as $M \to \infty$, and the inequality in Assumption (ref)(ref) can be automatically satisfied by using Chernoff's inequality, which recovers Theorems 5.1 and 5.2 in lin2021estimation. Similar discussions also apply to Assumption (ref).
theoremAssume Assumption (ref) holds. We then have the following four are true. \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • Assumptions (ref), (ref)(ref)(ref), and (ref)(ref) hold; • Under Assumption (ref), Assumption (ref)(ref) holds; • Under Assumptions (ref) and (ref), Assumption (ref) holds; • Under Assumptions (ref) and (ref), Assumption (ref) holds with $p_1,p_2$ chosen to be $\infty$ and $1$. \end{enumerate}
corollary\phantomsection \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • (Double robustness of $\hat\tau_{\rm WNN}$) Suppose Assumption (ref) holds, and either Assumptions (ref)(ref) and (ref) or Assumption (ref)(ref) is true. We then have \begin{align*} \hat{\tau}_{{\rm WNN}} - \tau \stackrel{\sf p}{\longrightarrow} 0. \end{align*} • (Semiparametric efficiency of $\hat\tau_{\rm WNN}$) Under Assumptions (ref), (ref)(ref), (ref)(ref), (ref), (ref), (ref), we have \begin{align*} \sqrt{n} (\hat{\tau}_{{\rm WNN}} - \tau) \stackrel{\sf d}{\longrightarrow} N(0,\sigma^2). \end{align*} • (Double robustness of $\tilde{\tau}_{{\rm WNN},N} $) Suppose Assumption (ref) holds, and either Assumptions (ref)(ref) and (ref) or Assumption (ref)(ref) is true. We then have \begin{align*} \tilde{\tau}_{{\rm WNN},N} - \tau \stackrel{\sf p}{\longrightarrow} 0. \end{align*} • (Semiparametric efficiency of $\tilde{\tau}_{{\rm WNN},N}$) Under Assumptions (ref), (ref)(ref), (ref)(ref), (ref), (ref), (ref), 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}

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 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. heckman1997matching, heckman1998characterizing, heckman1998matching, and frolich2005matching.

In detail, for any unit $i \in \llbracket n\rrbracket$, local linear matching uses the local linear regression fan2018local to minimize

align[align omitted — 115 chars of source]

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 (ref) 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

align*[align* omitted — 97 chars of source]

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 ruppert1994multivariate). It can be satisfied by many kernels, e.g., the spherically symmetric kernels and product kernels based on symmetric univariate kernels simonoff2012smoothing.

assumption\phantomsection 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$.
theoremAssume Assumptions (ref) and (ref) hold. We then have the following four are true. \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • Assumptions (ref), (ref)(ref), (ref)(ref) hold; if $K(\cdot)$ is bounded with a compact support and is bounded away from zero, then Assumption (ref)(ref) holds; • Under Assumptions (ref), (ref), Assumption (ref)(ref) holds; • Under Assumptions (ref), (ref), (ref), Assumption (ref) holds; • Under Assumptions (ref), (ref), (ref), Assumption (ref) holds with $p_1,p_2$ chosen to be $\infty$ and $1$. \end{enumerate}
corollary\phantomsection \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • (Double robustness of $\hat{\tau}_{\rm LL}$) Suppose Assumptions (ref) and (ref) hold and either Assumptions (ref)(ref), (ref), (ref) hold or Assumption (ref)(ref) 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*} • (Semiparametric efficiency of $\hat{\tau}_{\rm LL}$) Under Assumptions (ref), (ref)(ref), (ref)(ref), (ref), (ref), (ref), (ref), (ref), we have \begin{align*} \sqrt{n} (\hat{\tau}_{\rm LL} - \tau) \stackrel{\sf d}{\longrightarrow} N(0,\sigma^2). \end{align*} • (Double robustness of $\tilde{\tau}_{{\rm LL,}N}$) Suppose Assumptions (ref) and (ref) hold and either Assumptions (ref)(ref), (ref), (ref) hold or Assumption (ref)(ref) 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*} • (Semiparametric efficiency of $\tilde{\tau}_{{\rm LL,}N}$) Under Assumptions (ref), (ref)(ref), (ref)(ref), (ref), (ref), (ref), (ref), (ref), 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}

Random forests

This section studies random forests as an imputation method to estimate the ATE. Since being invented by Leo Breiman breiman2001random, random forests have proven to be practically powerful in conducting regression and classification tasks; cf. the survey of biau2016random. However, it is not until very recent that some major advances were made towards using random forests for inferring causal effect athey2019machine; notable works include hill2011bayesian, athey2016recursive, athey2019machine, 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.

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 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:

align*[align* omitted — 262 chars of source]

and

align*[align* omitted — 262 chars of source]

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:

align*[align* omitted — 325 chars of source]

and

align*[align* omitted — 325 chars of source]

It is well known that random forests, formulable as a special type of weighted NN regressions, constitute linear smoothers 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

align*[align* omitted — 223 chars of source]

We then denote the corresponding regression-adjusted random forest-based imputation ATE estimator by $\hat{\tau}_{\rm RF}$.

Inference theory

In order to verify the conditions in Section (ref), the following assumptions are needed and were intentionally designed to be general.

assumption\phantomsection \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • 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. • 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} \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}
assumption\phantomsection The tree is honest, that is, the tree does not use the responses $Y_i$'s to choose the place to split.
assumption\phantomsection There exists a positive integer $k$ such that \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • Assumptions (ref)(ref), (ref) hold; • 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}) for all \ell \in \llbracket k-1\rrbracket. \end{align*} \end{enumerate}
remarkAssumption (ref) requires $s_n=O(n^{1/2})$. In the literature, mentch2016quantifying required a similar condition, $s_n = o(n^{1/2})$, for establishing asymptotic normality of random forests. 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) also requires $n/B=O(1)$, which echoes wager2014confidence, where the authors recommended a similar $B \asymp n$ condition. Conditions similar to the two leaf size conditions in (ref) have been discussed in multiple places. There the first requirement in (ref) regulates the smallest size of the terminal leaves, which echoes the discussions in lin2006random. The second requirement in (ref) is very related to wager2018estimation; we defer more discussions on it as well as those on Assumption (ref)(ref) to Lemma (ref) and Proposition (ref) ahead.
remarkThe “honesty” condition, Assumption (ref), corresponds to Definition 2 in wager2018estimation. This condition is usually achieved by implementing sample splitting as was suggested and also analyzed in wager2018estimation. It is also satisfied by a variety of alternatives to Breiman's original random forests, including the centered forest biau2008consistency,scornet2016asymptotics and the purely uniform random forests 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 scornet2015consistency, chi2020asymptotic, and kulowski2022. Unfortunately, our analysis hinges on a control of the leaf sizes that is seemingly hard to pursue without Assumption (ref).

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)(ref) also gives rise to a consistent random forests-based density ratio estimator, which can be of independent interest.

theoremAssume Assumption (ref) holds. We then have the following four are true. \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • Assumptions (ref), (ref)(ref), (ref)(ref) hold. • Under Assumptions (ref) and (ref), Assumption (ref)(ref) holds. • Under Assumption (ref), Assumption (ref)(ref) holds. • Under Assumptions (ref)-(ref), Assumption (ref) holds. \end{enumerate}
corollary\phantomsection \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • (Double robustness of $\hat\tau_{\rm RF}$) Suppose Assumption (ref) holds and either Assumptions (ref)(ref), (ref), (ref) or Assumptions (ref)(ref), (ref) hold. We then have \begin{align*} \hat{\tau}_{{\rm RF}} - \tau \stackrel{\sf p}{\longrightarrow} 0. \end{align*} • (Semiparametric efficiency of $\hat\tau_{\rm RF}$) Under Assumptions (ref), (ref)(ref), (ref)(ref), (ref), (ref)-(ref), we have \begin{align*} \sqrt{n} (\hat{\tau}_{{\rm RF}} - \tau) \stackrel{\sf d}{\longrightarrow} N(0,\sigma^2). \end{align*} \end{enumerate}

Balanced and regular random forests

The goal here is to decipher the second part of (ref) and Assumption (ref)(ref); cf. the discussions in Remark (ref). To this end, we leverage the technical proofs of wager2018estimation and 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.

assumptionWe consider the following type of trees: \begin{enumerate}[itemsep=-.5ex,label=(\roman*)] • 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; • 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}

Notably speaking, Assumption 3 in meinshausen2006quantile and Definitions 3 and 4 in wager2018estimation considered regular and random-split conditions that are similar to Assumption (ref). In practice, Assumption (ref) can always be satisfied by controlling how tree grows in the implementation.

For those trees that satisfy Assumption (ref), we have the following lemma, which controls arbitrary finite moment of the diameter of the terminal leaves.

lemmaLet $\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.

Lemma (ref) then yields sufficient conditions guaranteeing the validity of the second part of (ref) and Assumption (ref)(ref).

proposition[Sufficient conditions on the leaf sizes] 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} \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} \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}

Of note, in Proposition (ref) the requirements about $s$ and $\theta$ are much weaker for double robustness (corresponding to (ref)) than for semiparametric efficiency (corresponding to (ref)).

remarkLemma (ref) is key to our analysis and is a stronger version of Lemma 1 in wager2018estimation. In detail, Lemma 1 in 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) 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 wager2018estimation with Assumption (ref).
remarkIt is worth noting that Lemma (ref) and Proposition (ref) do not require the tree to be honest. This is in line with Lemma 2 in meinshausen2006quantile for quantile regression tree using the responses and Lemma 1 in wager2018estimation without assuming honesty. It indicates that the results in Lemma (ref) and Proposition (ref) can be applied to more general random forests, e.g., the tree based on CART criteria breiman2017classification with consistency analyzed in scornet2015consistency. However, the “regular” and “random-split” conditions enforced in Assumption (ref) seem inevitable to our analysis. Later, we require honesty for the double robustness and semiparametric efficiency of $\hat\tau_{\rm RF}$.

Acknowledgement

We thank helpful discussions with Peng Ding, Kevin Guo, and Elizabeth Stuart on the matching procedure, and Yingying Fan on the random forest.