Finite-Sample Optimal Estimation and Inference on Average Treatment Effects Under Unconfoundedness

Finite-Sample Optimal Estimation and Inference on Average Treatment Effects Under Unconfoundedness

Introduction

To estimate the \ac{ATE} of a binary treatment in observational studies, it is often assumed that the treatment is unconfounded given a set of pretreatment covariates. This assumption implies that systematic differences in outcomes between treated and control units with the same values of the covariates are attributable to the treatment. When the covariates are continuously distributed, it is not possible to perfectly match the treated and control units based on their covariate values, and estimation of the \ac{ATE} requires nonparametric regularization methods such as kernel, series or sieve estimators, or matching estimators that allow for imperfect matches.

Many of these estimators are $\sqrt{n}$-consistent, asymptotically unbiased and normally distributed, provided that, in addition to unconfoundedness, one also assumes overlap of the covariate distributions of treated and untreated subpopulations, as well as enough smoothness of either the propensity score or the conditional mean of the outcome given the treatment and covariates hahn_role_1998,heckman_matching_1998,hirano_efficient_2003,chen_semiparametric_2008. The standard approach to constructing \acfp{CI}, which we refer to as $\sqrt{n}$-inference, is to take any such estimator, and add and subtract its standard deviation times the conventional 1.96 critical value (for nominal 95% \acp{CI}).

However, in many applications, this approach may not perform well for two reasons. First, often the overlap is limited, which may lead to failure of asymptotic normality and slower convergence rates KhTa10,KhNe13,busso_new_2014,rothe_robust_2017.\footnote{To deal with limited overlap, one can redefine the estimand to a treatment effect for a subset of the population for which overlap holds, as in heckman_matching_1997, galiani_water_2005, bailey_war_2015 or crump_dealing_2009. However, this estimand is typically less policy relevant.} Second, even under perfect overlap, asymptotic unbiasedness requires a lot of smoothness: one typically needs continuous differentiability of the order $p/2$ at minimum chen_semiparametric_2008, and often of the order $p+1$ or higher hahn_role_1998,heckman_matching_1998,hirano_efficient_2003, where $p$ is the dimension of the covariates. Unless $p$ is very small, such assumptions are hard to evaluate, and may be much stronger than the researcher is willing to impose. Furthermore, even if the asymptotic bias is negligible, the actual finite-sample bias may affect coverage of standard \acp{CI} robins_toward_1997.

In this paper, we instead treat the smoothness and/or shape restrictions on the conditional mean of the outcome given the treatment and covariates as given and determined by the researcher. Given these restrictions, we show how to construct finite-sample valid \acp{CI} based on any estimator $\sum_{i=1}^{n}k_{i}Y_{i}$ that is linear in the outcomes $Y_{i}$, where the weights $\{k_{i}\}_{i=1}^{n}$ depend on the covariates and treatments for the entire sample.\footnote{If the restrictions are asymmetric, we allow for the estimators to be affine rather than linear, taking the form $a+\sum_{i=1}^{n}k_{i}Y_{i}$ where the intercept $a$ may also depend on the covariates and treatment. See (ref).} To do so, we assume that the regression errors are normal with known variance, and we view the treatment and covariates as fixed. The latter allows us to explicitly calculate the worst-case finite-sample bias of the estimator under the maintained restrictions on the conditional mean. Our \acp{CI} are constructed by simply adding and subtracting the estimator's standard error times a critical value that is larger than the usual $1.96$ critical value, and takes into account the potential bias of the estimator.\footnote{The worst-case bias calculations require the researcher to fully specify the restrictions on the conditional mean, including any smoothness constants. The results in (ref) imply that a priori specification of the smoothness constants is unavoidable, and we therefore recommend reporting \acp{CI} for a range of smoothness constants as a form of sensitivity analysis.} We then show how to choose the weights $\{k_{i}\}_{i=1}^{n}$ optimally, in order to minimize the \ac{RMSE} of the estimator, or the length of the resulting \ac{CI}: one needs to solve a finite-sample bias-variance tradeoff problem, which can be cast as a convex programming problem. Furthermore, we show that, once the weights are optimized, such estimators and \acp{CI} are highly efficient among all procedures.

To make further progress on characterizing the optimal weights, we focus on the case where the conditional mean is assumed to satisfy a Lipschitz constraint. We show that, for a given sample size, the optimal estimator (both for \ac{RMSE} and \ac{CI} length) reduces to a matching estimator with a single match when the Lipschitz constant is large enough. While the optimal estimator does not admit a closed form in general, we develop a computationally fast algorithm that traces out the weights as a function of the Lipschitz constant, analogous to the least angle regression algorithm for the LASSO solution path ehjt04.

When the assumption of normal errors and known variance is dropped, feasible versions of our \acp{CI} are valid asymptotically, uniformly in the underlying distribution li89honest. Importantly, this result obtains whether $\sqrt{n}$-inference is possible or whether it is impossible, be it due to low regularity of the regression function\footnote{We show that for $\sqrt{n}$-inference to be possible, one needs to assume a bound on the derivative of the conditional mean of order at least $p/2$. If one only bounds derivatives of lower order, the bias will asymptotically dominate the variance---in contrast to, say, estimation of a conditional mean at a point, it is not possible to “undersmooth”.}, limited overlap, or complete lack of overlap; we do not even require that the \ac{ATE} be point identified. As we further discuss in (ref), our efficiency theory is analogous to the classic theory in the linear regression model, where finite-sample optimality results rely on the strong assumption of normal homoskedastic errors, but asymptotic validity of \acp{CI} based on heteroskedasticity robust standard errors obtains under weak assumptions.

The key condition underlying the asymptotic validity of our \acp{CI} is that the estimator doesn't put too much weight $k_{i}$ on any individual observation: this ensures that the estimator is asymptotically normal when normalized by its standard deviation. However, since the weights solve a bias-variance tradeoff, no single observation can receive too much weight---otherwise, in large samples, regardless of the amount of overlap, a large decrease in variance could be achieved at a small cost to bias. On the other hand, asymptotic normality may fail under limited overlap for other estimators: we show that asymptotic normality for matching estimators generally requires strong overlap.

We illustrate the results in an application to the \ac{NSW} Demonstration. We find that our \acp{CI} are substantially different from those based on conventional $\sqrt{n}$-asymptotic theory, with bias determining a substantial portion of the \ac{CI} width. In line with the theoretical results above, we also find evidence that, in contrast to the estimator we propose, the weights for the matching estimator are too large for a normal approximation to its sampling distribution to be reliable.

