EconBase
← Back to paper

Nonparametric estimation of causal heterogeneity under high-dimensional confounding

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.

52,745 characters · 16 sections · 81 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.

\thispagestyle{empty}

center[center omitted — 1,485 chars of source]

JEL classification: C14, C21\\ Keywords: causal machine learning, effect heterogeneity, group average treatment effects, semiparametric efficiency, ensemble learning\\

\setcounter{footnote}{0} \setcounter{page}{1} \doublespacing \allowdisplaybreaks

Introduction

Recently, new machine learning based estimators showed immense potential to systematically uncovering causal effect heterogeneity so that there is now a rapidly growing literature on this topic (e.g., see the overviews in Athey_Imbens_2017,Athey_Imbens_2019 and Knaus_Lechner_Strittmatter_2018). In the context of heterogeneity, the respective aggregation levels for which the heterogeneity is estimated is playing an important role. Most papers of this literature focus on a selection-on-observable framework and investigate estimators for the heterogeneity at the lowest aggregation level to uncover possible heterogeneities to the largest extent possible. While this finest level of causal granularity is obviously of interest, Chernozhukov_FernandezVal_Luo_2018 and Lechner_2018 argue to analyse heterogeneity at higher levels, so called `Group Average Treatment Effects' (GATEs). Such aggregates can be estimated more precisely, may be far more easily interpretable by researchers in substantive terms, and are more useful for decision makers. In particular, some subgroup heterogeneities are of limited value per se because it is hard to justify a decision or policy based on certain characteristics (race, gender etc.). Therefore, decision-makers are often only interested in effect heterogeneities based on a rather small subset of available covariates. This paper suggests an approach that is based on statistical-learning assisted estimation of the GATEs for the various discrete and continuous variables of interest, and subsequent non-parametric aggregation of the GATEs to obtain `Average Treatment Effects' (ATEs).\footnote{Lechner_2018 also proposed this aggregation idea. However, that paper considered only a version of a Causal Forest while here we are in principle agnostic with respect to the machine learning method used. Furthermore, it considered only GATEs based on discrete variables and thus GATEs were obtained as unweighted within-cell means.}\\ More technically speaking, in effect heterogeneity analysis covariates do not (only) serve the purpose of making identifying assumptions credible. They become part of the outcome analysis by discriminating different subgroups of units for which the effect is of interest. Further, whenever new observations enter the sample the covariate realizations could be used to predict a causal effect. The set of covariates to be included in the statistical model to explore effect heterogeneity is therefore not a statistical but rather a substantive decision.\\ The estimation of subgroup specific effects is a tedious task when there is confounding. In such settings, causal effects are typically only identified if the researcher includes the confounding covariates in the statistical model as well. Hence, the identifying assumptions dictate the inclusion of the set of covariates required. In empirical research based on selection-on-observables, the credibility of causal effects estimation often depends on a very large set of possible covariates with very many possible functional forms. Qualitatively assessing which covariates should ultimately enter the model in which specific form in a non-systematic fashion is prone to be flawed.\\ There are currently several suggestions to estimate heterogeneous effects when there is confounding. The general concept of estimating causal effects conditional on covariates\footnote{We avoid the imprecise term `Conditional Average Treatment Effect' because it is unclear which conditioning set is actually meant.} already dates back to Hahn_1998. He suggested estimating a nonparametric outcome regression on the set of covariates that needs to be controlled for. Averaging over the conditional means leads to estimators of ATE that attain the semiparametric efficiency bound. In practice, however, nonparametric regression with many covariates is hardly feasible because the convergence rate of nonparametric methods exponentially decreases with the number of covariates included. Recently, Wager_Athey_2018 follow the same ideas as in Hahn_1998 but use Causal Forests instead of standard nonparametric regression. Athey_Tibshirani_Wager_2019 and Lechner_2018 modify the Random Forest algorithm to better adjust for confounding and improve precision. Outcome-based models adjust for confounding and infer heterogeneous effects in a single estimation step. Therefore, in all of these approaches inference for effect heterogeneity relies on a dimension of the covariate space that is fixed. Given the previous discussion, this might be a very strong assumption.\\ In this paper we follow an alternative approach in the literature. The two distinct roles of the covariates -- adjusting for confounding and estimating heterogeneous effects -- are explicitly reflected in a two-step estimation procedure for the GATEs. This idea is conceptionally not new in the literature.\footnote{A few days before this work appeared first on arXiv, Fan_Hsu_Lieli_Zhang_2019 published their independent work on arXiv (up to that moment unknown to us) that uses similar ideas about aggregation and machine learning.} In the context of difference-in-differences estimation, Abadie_2005 shows that propensity scores weighted outcomes can be used as a dependent variable in a second stage regression on the covariates that are of interest for heterogeneous effects. Abrevaya_Hsu_Lieli_2015 use a similar idea in the standard selection-on-observables setting. They provide inferential results for nonparametric and parametric propensity score first stages with nonparametric second stages. In line with other results in the literature on average effects (Hirano_Imbens_Ridder_2003, Robins_Rotnitzky_Zhao_1994, Lunceford_Davidian_2004), they show that the variance for Inverse Probability Weighting (IPW) estimators can be substantially decreased when the propensity score is estimated nonparametrically. Since their second stage also relies on nonparametric regression, the validity of their asymptotic results requires jointly choosing two kernel bandwidths which have to be in a rather small feasible interval. Lee_Okui_Whang_2016 augment the model of Abrevaya_Hsu_Lieli_2015 by including outcome projections (Augmented IPW, AIPW) and show that when both the propensity score and the outcome projections are estimated parametrically, one can treat the nuisance parameters as if they were known. Their asymptotic results with parametric first stages are then equivalent\footnote{While Abrevaya_Hsu_Lieli_2015 use local constant nonparametric regression, Lee_Okui_Whang_2016 show their results with local linear nonparametric regression.} to those of Abrevaya_Hsu_Lieli_2015 with nonparametric propensity score estimation. Still the parametric nuisance models used could lead to substantial misspecification bias if the functional form does not coincide with the unknown data generating process (DGP). Moreover, if there are more potential regressors than observations even parametric estimators collapse.\\ In their contributions to ATE estimation Belloni_Chernozhukov_Hansen_2014 and Chernozhukov_Chetverikov_Demirer_Duflo_Hansen_Newey_2017 show how AIPW type estimators can be adapted to settings where the dimension of the relevant confounders grows with the sample size. They use various machine learning methods to estimate the propensity score and the outcome projections. Recently, Chernozhukov_Semenova_2017 use their framework to estimate effect heterogeneity based on linear models. They provide conditions under which their second-stage linear model can become increasingly flexible.\\ Postulating nonparametric second stages, we do not assume any specific functional form of the GATE. We contribute to the literature on GATE estimation by allowing for flexible functional forms as well as a high-dimensional\footnote{In general the term `high-dimensional' refers to the fact that the dimension of the model can grow with the sample size. We will provide specific rate conditions in the main part of the paper.} confounder space. This enables our proposed estimator to be robust against functional form misspecification and to remain consistent even if the number of covariates relative to the sample size is large. In particular, we provide a generic statistical framework such that the convergence rate requirements of the first stage nuisance estimation are coupled with the kernel bandwidth second stage nonparametric convergence.\\ Additionally, we link our identification and estimation result to semiparametric efficiency theory by providing a new estimator for the ATE that can be estimated as a by-product of the GATEs. The estimator aggregates over all point estimates of the GATEs. We show that under certain convergence conditions for the kernel bandwidth, asymptotically it hits the variance lower bound of the semiparametric estimation problem. We therefore also contribute to the small literature on three-step semiparametric ATE estimation. Specifically, Hahn_Ridder_2013 (for an alternative theoretical development see also Mammen_Rothe_Schienle_2012) investigate a related set-up showing that nonparametric regression on an estimated propensity score can lead to efficient estimation of ATE. To the best of our knowledge, this paper is, however, the first that analyses the asymptotic properties of averaging a transformed outcome projection instead of the outcome projection on the propensity score. Like propensity score matching, our three-step estimator might have better finite-sample properties than two-step estimators that share the same first-order asymptotic properties (Robins_Rotnitzky_Zhao_1994, Hirano_Imbens_Ridder_2003, Chernozhukov_Chetverikov_Demirer_Duflo_Hansen_Newey_2017) because the propensity score weights are subject to an additional smoothing step. Unlike propensity score matching, IPW or Hahn_1998's (\citeyear*{Hahn_1998}) estimator, the proposed ATE estimator remains feasible when the dimension of the confounders entering the model is high.\\ After providing some more information on the theoretical background in Section (ref), we present the details of our main asymptotic results in Section (ref). An empirical example in Section (ref) compares different alternative estimators for GATE and ATE and illustrates the applicability and usefulness of the new methods. The last section concludes. The formal proofs of our theorems as well as some details on the empirical implementation are relegated to the Appendix.

Methodology

Notation

Suppose that we observe an independent and identically distributed random sample $\{w_i\}_{i=1}^N$ with sample size $N$ where $w_i=(y_i,d_i,x_i,z_i)$. Denote with uppercase letters a variable and with lowercase letters its realizations. Then $Y$ is the outcome variable and $D$ is the binary treatment of interest. To describe causal effects, we use Rubin_1974's (\citeyear*{Rubin_1974}) potential outcome notation such that $Y^d$ is the outcome that would have been observed under treatment $D=d$. Further, $X$ is a matrix of observed covariates with support $\mathcal{X}$ and $Z\subseteq X$ as a set of predefined variables where the researcher is interested in effect heterogeneity with support $\mathcal{Z}$. Also let $X\in\mathbb{R}^{\dim{X}}$ and $Z\in\mathbb{R}^{\dim{Z}}$ and denote $\lambda_X=\dim{X}$ and $\lambda_Z=\dim{Z}$. Potentially we have that $\lambda_X\rightarrow\infty$ when $N\rightarrow\infty$ whereas $\lambda_Z$ is fixed.\footnote{The concrete growth rates of $\lambda_X$ in relation to $N$ will be discussed in Section (ref).} Hence, we explicitly allow for models where the dimension of $X$ is high-dimensional but the dimension of the subset of covariates that is of interest for the heterogeneity analysis does not grow with the sample size. We remain agnostic about the underlying cumulative distribution from which the sample of $W=(Y,D,X,Z)$ is drawn $F=F(W)$ and just assume that it exists with density $f=f(W)$.

Semiparametric efficiency theory

The main parameter of interest in this study is the GATE defined as

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

Since we want to avoid usually unrealistic parametric assumptions on the underlying DGP, we allow for a flexible function $\psi(W,\cdot)$ such that GATE is identified as

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

using the variables observed or functions of them. Many possible transformations of the outcome exists (see the references mentioned in the Introduction). It is, however, a priori unclear which one is a `good' transformation in the sense that it achieves the variance lower bound for the problem. Hence, ideally our exposition would start by deriving the semiparametric efficiency bound for the problem at hand and then use the moment condition implied as an estimand for GATE. However, since parameters using `last stage' nonparametric projections are not pathwise differentiable, standard semiparametric efficiency bounds cannot be derived following established theory (e.g. Bickel_Klaassen_Ritov_Wellner_1993, Newey_1994, Hahn_1998, Tsiatis_2006). This was also noted for different problems in Rubin_vanderLaan_2007 and Kennedy_Ma_McHugh_Small_2017. We follow their approaches. Instead of directly relying on an efficiency result for GATE, we use that

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

implying Hahn_1998's (\citeyear*{Hahn_1998}) efficient score function for ATE

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

where $p(X)=\mathbb{E}\left[D|X\right]$ denotes the propensity score and $m_d(X)=\mathbb{E}\left[Y|X,D=d\right]$ for $d\in\{0,1\}$ denotes the conditional expectations of the outcome in the treatment-specific subpopulations.

Parameter identification

For identification of GATE and ATE we make the following assumptions.\\

assumption[Conditional independence] \begin{align*} Y^0,Y^1\perp D|X=x \quad \forall x\in\mathcal{X}\\ \end{align*}
assumption[Stable Unit Treatment Value Assumption (SUTVA)] \begin{align*} Y=DY^1+(1-D)Y^0\\ \end{align*}
assumption[Exogeneity of confounders] \begin{align*} X^1=X^0\\ \end{align*}
assumption[Common support] \begin{align*} c<p(X)<1-c \end{align*} for some small positive constant $c$.\\

Assuming that appropriate moments exist, then for GATE we have

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

The exposition shows that IPW and outcome based estimands are embedded in the estimand based on $\psi(W,p,m_0,m_1)$. Finally, by noticing that

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

identification of ATE trivially follows from these considerations.

Main results

GATE estimation

Proposed estimator

The identification results from the preceding section suggest a two-step estimation strategy. The details of our proposed estimator are described in Procedure 1.

figure[figure omitted — 1,291 chars of source]

In a first step a sample plug-in versions of $\psi(W,p,m_0,m_1)$ can be obtained by estimating the nuisance parameters. In a second step the $\psi$-vector can be projected on $Z$. Our goal is to estimate both stages as flexible as possible and to avoid parametric assumptions. Further, our estimator can cope with settings where $\lambda_X$ is very large which precludes classical nonparametric and parametric methods to estimate the first stage nuisances $p(x)$, $m_0(x)$ and $m_1(x)$. However, we can use a large class of supervised machine learning algorithms that have been shown to be very effective predictors for such types of tasks. Following the suggestions of Chernozhukov_Chetverikov_Demirer_Duflo_Hansen_Newey_2017 we apply a cross-fitting algorithm for the nuisance parameter estimation step in order to guarantee that the resulting estimator of $\psi(W,p,m_0,m_1)$ consists of independent observations. The requirements for the second stage estimation step are more sophisticated as this estimator should allow for valid inference. To estimate GATE flexibly, we apply nonparametric local constant regression in the second step.

Asymptotic results

We now investigate the theoretical properties of our proposed estimation procedure. To ease the notational burden, we start with some definitions.\\

definition[Norms] Denote by $\lVert g(X)\rVert_p$ the $L_p$ norm of the generic function $g(\cdot)$. Further denote the supremum norm by $\sup_{X\in\mathcal{X}}\lvert g(X)\rvert=\lVert g(X)\rVert_{\infty}$.\\
definition[Rates] The nuisance parameter first stage estimates $\hat{p}$, $\hat{m}_0$ and $\hat{m}_1$ obtained by the sample splitting procedure described above belong to the realization sets $\mathcal{P}$, $\mathcal{M}_0$ and $\mathcal{M}_1$ with probability $1-o(1)$. For any realization $p^*$, $m_0^*$ and $m_1^*$ in the sets define the rates \begin{align*} \epsilon_{m_d}&=\sup_{m_d^*\in\mathcal{M}_d}\lVert m_d^*(X)-m_d(X)\rVert_2\\ \epsilon_{p}&=\sup_{p^*\in\mathcal{P}}\lVert p^*(X)-p(X)\rVert_2\\ \epsilon_{\max}&=\max\{\epsilon_{m_0},\epsilon_{m_1},\epsilon_p\}.\\ \end{align*}
definition[Scaling factor] For any function $g$ define a scaling parameter $\delta_g$ that determines $g=O\left(N^{-\delta_g}\right)$.\\

We then make the following standard assumptions on the kernel regression step (see for example Pagan_Ullah_1999).\\

assumption[Kernel regression] \begin{enumerate} • $Z=z$ is a point in the interior of the support $\mathcal{Z}$. • The density function estimator is uniformly bounded away from zero such that $\inf_{z\in\mathcal{Z}}\hat{f}(z)\geq C$ where $C>0$ is a generic constant. • The Kernel function $K(u)$ is $r$ times continuously differentiable, symmetric and of order $r$ in the sense $\int u^{r-1}K(u)du=0$ and $\int u^{r}K(u)du=O(1)$ for $r\in\mathbb{N}$. • $f(z)$ and $\mathbb{E}\left(\psi(W,p,m_0,m_1)|Z=z\right)$ are $r$ times continuously differentiable. • Further the Kernel function satisfies (i) $\int K(u)du=1$, (ii) $\int \lvert K(u)\rvert^{2+C} du=O(1)$ for any $C>0$, (iii) $\lvert u\rvert\lvert K(u)\rvert\rightarrow 0$ as $\lvert u \rvert\rightarrow\infty$, (iv) $\left\lVert K(u)\right\rVert_{\infty}=O(1)$ and (v) $\int K^2(u)du=O(1)$.\\ \end{enumerate}

Assumption (ref) comprises the standard nonparametric local constant regression assumptions allowing for multiple covariates and higher-order kernels. For illustrative purposes multivariate regression results are derived assuming the same bandwidth for every regressor. Further, we have to make boundedness assumptions on the second moment of the sample error of the outcome model and on the nuisance prediction errors.\\

assumption[Boundedness of conditional variances] The conditional variances of the outcome models are bounded such that they obey \begin{align*} \mathbb{E}\left[\left(DY-m_1(X)\right)^2|X\right]=O(1)\quad and \quad \mathbb{E}\left[\left((1-D)Y-m_0(X)\right)^2|X\right]=O(1).\\ \end{align*}
assumption[Boundedness of convergence rates] The nuisance parameter prediction errors are bounded such that they obey \begin{align*} \sup_{m_d^*\in\mathcal{M}_d}\lVert m_d^*(X)-m_d(X)\rVert_\infty=O(1)\quad for $d\in\{0,1\}$ and \quad \sup_{p^*\in\mathcal{P}}\lVert p^*(X)-p(X)\rVert_\infty=O(1).\\ \end{align*}

Additionally, the convergence rates of our first stage nuisance parameter prediction and the second stage nonparametric regression are assumed to be as follows:\\

assumption[Coupled convergence (GATE)] The bandwidth $h$ and the sample size $N$ jointly converge such that \begin{itemize} • $h=o(1)$, $Nh^{\lambda_Z}\rightarrow\infty$ as $N\rightarrow\infty$ and • $N^{\frac{1}{2}}h^{\frac{1}{2}\lambda_Z}h^r=o(1)$. \end{itemize} Further, $N$ and $h$ satisfy the joint convergence conditions with the nuisance parameter convergence rates \begin{itemize} • $h^{-\frac{1}{2}\lambda_Z}\epsilon_{\max}=o(1)$$N^{\frac{1}{2}}h^{-\frac{1}{2}\lambda_Z}\epsilon_{m_0}\epsilon_p+N^{\frac{1}{2}}h^{-\frac{1}{2}\lambda_Z}\epsilon_{m_1}\epsilon_p=o(1)$.\\ \end{itemize}

Assumption (ref) comprises the coupled convergence rate assumptions that are at the centre of our theoretical results. Conditions (i) and (ii) quantify how the bandwidth has to converge to zero in relation to the sample size $N$ and the number of regressors $\lambda_Z$. As usual, the bandwidth has to go to zero but slower than the sample size grows to infinity. Also the bandwidth has to be chosen such that the asymptotic bias term vanishes faster than the variance. This allows to apply the Central Limit Theorem and makes the estimator asymptotically unbiased. In particular, condition (ii) requires undersmoothing in the sense that the bandwidth has to be below the mean squared error (MSE) optimal rate. As discussed for example in Pagan_Ullah_1999, choosing a higher order kernel mitigates the problem.\\ Conditions (iii) and (iv) state that $h$ has to be chosen such that first stage convergence rates vanish fast enough. In particular, by condition (iv) the joint convergence rates from propensity score and outcome projection estimation have to vanish faster than $\sqrt{N}$ scaled with the kernel bandwidth. Since $h\rightarrow 0$, this is a more restrictive assumption compared to the rate conditions usually required for average effects estimation (see for example Chernozhukov_Chetverikov_Demirer_Duflo_Hansen_Newey_2017). In contrast to average effects, one is only interested in estimating the effect at a prespecified point $Z=z$. Thus, since sample observations enter the estimator in a weighted form, the prediction precision needed is for the lower effective sample size around $Z=z$. Hence, the first stage prediction guarantees need to adapt to this smaller sample conditions and therefore achieve a faster joint rate of convergence in terms of the sample size $N$. Condition (iii) additionally prevents the worst rate from becoming arbitrarily slow especially when $\lambda_Z$ is larger than one. Still, our estimator has a `rate' double robustness feature in the sense that joint rates can vanish relatively slowly but all single first stage rates are restricted from converging very slowly. $L_2$ convergence rates of many supervised machine learning methods satisfy these properties under sparsity conditions. For example Belloni_Chernozhukov_2013 show that the predictive error of the Lasso is of order $O\left(\sqrt{\frac{s\log\max(\lambda_X,N)}{N}}\right)$ where $s$ the unknown number of true coefficients in the oracle model. Suppose that $s$ and $\lambda_X$ are equal in the outcome and the propensity score models then we require $\frac{s^2\log^2\max(\lambda_X,N)}{Nh^{\lambda_Z}}\rightarrow 0$. Hence the dimension of the confounding variables $\lambda_X$ can grow with the effective sample size $Nh^{\lambda_Z}$. Similar rates can be shown for $L_2$ boosting (Luo_Spindler_2016) and nonlinear models like Random Forests (Wager_Walther_2015) or forms of Deep Neural Nets (Farrell_Liang_Misra_2018).\footnote{For the concrete dependence of sparsity conditions on the parameters of the predictors see the references mentioned.}\\ A natural question is then if a bandwidth exists that satisfies the rate conditions in Assumption (ref). Indeed, one can show (for more details see Appendix (ref)) that the theoretical range of possible bandwidth choices can be described by

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

and we achieve a condition for the order of the kernel

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

For example if we restrict ourselves on second order kernel functions then for $\lambda_Z=1$ we require $\delta_{\epsilon_p}+\delta_{\epsilon_{m_d}}=\frac{3}{5}$. Similarly, for $\lambda_Z=2$ and $\lambda_Z=3$, $\delta_{\epsilon_p}+\delta_{\epsilon_{m_d}}=\frac{2}{3}$ and $\delta_{\epsilon_p}+\delta_{\epsilon_{m_d}}=\frac{5}{7}$ are required respectively. Thus, for a growing dimension $\lambda_Z$ the joint rate condition for the first stage nuisance parameters approaches the parametric rate.\\ The discussion indicates that given one has chosen an appropriate order of the kernel function, the researcher can choose the bandwidth somehow below but not too much below the MSE optimal rate. One could therefore simply use a certain fraction (e.g. 0.9) of the cross-validation bandwidth choice. Hence, from a practical perspective our bandwidth choice problem is equivalent to nonparametric regression with undersmoothing.\\ Given these assumptions we can then derive the first main theoretical result.\\

theoremUnder Assumptions (ref)-(ref) our proposed estimation procedure for GATE obeys \begin{align*} \sqrt{Nh^{\lambda_Z}}\left(\hat{\tau}-\tau\right)=\frac{1}{\sqrt{Nh^{\lambda_Z}}}\sum_{i=1}^N\frac{K\left(\frac{z_i-z}{h}\right)}{\frac{1}{Nh^{\lambda_Z}}\sum_{i=1}^NK\left(\frac{z_i-z}{h}\right)}\left(\psi(W_i,p,m_0,m_1)-\tau\right)+o(1) \end{align*} and \begin{align*} \sqrt{Nh^{\lambda_Z}}\left(\hat{\tau}-\tau\right)\rightarrow_d N(0,\sigma_{GATE}^2) \end{align*} with $\sigma_{\text{GATE}}^2=\frac{\int K(u)^2du\times\mathbb{E}\left[\left(\psi(W_i,p,m_0,m_1)-\tau\right)^2|Z=z\right]}{f(z)}$.\\

Theorem (ref) shows that under the assumptions discussed above the speed of convergence is determined only by the nonparametric regression step. In particular, it does not depend on the first stage estimation steps. An equivalent result can also be achieved by using IPW with nonparametric first stages (see Abrevaya_Hsu_Lieli_2015). However, this requires an additional bandwidth choice for the first stage propensity score regression and is limited to the case when also $\lambda_X$ is very small. The dimension of $X$ can be increased under functional form assumptions for the first stage. However, as shown by the authors the price to pay is an increase in the asymptotic variance. This is not the case for the estimator proposed in this paper. Also Theorem (ref) is valid under generally weaker conditions compared to the results in Lee_Okui_Whang_2016 for parametric first stages. Heuristically\footnote{In contrast, to Lee_Okui_Whang_2016 we use local constant instead of local linear regression and introduce cross-fitting for nuisance parameter estimation. This should, however, not be a concern for the intuitive argument made.}, if the first stage estimators converge at $\sqrt{N}$ then our conditions on the bandwidth are satisfied and our asymptotic results continue to apply. To this extent, our results comprise the result of Lee_Okui_Whang_2016 as a special case.

Joint estimation of GATE and ATE

Proposed estimator

Given the considerations so far, it appears `naturally' to estimate ATE in three steps as an average of GATEs in the sample. The details of the proposed estimator are described in Procedure 2.

figure[figure omitted — 392 chars of source]

As suggested in Chernozhukov_Chetverikov_Demirer_Duflo_Hansen_Newey_2017 one could also directly estimate ATE as the average of the vector $\psi(W,p,m_0,m_1)$ using the first stage nuisance parameter predictions. However, as a by-product of GATE estimation, using an additional kernel smoothing step may lead to an ATE estimator with better finite sample properties. In particular, the propensity score weights do not enter the last step of our estimator directly and our hope is that small misspecification errors of propensity scores close to zero or one are therefore smoothed out. The sensitivity of estimators incorporating inverse propensity score weights directly is the subject of many Monte Carlo experiments (e.g. Huber_Lechner_Wunsch_2013, and Froelich_2004b). We notice that similar reasoning is also behind three-step estimators that apply nonparametric regression on an estimated propensity score often used in practice.

Asymptotic results

To obtain our theoretical results we have to modify Assumption (ref) slightly.\\

assumptionprime[Coupled convergence (ATE)] The bandwidth $h$ and the sample size $N$ jointly converge such that \begin{itemize} • $h=o(1)$, $Nh^{\lambda_Z}\rightarrow\infty$ as $N\rightarrow\infty$, • $N^{\frac{1}{2}}h^{\frac{1}{2}\lambda_Z}h^r=o(1)$ and • $Nh^{4r}=o(1)$ and $Nh^{2\lambda_Z}\rightarrow\infty$. \end{itemize} Further, $N$ and $h$ satisfy the joint convergence conditions with the nuisance parameter convergence rates \begin{itemize} • $h^{-\lambda_Z}\epsilon_{\max}=o(1)$ and • $N^{\frac{1}{2}}h^{-\lambda_Z}\epsilon_{m_0}\epsilon_p+N^{\frac{1}{2}}h^{-\lambda_Z}\epsilon_{m_1}\epsilon_p=o(1)$.\\ \end{itemize}

We notice that averaging over the estimated projection $\mathbb{E}\left[\psi(W,p,m_0,m_1)|X\right]$ is a partial mean problem in the sense of Newey_1994b, Newey_1994. While parts (i) and (ii) of Assumption (ref) remain unchanged, the additional condition (iii) is necessary in order to guarantee that the MSE of the kernel regression estimator scaled with $N^{\frac{1}{4}}$ converges to zero. In this way we guarantee the applicability of Newey_1994's (\citeyear*{Newey_1994}) framework. We could have also assumed uniform convergence rates for the kernel regression step. However, this would involve a unnecessarily strong condition (for a discussion see Newey_1994 and also Newey_McFadden_1994).\\ Since we want ATE to converge with a rate of $\sqrt{N}$, the requirements on the first stage convergence rates in condition (v) are more restrictive than those in the respective condition of Assumption (ref). The range of theoretically feasible bandwidth choices reduces to

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

Assuming $\frac{1}{4r}<\frac{1}{\lambda_Z+2r}$ we get a modified condition for the order of the kernel function

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

In general this result indicates that one relies on a higher-order kernel function whenever $\lambda_Z>1$ when GATE and ATE are estimated jointly.\\ Under the stronger Assumption (ref) we can then derive the following efficiency result.\\

theoremUnder Assumptions (ref)-(ref) and (ref) and the regularity conditions on the nonparametric second step as in Newey_1994 our proposed estimation procedure for ATE has the influence function \begin{align*} \frac{D(Y-m_1(X))}{p(X)}-\frac{(1-D)(Y-m_0(X))}{1-p(X)}+m_1(X)-m_0(X)-\theta \end{align*} and therefore obeys \begin{align*} \sqrt{N}(\hat{\theta}-\theta)\rightarrow_d N(0,\sigma_{ATE}^2) \end{align*} where $\sigma_{ATE}^2$ is the semiparametric efficiency bound of Hahn_1998.\\

Conceptually, Theorem (ref) underpins our intuition from semiparametric theory outlined in Section (ref). The result shows that indeed every estimator that involves a nonparametric projection of the AIPW modified outcome on any low-dimensional subset of $X$ is consistent, asymptotically normal and achieves the semiparametric efficiency bound. This asymptotic result has also been shown for other estimators already discussed in Section (ref). In contrast to Hirano_Imbens_Ridder_2003's (\citeyear*{Hirano_Imbens_Ridder_2003}) estimator, Hahn_1998's (\citeyear*{Hahn_1998}) estimator and matching on the propensity score (Hahn_Ridder_2013), we do not rely on nonparametrically estimated first stages. Due to the fact that $\lambda_X\rightarrow\infty$ these estimators are of no practical use in our setting. Further, unlike AIPW with machine learning nuisance parameter estimation (Chernozhukov_Chetverikov_Demirer_Duflo_Hansen_Newey_2017), our estimator involves an additional step. Therefore the inverse propensity score does not directly enter our estimator but is smoothed through the additional nonparametric step. Asymptotically, this does not make any difference as the result in Theorem (ref) shows. However, in finite sample this could be a major advantage over the usual AIPW estimator.

Illustrative example

We investigate the applicability of our methods using Cattaneo_2010's (\citeyear*{Cattaneo_2010}) dataset on the effect of cigarette smoking on birthweight available from the Stata website.\footnote{The original dataset can be retrieved from \href{http://www.stata-press.com/data/r13/cattaneo2.dta}{here}.} The dataset contains the outcome variable birthweight in grams ($Y$), whether the mother smoked during pregnancy ($D=1$) and several covariates on the mother's health and socio-economic background ($X$). A detailed description of all covariates in the dataset can be found in Appendix (ref). Applied studies with different estimation approaches unambiguously find negative average effects (see Abrevaya_2006, daVeiga_Wilder_2008, Walker_Tekin_Wallace_2009). Conditional average treatment effects were investigated by Abrevaya_Hsu_Lieli_2015 and Lee_Okui_Whang_2016 who find that mother's age is associated with increasingly negative effects of smoking. We replicate their results and compare their estimators with ours. Clearly, this type of analysis is limited in its scope since the true DGP remains unknown. However, the dataset has the particular advantage that some strong hypothesis about the estimation results are plausible. (i) The effect of smoking on birthweight should be either negative or zero. (ii) The effect should be increasingly negative with mother's age.\\ As a second example we consider how the effect changes with the number of prenatal care visits. On the one hand a very low number of care visits could indicate the mother's insufficient access to medical infrastructure and therefore could be associated with particularly negative effects. On the other hand a very high number of care visits could indicate a poor health situation. Hence, it is a priori unclear how the treatment effect and health care visits are exactly related.

Empirical results

figure[figure omitted — 941 chars of source]
figure[figure omitted — 1,506 chars of source]

Figure (ref) depicts the main results of our empirical analysis. We estimate the GATEs as described in Procedure 1 using an ensemble learner comprising Lasso, Ridge, Elastic Net and a Random Forest. The weights of the ensemble are obtained by cross-validating the out-of-sample MSE of the procedure. $X$ in our specification is an extended variable set (`alldata') and is exactly documented in Appendix (ref). For example in contrast to Lee_Okui_Whang_2016 we also include the available characteristics for the father of the child, since they could be a good predictor for the smoking behaviour of the mother. The covariates enter our model very flexibly. For the penalized regression predictors we allow for polynomials up to order four and all two way interactions. The Random Forest has the particular advantage of being an ensemble of trees itself and is therefore very flexible by construction. The results are generally in line with the hypothesis made. In particular, the effect of smoking is unambiguously negative over the whole support of mother's age and prenatal care visits. As expected the effect increases with age. Interestingly, a higher number of prenatal care visits seems to be associated with higher negative effects.\\ We estimate all our results with second-order Gaussian kernel functions. The same analysis using higher order Gaussian kernel functions as proposed by Li_Racine_2007 yields similar results (see Appendix (ref)). In practice the biggest challenge is to determine the bandwidth for the nonparametric regression. To achieve undersmoothing, we multiply the bandwidth obtained by leave-one-out cross-validation with 0.9. Since this choice is arbitrary, the stability of our results towards this choice is a particular concern. Figure (ref) shows that our estimator is relatively robust regarding this choice. A major change in the shape of the function only appears for massive oversmoothing. An equivalent analysis for prenatal care visits yielding the same conclusion is relegated to Appendix (ref).

table[table omitted — 1,012 chars of source]

\\ Finally, Table (ref) shows the results for ATE estimation as described in Procedure 2. In line with the previous literature mentioned above, the average effect of smoking is estimated to be negative. Crucially, the estimated effect turns out to be very robust regarding the choice of the smoothing variable.

Comparison with other estimators

figure[figure omitted — 1,091 chars of source]

A `fair' comparison with other estimators is hardly feasible because our approach does not require specific functional form assumptions.\footnote{We do not consider nonparametric propensity score estimation as suggested in Abrevaya_Hsu_Lieli_2015 because it most likely does not allow to include all potential confounders in order to make Assumption (ref) credible.} In other words, the related estimators of Abrevaya_Hsu_Lieli_2015 and Lee_Okui_Whang_2016 suppose that they know the true propensity score or outcome projection specifications. Since we cannot compare our estimator against every possible parametric specification, we use the specification selected by Lee_Okui_Whang_2016 as a benchmark. Figure (ref) depicts GATE estimation results using the benchmark models. Strikingly, the IPW based estimator gives implausible results. For mother's age positive effects of smoking can almost nowhere be excluded. For care visits we do not obtain a GATE estimation result for our bandwidth choice. In fact, ATE is estimated indicating that the bandwidth is too large in order to obtain GATE estimates. In line with the theoretical result, standard errors are inflated compared to our estimation procedure. The AIPW based estimator with parametric models for the propensity score and the outcome projections gives plausible results for GATE with respect to mother's age. Slight differences arise when comparing the GATE curves regarding the number of prenatal care visits.

table[table omitted — 1,015 chars of source]

\\ Table (ref) shows the results for ATE estimation. As expected the results of Procedure 2 in Table (ref) are roughly in line with the standard AIPW based ATE estimator with ensemble first stages.\footnote{This might also be seen as a implicit test for the credibility of the stronger conditions required for the smoothed estimator compared to the averaged efficient score as in Chernozhukov_Chetverikov_Demirer_Duflo_Hansen_Newey_2017.} The relative bad performance of IPW based estimation for the GATEs is also reflected in the estimation of ATE. In particular, the standard error nearly doubles compared to AIPW based estimators and point estimates are reduced. Interestingly, for average effects there seems to be only little value-added for the flexible machine learning based estimators compared to the parametric specification.

Conclusion

In this study we propose new estimators for specific conditional and average causal effects when the dimension of the covariate space is high. In particular, by discriminating the different roles of covariates (adjusting for confounding vs. measuring causal heterogeneity of interest) in our approach, they can be included very flexibly -- not relying on any functional form assumptions. Rather, we show coupled convergence conditions for the different steps involved. The procedures suggested are based on semiparametric efficiency theory. In this sense, our proposed three-step estimator for ATE estimation is shown to reach the semiparametric efficiency bound. A widely used empirical example shows that our estimators are useful in practice. Compared to other estimators their desirable theoretical properties and increased flexibility could lead to divergent empirical results.\\ The specific structure of the GATE problem should be easily applicable to related settings. For example efficient score based estimation can also be used for instrumental variables problems (Chernozhukov_Chetverikov_Demirer_Duflo_Hansen_Newey_2017), difference-in-differences estimation (Zimmert_2018) and continuous treatment settings (Kennedy_Ma_McHugh_Small_2017).\\ Some other interesting problems and refinements are beyond the scope of this study and have to be left for further research as well. For example, the nonparametric regression estimator could be refined to the extent that its bandwidth is chosen in a data-adaptive manner. As an alternative to classical nonparametric regression, one could also investigate using methods from the toolbox of supervised machine learning. This might help to get reliable estimators even for cases when the dimension of $Z$ is moderately higher than considered in this paper, while sacrificing only little flexibility.\\ Finally, it might be worth to investigate the finite sample properties of the proposed three-step estimators for ATE compared to averaging the efficient score vector directly. Here, we consider our ATE estimator as a by-product of the GATE procedure underpinning the theoretical motivation of our framework. While the smoothed three-step estimator is first order asymptotically equivalent to directly averaging the efficient score vector, it might posses better finite sample properties since it does not directly rely on propensity score weights. However, the finite sample performance may crucially rely on the bandwidth choice and the set of covariates in $Z$. We regard this as yet another interesting direction for further research. \printbibliography