Our results rely on the key insight that, if we condition on the treatment and covariates, the \ac{ATE} is a linear functional of a regression function. This puts the problem in the general framework of donoho94 and CaLo04 and allows us to apply the sharp efficiency bounds in ArKo18optimal. The form of the optimal estimator and \acp{CI} follows by applying the general framework. The rest of our finite-sample results, as well as all asymptotic results, are novel and require substantial further analysis. In particular, solving for the optimal weights $k_{i}$ in general requires solving an optimization problem over the space of functions in $p$ variables, which makes simple strategies, such as gridding, infeasible unless the dimension of covariates $p$ is very small. We show that under Lipschitz smoothness, the problem can be recast so that the computational complexity depends only on $n$ and not on $p$, and our solution algorithm uses insights from rosset_piecewise_2007 to further speed up the computation. In independent and contemporaneous work, kallus2017 computes optimal linear weights using a different characterization of the optimization problem.

In contrast, without conditioning on the treatment and covariates, the problem is more difficult: while upper and lower bounds for the rate of convergence have been developed robins_semiparametric_2009, efficiency bounds that are sharp in finite samples remain elusive. Whether one should condition on the treatment and covariates when evaluating estimators and \acp{CI} is itself an interesting question. A previous version of this paper ArKo18ate considers this question in the context of our empirical application, and AbImZh14,aaiw14 give a discussion in related settings. Since our \acp{CI} are valid unconditionally, they can be used in either setting.\footnote{While we focus on a treatment effect that conditions on realized covariates in the sample, our approach can also be used to construct \acp{CI} for the population \ac{ATE}\@; see (ref).}

The remainder of this paper is organized as follows. (ref) presents the model and gives the main finite-sample results. (ref) considers practical implementation issues. (ref) presents asymptotic results. (ref) discusses an application to the \ac{NSW}. Additional results, details, and proofs are collected in the appendices and the supplemental materials.

Setup and finite-sample results

This section sets up the model, and shows how to construct finite-sample optimal estimators and well as finite-sample valid and optimal \acp{CI} under general smoothness restrictions on the conditional mean of the outcome. We then specialize the results to the case with Lipschitz smoothness. Proofs and additional details are given in (ref).

Setup

We have a random sample of $n$ units. Each unit $i=1, \dotsc, n$ is characterized by a pair of potential outcomes $Y_{i}(0)$ and $Y_{i}(1)$ under no treatment and treatment, respectively, a covariate vector $X_{i}\in\mathbb{R}^{p}$, and a treatment indicator $D_{i}\in\{0,1\}$. Unless stated otherwise, we condition on the realized values $\{x_{i}, d_{i}\}_{i=1}^{n}$ of the covariates and treatment $\{X_{i}, D_{i}\}_{i=1}^{n}$, so that probability statements are with respect to the conditional distribution of $\{Y_{i}(0), Y_{i}(1)\}_{i=1}^{n}$. The realized outcome is given by $Y_{i}=Y_{i}(1)d_{i}+Y_{i}(0)(1-d_{i})$. Letting $f(x_{i}, d_{i})$ denote the conditional mean of $Y_{i}$, we obtain a fixed design regression model

We are interested in the \acf{CATE}\footnote{We note that the terminology varies in the literature. Some papers call this object the sample average treatment effect (SATE); other papers use the terms \ac{CATE} and SATE for different objects entirely.}, which, under the assumption of unconfoundedness, $(Y_{i}(1), Y_{i}(0))\protect\mathpalette{\protect\independenT}{\perp} D_{i}\mid X_{i}$, is given by

To obtain finite-sample results, we further assume that $u_{i}$ is normal,

and that the variance function $\sigma^2(x_i, d_i)$ is known. We relax both assumptions in (ref).

We assume that $f$ lies in a known function class $\mathcal{F}$, which formalizes the “regularity” or “smoothness” that we are willing to impose. We require that $\mathcal{F}$ be convex and centrosymmetric, i.e.\ that $f\in\mathcal{F}$ implies $-f\in\mathcal{F}$. While convexity is essential for most of our results, centrosymmetry can be relaxed---see (ref). Our setup covers classical nonparametric function classes, which place bounds on (possibly higher order) derivatives of $f$. As a leading example, we place Lipschitz constraints on $f(\cdot,0)$ and $f(\cdot,1)$:

where $\norm{\cdot}_{\mathcal{X}}$ is a norm on $x$, and $C$ denotes the Lipschitz constant, which for simplicity we take to be the same for both $f(\cdot, 1)$ and $f(\cdot, 0)$. As we discuss in (ref), the function class $\mathcal{F}$ needs to be specified ex ante: for the Lipschitz case, for instance, one cannot use data-driven procedures to estimate $C$, or pick the norm $\norm{\cdot}_{\mathcal{X}}$. We discuss how to make these choices in practice in the context of our empirical application in (ref).

Our goal is to construct estimators and confidence sets for the \ac{CATE} parameter $Lf$. Letting $P_{f}$ denote the probability computed under $f$, a set $\mathcal{C}$ is a $100\cdot (1-\alpha)\%$ confidence set for $Lf$ if

Linear estimators

We start by showing how to construct \acp{CI} based on estimators that are linear in the outcomes,

For now, we treat the weights $k$ as given---in (ref), we will show how to choose them optimally for a general class of criteria that include \ac{RMSE} and \ac{CI} length. In (ref) and (ref), we show that, if we choose the weights optimally, the resulting estimator and \acp{CI} are optimal or near optimal among all procedures, including nonlinear ones. The class of linear estimators covers many estimators that are popular in practice, such as series or kernel estimators, or various matching estimators.\footnote{Nonlinear estimators include those based on regression trees and artificial neural networks, as well as those using nonlinear thresholding to perform variable selection; see donoho_minimax_1998 for a discussion of cases where, in contrast to the present setting, nonlinear estimators outperform linear estimators.} For example, the matching estimator with $M$ matches that matches (with replacement) on covariates takes the form $L\hat{f}_{M}$, where $\hat{f}_{M}(x_{i}, d_{i})=Y_{i}$, and $\hat{f}_{M}(x_{i},1-d_{i})=\sum_{j=1}^{n}W_{M, ij} Y_{j}$. Here $W_{M, ij}=1/M$ if $j$ is among the $M$ observations with treatment status $d_{j}=1-d_{i}$ that are the closest to $i$ (using the norm $\norm{\cdot}_{\mathcal{X}}$), and zero otherwise. For this estimator, the weights take the form

where $K_{M}(i)=M\sum_{j=1}^{n}W_{M, ji}$ is the number of times observation $i$ is matched.

The estimator $\hat{L}_{k}$ is normally distributed with variance $\operatorname{sd}(\hat L_{k})^2=\sum_{i=1}^{n} k(x_i, d_i)^2\sigma^2(x_i, d_i)$ and maximum bias

By centrosymmetry of $\mathcal{F}$, the minimum bias is given by $-\operatorname{\overline{bias}}_{\mathcal{F}}(\hat{L}_{k})$.

We form one-sided $100\cdot (1-\alpha)\%$ \acp{CI} based on $\hat L_{k}$, as

where $z_{1-\alpha}$ is the $1-\alpha$ quantile of a standard normal distribution. Subtracting the maximum bias, in addition to subtracting $z_{1-\alpha}\operatorname{sd}(\hat L_{k})$, is necessary to prevent undercoverage.

One could form a two-sided \ac{CI} centered at $\hat{L}_{k}$ by adding and subtracting $\operatorname{\overline{bias}}_{\mathcal{F}}(\hat{L}_{k})+z_{1-\alpha/2}\operatorname{sd}(\hat{L}_{k})$. However, this is conservative since the bias cannot be equal to $\operatorname{\overline{bias}}_{\mathcal{F}}(\hat{L}_{k})$ and to $-\operatorname{\overline{bias}}_{\mathcal{F}}(\hat{L}_{k})$ at once. Instead, observe that under any $f\in\mathcal{F}$, the $z$-statistic $(\hat L_{k}-Lf)/\operatorname{sd}(\hat L_{k})$ is distributed $N(t,1)$, where $t=E_{f}(\hat{L}_{k}-Lf)/\operatorname{sd}(\hat L_{k})$, and $t$ is bounded in absolute value by $b=\operatorname{\overline{bias}}_{\mathcal{F}}(\hat{L}_{k})/\operatorname{sd}(\hat L_{k})$, the ratio of the worst-case bias to standard deviation. Thus, denoting the $1-\alpha$ quantile of a $\abs{N(b,1)}$ distribution by $\operatorname{cv}_{\alpha}(b)$, a two-sided \ac{CI} can be formed as

If $\hat{L}_{k}$ is unbiased, the critical value reduces to the conventional critical value: $\operatorname{cv}_{\alpha}(0)=z_{1-\alpha/2}$. If $b>0$, it will be larger: for $b\geq 1.5$ and $\alpha\leq 0.2$, $\operatorname{cv}_{\alpha}(b)\approx b+z_{1-\alpha}$ up to three decimal places,\footnote{The critical value $\operatorname{cv}_{1-\alpha}(b)$ can be computed as the square root of the $1-\alpha$ quantile of a non-central $\chi^{2}$ distribution with $1$ degree of freedom and non-centrality parameter $b^{2}$.} Following donoho94, we refer to the \ac{CI} in (ref) as a \ac{FLCI}, since its length does not depend on the realized outcomes, only on the known variance function $\sigma^{2}(\cdot, \cdot)$ and on $\{x_{i}, d_{i}\}_{i=1}^{n}$.

Optimal linear estimators and \texorpdfstring{\acp{CI}}{CIs}

We now show how to choose the weights $k$ optimally. To that end, we need to define the criteria that we wish to optimize. To evaluate estimators, we consider their maximum \acf{RMSE},

One-sided \acp{CI} can be compared using quantiles of excess length (see (ref)). Finally, to evaluate \acp{FLCI}, we simply consider their length, $2\operatorname{cv}_{\alpha}(\operatorname{\overline{bias}}_{\mathcal{F}}(\hat{L}_{k})/\operatorname{sd}(\hat{L}_{k}))\cdot \operatorname{sd}(\hat{L}_{k})$. Since the length is fixed---it doesn't depend on the data $\{Y_{i}\}_{i=1}^{n}$---choosing the weights $k$ to minimize the length does not affect the coverage properties of the resulting \ac{CI}.

Both performance criteria---\ac{FLCI} and \ac{RMSE}---depend on $k$ only through $\operatorname{\overline{bias}}_{\mathcal{F}}(\hat L_{k})$, and $\operatorname{sd}(\hat L_{k})$, and they are increasing in both quantities (this is also true for performance of one-sided \acp{CI}; see (ref)). Therefore, to find the optimal weights, it suffices to first find weights that minimize the worst-case bias $\operatorname{\overline{bias}}_{\mathcal{F}}(\hat{L}_{k})$ subject to a bound on variance. We can then vary the bound to find the optimal bias-variance tradeoff for a given performance criterion (\ac{FLCI} or \ac{RMSE}). It follows from donoho94 and low95 that this bias-variance frontier can be traced out by solving a convex optimization indexed by a parameter $\delta$ that plays a role analogous to that of a bandwidth; it can be thought of as indexing the relative weight on the variance. Varying $\delta$ then traces out the optimal bias-variance frontier.

For a simple statement of the Donoho-Low result, assume that the parameter space $\mathcal{F}$, in addition to being convex and centrosymmetric, does not restrict the value of \ac{CATE} in the sense that the function $\iota_{\kappa}(x, d)=\kappa d$ lies in $\mathcal{F}$ for all $\kappa\in\mathbb{R}$ (see (ref) for a general statement)\footnote{We also assume the regularity condition that if $\lambda f+\iota_{\kappa}\in\mathcal{F}$ for all $0\leq \lambda< 1$, then $f+\iota_{\kappa}\in \mathcal{F}$. Since $L\iota_{\kappa}=\kappa$, $\{\iota_{\kappa}\}_{\kappa\in\mathbb{R}}$ is the smoothest set of functions that span the potential values of the \ac{CATE} parameter, so this assumption will typically hold if the possible values of $Lf$ are unrestricted.}. Given $\delta>0$, let $f^{*}_{\delta}$ solve

With a slight abuse of notation, define

Then the maximum bias of $\hat{L}_{\delta}$ occurs at $-f^{*}_{\delta}$, and the minimum bias occurs at $f^{*}_{\delta}$, so that

Also, $\hat{L}_{\delta}$ minimizes the worst-case bias among all linear estimators with variance bounded by $\operatorname{sd}(\hat{L}_{\delta})^{2}=\delta^{2}/(2\sum_{j=1}^{n}d_{j}f^*_\delta(x_{j}, d_{j})/\sigma^2(x_{j}, d_{j}))^{2}$. Thus, the estimators $\{\hat{L}_{\delta}\}_{\delta>0}$ trace out the optimal bias-variance frontier.

The weights leading to the shortest possible \ac{FLCI} are given by $k^{*}_{\delta_{\textnormal{FLCI}}}$, where $\delta_{\textnormal{FLCI}}$ minimizes $\operatorname{cv}_{\alpha}(\operatorname{\overline{bias}}_{\mathcal{F}}(\hat{L}_{\delta})/\operatorname{sd}(\hat{L}_{\delta}))\cdot\operatorname{sd}(\hat{L}_{\delta})$ over $\delta$. Similarly, the optimal weights for estimation are given by $k^{*}_{\delta_{\textnormal{RMSE}}}$, where $\delta_{\textnormal{RMSE}}$ minimizes $\operatorname{\overline{bias}}_{\mathcal{F}}(\hat{L}_{\delta})^{2}+\operatorname{sd}(\hat{L}_{\delta})^{2}$.

Estimators and \texorpdfstring{\acp{CI}}{CIs} under Lipschitz smoothness

Computing a \ac{FLCI} based on a linear estimator $\hat{L}_{k}$ with a given set of weights $k$ requires computing the worst-case bias (ref). Computing the \ac{RMSE}-optimal estimator, and the optimal \ac{FLCI} requires solving the optimization problem (ref), and then varying $\delta$ to find the optimal bias-variance tradeoff. Both optimization problems require optimizing over the set $\mathcal{F}$, which, in nonparametric settings, is infinite-dimensional. We now focus on the Lipschitz class $\mathcal{F}=\mathcal{F}_{\textnormal{Lip}}(C)$, and show that in this case, the solution to (ref) can be found by solving a finite-dimensional linear program. The optimization problem (ref) can be cast as a finite-dimensional convex program. Furthermore, if the program is put into a Lagrangian form, then the solution is a piecewise linear function of the Lagrange multiplier, and one can trace the entire solution path $\{\hat{L}_{\delta}\}_{\delta>0}$ using an algorithm similar to the LASSO/LAR algorithm of ehjt04.

We leverage the fact that in both optimization problems, we can identify $f$ with the vector $(f(x_1,0), \dotsc, f(x_n,0), f(x_1,1), \dotsc, \allowbreak f(x_n,1))' \in \mathbb{R}^{2n}$, and replace the functional constraint $f\in\mathcal{F}=\mathcal{F}_{\textnormal{Lip}}(C)$ with $2n(n-1)$ linear inequality constraints

This fact follows from the observation that both in (ref) and in (ref), the objective and constraints depend on $f$ only through its value at these $2n$ points, and from the result that if the Lipschitz constraints hold at these points, it is always possible to find a function $f\in\mathcal{F}_{\textnormal{Lip}}(C)$ that interpolates these points beliakov_interpolation_2006.

The assumption that $\hat{L}_{k}$ satisfies (ref) is necessary to prevent the bias from becoming arbitrarily large at multiples of $f(x, d)=d$ and $f(x, d)=1-d$. (ref) implies that the formulas for one-sided \acp{CI} and two-sided \acp{FLCI} given in (ref) hold with $\operatorname{\overline{bias}}_{\mathcal{F}_{\textnormal{Lip}}(C)}(\hat{L}_{k})$ given in (ref). The last part of the theorem says that it suffices to impose at most $2n_{0}n_{1}$ of the constraints in (ref), where $n_{d}$ is the number of observations with $d_{i}=d$. The condition on the weights $k$ holds, for example, for the matching estimator given in (ref). Since for the matching estimator $k(x_{i}, d_{i})=(2d_{i}-1)/n$ if observation $i$ is not used as a match, the theorem says that one only needs to impose (ref) for pairs of observations with opposite treatment status for which one of the observations is used as a match.

For \ac{RMSE}-optimal estimators and optimal \acp{FLCI}, we have the following result:

(ref) reduces the infinite-dimensional program (ref) to a quadratic optimization problem in $\mathbb{R}^{2n}$ with one quadratic and $2n(n-1)$ linear constraints, and a linear objective function. The computational difficulty can be shown to be polynomial in $n$, and it does not depend on the covariate dimension $p$. If the variance is homoskedastic for each treatment group, then the number of linear constraints can be reduced to $2n_{0}n_{1}$, and the entire solution path can be computed efficiently using the piecewise linear algorithm given in (ref). As a result, implementing the estimator is quite fast: the main specification in the empirical application in (ref) takes less than a minute on a laptop computer.

As we discuss in more detail in (ref), it follows from the algorithm that, similarly to the matching estimator (see (ref)), the optimal estimator takes the form $\hat{L}_{\delta}=L\hat{f}_{\delta}$, where $\hat{f}_{\delta}(x_{i}, d_{i})=Y_{i}$, and $\hat{f}_{\delta}(x_{i},1-d_{i})=\sum_{j=1}^{n}W_{\delta, ij}Y_{j}$, and the weights $W_{\delta, ij}$ correspond to the Lagrange multipliers associated with the constraints (ref) for $d=d_{i}$, scaled to sum to one, $\sum_{j=1}^{n}W_{\delta, ij}=1$. The weights are zero unless $d_{j}=1-d_{i}$ and $j$ is close to $i$ according to a matrix of “effective distances.” The “effective distance” between $i$ and $j$ increases with the total weight $\sum_{i=1}^{n}W_{\delta, ij}$ that we already put on $j$. Thus, we may interpret observations $j$ with non-zero weight $W_{\delta, ij}$ as being “matched” to $i$. The number of matches varies across observations $i$, increases with $\delta$, and depends on the number of units with opposite treatment status that are close to $i$ according to the matrix of effective distances. Observations for which there are more good matches receive relatively more matches, since this decreases the variance of the estimator at little cost in terms of bias. On the other hand, the weight $k^{*}_{\delta}(x_{j}, d_{j})=\frac{1}{n}(1-2d_{j})(1+\sum_{i=1}^{n}W_{\delta, ij})$ on $j$ increases with the number of times it has been used as a match, which increases the variance of the estimator. Using the “effective distance” matrix trades off this increase in the variance against an increase in the bias that results from using a lower-quality match instead.

If the constant $C$ is large enough, the increase in the bias from using more than a single match for each $i$ is greater than any reduction in the variance of the estimator, and the optimal estimator takes the form of a matching estimator with a single match:

The single closest match condition will hold with probability one if conditional on each treatment value, at least one of the covariates $x_{i}$ is drawn from a continuous distribution. In contemporaneous work, kallus2017 gives a similar result using a different method of proof. In the other direction, as $C\to 0$, $\hat L_{\delta_{\textnormal{RMSE}}}$ and $\hat L_{\delta_{\textnormal{FLCI}}}$ both converge to the difference-in-means estimator that compares the average outcomes for the treated and untreated units.

(ref) does not imply that one should choose a large $C$ simply to justify matching with a single match as an optimal estimator: the chosen value of $C$ should instead represent a priori bounds on the smoothness of $f$ formulated by the researcher. Nonetheless, if it is difficult to formulate such bounds, a conservative choice of $C$ may be appropriate. If the choice is conservative enough, (ref) will be relevant.

For the optimality result in (ref), it is important that the metric on $x$ used to define the matching estimator is the same as that used to define the Lipschitz constraint. This formalizes the argument in zhao_using_2004 that conditions on the regression function should be considered when defining the metric used for matching. In the supplemental materials, we illustrate the efficiency loss from matching with the “wrong” metric in the context of our empirical application presented in (ref).

A disadvantage of imposing the Lipschitz condition directly on $f$ is that it rules out the simple linear model for $f$, unless we impose a priori bounds on the magnitude of the regression coefficients. In (ref), we give analogs of (ref) under an alternative specification for $\mathcal{F}$ that imposes the Lipschitz condition on $f$ after partialling out the best linear predictor (and thus allows for unrestricted linear response). We show that in this case, if the constant $C$ is large enough, the optimal estimator takes the form of a regression-adjusted matching estimator with a single match.

Adaptation bounds and optimality among nonlinear procedures

The results in (ref) and (ref) show how to construct \ac{RMSE}-optimal linear estimators, and the shortest \ac{FLCI} based on a linear estimator. Are these results useful, or do they overly restrict the class of procedures?

For estimation, (ref) in (ref) shows that the estimator $\hat{L}_{\delta_{\textnormal{RMSE}}}$ achieving the lowest \ac{RMSE} in the class of linear estimators is also highly efficient among all estimators: one cannot reduce the \ac{RMSE} by more than 10.6% by considering non-linear estimators in general, and, in particular applications, its efficiency can be shown to be even higher.

For \acp{FLCI}, an even stronger result obtains, addressing two concerns. First, their length is determined by the least-favorable function in $\mathcal{F}$ (that maximizes the potential bias), which may result in \acp{CI} that are “too long” when $f$ turns out to be smooth. Consequently, one may prefer a variable-length \ac{CI} that optimizes its expected length over a class of smoother functions $\mathcal{G}\subset\mathcal{F}$ (while maintaining coverage over all of $\mathcal{F}$), especially if this leads to substantial reduction in expected length when $f\in\mathcal{G}$. When such a \ac{CI} also simultaneously achieves near-optimal length over all of $\mathcal{F}$, it is referred to as “adaptive.” A related second concern is that implementing our \acp{CI} in practice requires the user to explicitly specify the parameter space $\mathcal{F}$, which in the case $\mathcal{F}=\mathcal{F}_{\textnormal{Lip}}(C)$ includes the Lipschitz constant $C$ and the norm $\norm{\cdot}_{\mathcal{X}}$. This rules out data-driven procedures that try to implicitly or explicitly pick $C$ or the norm using the data.

To address these concerns, we show in (ref) that attempts to form adaptive \acp{CI} cannot substantively improve upon the \acp{FLCI} we propose. In particular, (ref) gives a sharp bound on the length of a confidence set that optimizes its expected length at a smooth function of the form $g(x, d)=\kappa_{0}+\kappa_{1}d$, while maintaining coverage over the original parameter space $\mathcal{F}$. The sharp bound follows from general results in ArKo18optimal, and it gives a benchmark for the scope for improvement over the \ac{FLCI} centered at $\hat{L}_{\delta_{\textnormal{FLCI}}}$ (the \namecref{theorem:adaptation-theorem} also gives an analogous result for one-sided \acp{CI}). The \namecref{theorem:adaptation-theorem} also gives a universal lower bound for this sharp bound, which evaluates to 71.7% when $1-\alpha=0.95$. The sharp bound depends on the realized values of $\{x_{i}, d_{i}\}_{i=1}^{n}$ and the variance function $\sigma^{2}(\cdot, \cdot)$, and can be explicitly computed in a given application. We find that it typically is much higher than this lower bound. For example, in our empirical application in (ref), the \ac{FLCI} efficiency is over 97% at such smooth functions $g$ in our baseline specification. This implies that there is very little scope for improvement over the \ac{FLCI}\@.

Consequently, data-driven or adaptive methods for constructing \acp{CI} must either fail to meaningfully improve over the \ac{FLCI}, or else undercover for some $f\in\mathcal{F}$. It is thus not possible to, say, estimate the order of differentiability of $f$, or to estimate the Lipschitz constant $C$ for the purposes of forming a tighter \ac{CI}---the parameter space $\mathcal{F}$, including any smoothness constants, must be specified ex ante by the researcher. Because of this, by way of sensitivity analysis, since the \ac{CI} may undercover if $C$ is chosen too small in the sense that $\mathcal{F}_{\textnormal{Lip}}(C)$ excludes the true $f$, we recommend reporting estimates and \acp{CI} for a range of choices of $C$ to see how assumptions about the parameter space affect the results (the estimates and \acp{CI} vary with $C$ because $C$ affects both the optimal tuning parameters $\delta_{\textnormal{RMSE}}$ and $\delta_{\textnormal{FLCI}}$, and the critical value, via the worst-case bias). We adopt this approach in the empirical application in (ref). This mirrors the common practice of reporting results for different specifications of the regression function in parametric regression problems.

The key assumption needed for these efficiency bounds is that $\mathcal{F}$ be convex and centrosymmetric. This holds for $\mathcal{F}_{\textnormal{Lip}}(C)$, and, more generally, for parameter spaces that place bounds on derivatives of $f$. If additional restrictions such as monotonicity are used that break either convexity or centrosymmetry, then some degree of adaptation may be possible. While we leave the full exploration of this question for future research, we note that the approach in (ref) can still be used when the centrosymmetry assumption is dropped. As an example, a previous version of this paper ArKo18ate shows to implement our approach when $\mathcal{F}$ imposes Lipschitz and monotonicity constraints.

Practical implementation

We now discuss implementation of feasible versions of the estimators and \acp{CI} given in (ref) when the variance function $\sigma^{2}(x, d)$ is unknown and the errors $u_{i}$ may be non-normal. We discuss the optimality and validity of these feasible procedures. In (ref), we show how to use our approach to form \acp{CI} for the \ac{PATE}.

Baseline implementation

As a baseline, we propose the following implementation of our procedure:\footnote{An R package implementing this procedure, including an implementation of the piecewise linear algorithm, is available at \url{https://github.com/kolesarm/ATEHonest}.}

The conditional variance, can, for instance, be estimated with the nearest-neighbor variance estimator of AbIm06match $\hat{u}_{i}=J/(J+1)\cdot (Y_{i}-\hat{f}(x_i, d_i))^{2}$, where $\hat{f}(x_{i}, d_{i})$ is the average outcome of $J$ observations (excluding $i$) with treatment status $d_{i}$ that are closest to $i$ according to some distance. Note that since the length of the feasible \ac{CI} in (ref) depends on the variance estimates, it is no longer fixed, in contrast to the infeasible \ac{FLCI}.

In general, the point estimate $\tilde{L}_{\tilde{\delta}_{\textnormal{RMSE}}}$ will differ from the estimate $\tilde{L}_{\tilde{\delta}_{\textnormal{FLCI}}}$ used to form the \ac{CI}. Since reporting multiple estimates can be cumbersome, one can simply compute the \ac{CI} (ref) at $\delta=\tilde{\delta}_{\textnormal{RMSE}}$. The \ac{CI} will then be based on the same estimator reported as a point estimate. While this leads to some efficiency loss, in our main specification in the empirical application in (ref), we find that the resulting \ac{CI} is only 1.2% longer than the one that reoptimizes $\delta$ for \ac{CI} length.

\texorpdfstring{\ac{CI}}{CI} for the population average treatment effect

We now show how our approach can be adapted to construct \acp{CI} for the \ac{PATE} based on linear estimators $\hat{L}_{k}$ of the form (ref). To do so, we treat the covariates and treatment as random. Under random sampling, the \ac{PATE} is given by $\tau=E[Y_{i}(1)-Y_{i}(0)]=E[Lf]$, where $Lf$ is the \ac{CATE} as in (ref).

If we view $\hat{L}_{k}$ as an estimator of $\tau$, the quantity $\operatorname{\overline{bias}}_{\mathcal{F}}(\hat{L}_{k})$ now represents its worst-case bias conditional on $\{X_{i}, D_{i}\}_{i=1}^{n}$, and is therefore random under i.i.d.\ sampling. We thus cannot use the arguments underlying the construction in (ref). Instead, we simply add and subtract $\operatorname{\overline{bias}}_{\mathcal{F}}(\hat{L}_{k})$, in addition to adding and subtracting the usual critical value times a standard error based on the marginal, rather than conditional, variance of $\hat{L}_{k}$,

Here $\operatorname{se}_{\ensuremath\tau}(\hat{L}_{k})^{2}=\operatorname{se}(\hat{L}_{k})^{2}+\operatorname{se}(Lf)^{2}$, where $\operatorname{se}(\hat{L}_{k})^{2}$ is the conditional variance of the estimator using the weights $k(\cdot)$ in (ref), and $\operatorname{se}(Lf)^{2}$ is a consistent estimator of the variance of $Lf$, $\frac{1}{n}E[(f(X_{i}, 1)-f(X_{i}, 0)-\tau)^{2}]$. For the latter, if we use the nearest neighbor variance estimator $\hat{u}_{i}^{2}$ in (ref), and the linear estimator takes the form $\hat{L}_{k}=L\hat{f}$, where $\hat{f}(x_{i}, d_{i})=Y_{i}$ and $\hat{f}(x_{i},1-d_{i})=\sum_{j=1}^{n}W_{ij}Y_{j}$, where $\sum_{j=1}^{n}W_{ij}=1$ and $k_{j}=(1-2d_{j})(1+\sum_{i=1}^{n}W_{ij})/n$ (as discussed in (ref) and following (ref), this includes matching estimators, as well as the optimal estimator), one may use the nearest neighbor estimator suggested by AbIm06match, $\frac{1}{n^{2}}\sum_{i=1}^{n}(\hat{f}(x_{i},1)-\hat{f}(x_{i},0)-\hat{L}_{k})^{2}-\frac{1}{n^{2}} \sum_{i=1}^{n}(1+ \sum_{j=1}^{n}W_{ji}^{2})\hat{u}_{i}^{2}$. We note that the optimality results discussed in (ref) do not apply to the \ac{CI} in (ref). In (ref), we provide formal asymptotic coverage results for this \ac{CI} and its one-sided analog.

Asymptotic results

In (ref), we show that $\sqrt{n}$-inference is impossible when the dimension of covariates is high enough relative to the order of smoothness imposed by $\mathcal{F}$, regardless of how the \ac{CI} is formed. (ref) gives conditions for asymptotic validity of the feasible \acp{CI} when the variance function is unknown and the errors $u_{i}$ may be non-normal. (ref) discusses conditions for asymptotic validity and optimality of \acp{CI} based on matching estimators.

Impossibility of \texorpdfstring{$\sqrt{n}$}{root-n}-inference under low smoothness

As discussed in the introduction, the standard approach to inference, which we refer to as $\sqrt{n}$-inference, is based on estimators that are $\sqrt{n}$-consistent and asymptotically normal, with asymptotically negligible bias. We now show that if the dimension of the (continuously distributed) covariates $p$ is high enough relative to the smoothness of $\mathcal{F}$, this approach is infeasible.

To state the result, let $\Sigma(\gamma, C)$ denote the set of $\ell$-times differentiable functions $f$ such that, for all integers $k_1, k_2, \dotsc, k_{p}$ with $\sum_{j=1}^{p}k_j=\ell$, $\abs*{\frac{\partial^{\ell}f(x)}{\partial x_1^{k_1}\cdots \partial x_{p}^{k_{p}}} -\frac{\partial^\ell f(x')}{\partial x_1^{k_1}\cdots \partial x_{p}^{k_{p}}}} \le C\norm{x-x'}^{\gamma-\ell}$, where $\ell$ is the greatest integer strictly less than $\gamma$ and $\norm{x}^{2}=\sum_{j=1}^{p}x_{j}^{2}$. Note that $f\in\mathcal{F}_{\textnormal{Lip}}(C)$ is equivalent to $f(\cdot,1), f(\cdot,0)\in\Sigma(1,C)$.

The theorem shows that the excess length of a \ac{CI} with conditional coverage in the class with $f(\cdot,0), f(\cdot,1)\in\Sigma(\gamma, C)$ must be of order at least $n^{-\gamma/p}$, even at the “smooth” function $f(x, d)=0$. The Lipschitz case corresponds to $\gamma=1$, so that $\sqrt{n}$-inference is possible only when $p\le 2$. While (ref) considers a setting with normal errors, the same bound applies if the normality assumption is dropped (so long as the class of possible distributions for $u_{i}$ includes the normal distribution), since including other distributions only makes the problem more difficult. (ref) requires coverage conditional on the realizations of the covariates and treatment. For unconditional coverage, the results of robins_semiparametric_2009 imply that if $e\in\Sigma(\gamma_{e}, C)$, where $e(x)=P(D_{i}=1\mid X_{i}=x)$ denotes the propensity score, and if $f(\cdot,0), f(\cdot,1)\in\Sigma(\gamma, C)$, then $\sqrt{n}$-inference is impossible unless $\gamma_{e}+\gamma\ge p/2$. Thus, conditioning effectively takes away any role of smoothness of the propensity score.

Asymptotic validity of feasible \texorpdfstring{\acp{CI}}{CIs}

The following theorem gives sufficient conditions for the asymptotic validity of the feasible \acp{CI} given in (ref) based on the estimator $\tilde L_{\delta}$ when $\mathcal{F}=\mathcal{F}_{\textnormal{Lip}}(C)$. To allow us to capture cases in which the finite-sample bias of the estimator is non-negligible even though the asymptotic bias is negligible under standard asymptotics where the parameter space is fixed with the sample size, we allow $C=C_{n}\to\infty$ as $n\to\infty$. For concreteness, we restrict attention to standard errors based on nearest neighbor estimates, and for simplicity, we focus on the case where the preliminary variance estimate $\tilde{\sigma}^{2}(x, d)$ is non-random, leaving the extension to random $\tilde{\sigma}(x, d)$ to future research.

The main regularity condition, (ref), only requires that the covariate distribution for the treated and control units is well-behaved: we do not require overlap between these two distributions. If $C_n=C$ does not change with $n$, a bounded support condition suffices:

Thus, by (ref), feasible \acp{CI} are asymptotically valid in settings in which $\sqrt{n}$-inference is impossible, including in settings in which the covariate dimension $p$ is high enough so that (ref) applies, and settings with imperfect overlap KhTa10 including set identification due to complete lack of overlap. The estimator $\tilde{L}_{\delta}$ remains asymptotically normal, when normalized by its standard deviation: the “irregular” nature of the setting only shows up through non-negligible asymptotic bias, which is captured by the critical value $\operatorname{cv}_{\alpha}$ when constructing the \ac{CI}, and through a slower rate of convergence of the estimator (so the impossibility result of (ref) is not contradicted).

Asymptotic properties of \texorpdfstring{\acp{CI}}{CIs} based on matching estimators

For feasible \acp{CI} based on matching estimators, we obtain the following result:

(ref) is related to results of AbIm06match on asymptotic properties of matching estimators with a fixed number of matches. AbIm06match note that, when $p$ is large enough, the bias term will dominate, so that conventional \acp{CI} based on matching estimators will not be valid. In contrast, the \acp{CI} in (ref) remain valid even when $p$ is large, since they are widened to take into account the potential bias of the estimator. Alternatively, one can attempt to restore asymptotic coverage by subtracting an estimate of the bias based on higher-order smoothness assumptions. While this can lead to asymptotic validity when additional smoothness is available abadie_bias-corrected_2011, it follows from (ref) that such an approach will lead to asymptotic undercoverage under some sequence of regression functions in the Lipschitz class $\mathcal{F}_{\textnormal{Lip}}(C)$.

Relative to (ref), (ref) requires the additional condition (ref). This condition holds almost surely if $(X_{i}, D_{i})$ are drawn i.i.d.\ from a distribution where $\underline G(a)$ and $\overline G(a)$ are lower and upper bounds (up to constants) for $P(\norm{X_{i}-x}_{\mathcal{X}} \le a, \, D_{i}=d)$ for $x$ on the support of $X_{i}$ pollard_convergence_1984. The condition $\lim_{t\to 0}\overline G(\underline G^{-1}(t))^2/[t/\log t^{-1}]=0$ can thus be interpreted as an overlap condition. In particular, if the density of $X_{i}$ is bounded away from zero and infinity on a sufficiently regular support, then the condition holds if the propensity score $P(D_{i}=1\mid X_{i})$ is bounded away from zero and one.

In the supplemental materials, we strengthen the conclusions of (ref) and give an asymptotic analog of (ref) showing that, if there is sufficient overlap and $p> 2$, matching with $M=1$ matches is asymptotically efficient under the \ac{RMSE} and \ac{FLCI} criteria. Thus, while it may at first appear that, as argued in AbIm06match, the matching estimator is inefficient due to its slower than $\sqrt{n}$ rate of convergence, by (ref), the $\sqrt{n}$-rate is not feasible in this setting: the matching estimator in fact achieves the fastest possible rate, and, when $M=1$, the constant is also asymptotically optimal.

On the other hand, if there is not sufficient overlap, asymptotic normality may fail for the matching estimator. As an extreme example, suppose that $p=1$ and that $x_{j}<x_{i}$ for all observations where $d_{j}=0$ and $d_{i}=1$. Then each untreated observation will be matched to the same treated observation, the one with the smallest value of $x_{i}$ among treated observations. Consequently, the Lindeberg weight defined in (ref) will be bounded away from zero for this observation, and the \ac{CLT} will fail. In contrast, by (ref), the estimator $\tilde{L}_{\delta}$ will be asymptotically normal when scaled by its standard deviation even when there is no overlap between the distribution of $x_{i}$ for treated and untreated observations.

rothe_robust_2017 argues that, in settings with limited overlap, estimators of the \ac{CATE} may put a large amount of weight on a small number of observations. As a result, standard approaches to inference that rely on normal asymptotic approximations to the distribution of the $t$-statistic will be inaccurate in finite samples. Our results shed light on when such concerns are relevant. The above example shows that such concerns may indeed persist---even in large samples---if one uses a matching estimator with a fixed number of matches. Similarly to the discussion in (ref), in finite samples, one can assess these concerns directly by computing the maximum Lindeberg weight $\operatorname{Lind}(k_{\text{match}, M})$. Furthermore, it follows from the proof of (ref) that when $p>2$, bias will dominate variance asymptotically even if one attempts to “undersmooth” by using a matching estimator with a single match. In such settings, it is important to widen the \acp{CI} to take the bias into account, in addition to accounting for the potential inaccuracy of the normal asymptotic approximation.

Empirical application

We now illustrate our methods with an application to the \acf{NSW} demonstration. We use the same dataset as dw99 and abadie_bias-corrected_2011.\footnote{We use the data from Rajeev Dehejia's website, \url{http://users.nber.org/ rdehejia/nswdata2.html}.} In particular, the treated sample corresponds to 185 men in the \ac{NSW} experimental sample with non-missing prior earnings data who were randomly assigned to receive job training after December 1975, and completed it by January 1978; the sample with $d_i=0$ is a non-experimental sample of $2,490$ men taken from the PSID\@. The outcome $Y_{i}$ corresponds to earnings in 1978 in thousands of dollars, and the covariate vector $x_i$ contains the variables: age, education, indicators for black and Hispanic, indicator for marriage, earnings in 1974, earnings in 1975, and employment indicators for 1974 and 1975.\footnote{Following abadie_bias-corrected_2011, the no-degree indicator variable is dropped, and the employment indicators are defined as an indicator for nonzero earnings.} We assume that the unconfoundedness assumption holds given this covariate vector.

We are interested in the \ac{CATT},

The analysis in (ref) goes through essentially unchanged with this definition of $Lf$ (see (ref)).

Implementation details

To construct feasible versions of our estimators and \acp{CI}, we follow the baseline implementation in (ref). Implementing the procedure requires us to fix the norm $\norm{\cdot}_{\mathcal{X}}$ and the smoothness constant $C$ in the definition (ref). We consider weighted $\ell_{q}$ norms of the form

where $A$ is a diagonal matrix. Let us now describe our main specification; we discuss other specifications in the supplemental materials. To make the restrictions on $f$ implied by the choice of norm and $C$ interpretable, in our main specification, we set $C=1$, and $A=A_{\text{main}}$, with the diagonal elements of $A_{\text{main}}$ given in (ref). The $j$th diagonal element $A_{jj}$ then gives the a priori bound on the derivative of the regression function with respect to $x_{j}$ (i.e.\ the partial effect of increasing $x_{j}$ by one unit). We set $q=1$, so that the cumulative effect of changing multiple elements of $x_{j}$ by one unit is bounded by the sum of the corresponding elements $A_{jj}$.

The elements of $A_{\text{main}}$ are chosen to give restrictions on the conditional mean $f$ that are plausible when $C=1$; we report results for a range of choices of $C$ as a form of sensitivity analysis. While the outcome is measured in levels, it is easier to interpret the bounds in terms of percentage increase in expected earnings. As a benchmark, consider deviations from expected earnings when the expected earnings are \$10,000, that is $f(x_{i}, d_{i})=10$. Since the average earnings for the $d_{i}=1$ sample are \$6,400, with 78% of the treated sample reporting income below 10 thousand dollars, the implied percentage bounds for most people in the treated sample will be even more conservative than this benchmark. We set the coefficients on black, Hispanic, and married to $2.5$, implying that the wage gap due to race and marriage status is at most 25% at this benchmark. We set the coefficients on 1974 and 1975 earnings so that increasing earnings in each of the these years by $x$ units leads to at most $(0.5+0.5)x$ increase in 1978 earnings, that is at most a one-to-one increase. Including the employment indicators allows for a small discontinuous jump in addition for people with zero previous years' earnings. The implied bounds for the effect of age and education on expected earnings at 10 thousand dollars are 1.5% and 6%, respectively, in line with the 1980 census data.

With this choice of norm, $\norm{x}_{\mathcal{X}}=\norm{x}_{A_{\text{main}}, 1}=\sum_{j=1}^{p}\abs{A_{\text{main}, jj}x_{j}}$, we follow the baseline implementation in (ref). In step (ref), we use the homoskedastic estimate $\tilde \sigma^2(x, d)=\hat{\sigma}^{2}$. Here $\hat{\sigma}^{2}=\frac{1}{n}\sum_{i=1}^n \hat{u}_{i}^2$, where $\hat{u}^{2}_{i}$ is the nearest-neighbor variance estimator with $J=2$ neighbors, using Mahalanobis distance (using the metric $\norm{\cdot}_{A_{\text{main}},1}$ leads to very similar results). We compare our estimators and \acp{CI} to those based on matching estimators: here we again use the norm $\norm{\cdot}_{A_{\text{main}}, 1}$ to define distance. While all reported confidence intervals and standard errors use the heteroskedasticity-robust formula (ref) (using the nearest-neighbor estimate $\hat{u}_{i}^{2}$ above), when making efficiency comparisons, we use homoskedastic standard errors, so that \ac{RMSE} and \ac{CI} length are the same as those used to optimize the choice of the smoothing parameter $\delta$, or the number of matches $M$.

Results

(ref) reports the point estimates and \acp{CI} at values of $\delta$ and $M$ that optimize the \ac{RMSE} and \ac{CI} length criteria. There are three aspects of the results worth highlighting. First, in all cases, the worst-case bias is non-negligible relative to the standard error. This is consistent with (ref) on impossibility of $\sqrt{n}$-inference: under the Lipschitz smoothness assumption it is not possible to construct asymptotically unbiased estimators in this application given that the dimension of continuously distributed covariates is $4$. Our \acp{CI} reflect this by explicitly taking the bias into account. Second, in line with the predictions of (ref), the maximum Lindeberg weights (ref) are low for the feasible estimators $\tilde{L}_{\delta}$ in columns (1) and (2). However, due to imperfect overlap, for the matching estimator with $M=1$ match in column (3), it is considerably higher, and above the 0.075 cutoff suggested in NoRo19. Thus, its distribution may not be well-approximated by a normal distribution, in line with the discussion in (ref). Third, in Panel B, the table also reports \acp{CI} for the \ac{PATT}, constructed using the formula (ref), using the nearest neighbor variance estimator to estimate the marginal variance. These \acp{CI} are only slightly longer than those for the \ac{CATT}: depending on the column, the increase in length is 10.3--13.6%.

To examine sensitivity of the results to the specification of the parameter space $\mathcal{F}$, (ref) plots the estimator $\tilde{L}_{\delta}$ at the \ac{RMSE}-optimal choice of $\delta$, as well as \ac{CI} for the \ac{CATT} and \ac{PATT} when the smoothing parameter $\delta$ is chosen to optimize \ac{CI} length. For very small values of $C$---smaller than $0.1$---the Lipschitz assumption implies that selection on pretreatment variables does not lead to substantial bias, and the estimator and \acp{CI} incorporate this by tending toward the raw difference in means between treated and untreated individuals, which in this data set is negative. For $C\geq 0.2$, the point estimate is positive and remarkably stable as a function of $C$, ranging between $0.94$ and $1.14$, which suggests that the estimator and \acp{CI} are accounting for the possibility of selection bias by controlling for observables. The two-sided \acp{CI} become wider as $C$ increases, which, as can be seen from the figure, is due to greater potential bias resulting from a less restrictive parameter space.

According to (ref), matching with $M=1$ is efficient when $C$ is “large enough”. In our application, for $C\geq 2.8$, the efficiency of the matching estimator is at least 95% for both \ac{RMSE} and \ac{CI} length. Matching with $M=1$ leads to a modest efficiency loss in our main specification, where $C=1$: its efficiency is $90.4\%$ for \ac{RMSE}, and $86.0\%$ for the construction of two-sided \acp{CI}. However, inference results based on the matching estimator should be taken with a grain of caution due to the concerns with the accuracy of the \ac{CLT} approximation discussed above.

Comparison with experimental estimates

The present analysis follows, among other, LaLonde1986, dw99, smith_reconciling_2001, smith_does_2005 and abadie_bias-corrected_2011 in using a non-experimental sample to estimate treatment effects of the \ac{NSW} program. A major question in this literature has been whether the non-experimental sample can be used to obtain results that are in line with the estimates based on the original experimental sample of individuals who were randomized out of the \ac{NSW} program. In the experimental sample, the difference in means between the outcome for the treated and untreated individuals is $1.79$. Treating this estimator as an estimator of the \ac{CATT}, the (unconditional) robust standard error is $0.64$; treating it as an estimator of the \ac{PATT} (which also coincides with the \ac{PATE}), it is $0.67$.

The estimates in columns (1) and (2) of (ref) are slightly lower, although the difference between them and the experimental estimate is much smaller than the worst-case bias. Consequently, all the difference between the estimates can be explained by the bias alone. The large value of the worst-case bias also suggests that the goal of recovering the experimental estimates using the current non-experimental dataset is too ambitious, unless one is willing to impose substantially stronger smoothness assumptions. Furthermore, differences between the estimates reported here and the experimental estimate may also arise from (1) failure of the selection on observables assumption; and (2) the sampling error in the experimental and non-experimental estimates.