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.
83,235 characters · 10 sections · 77 citation commands
Estimation of Conditional Average Treatment Effects with High-Dimensional Data
In settings with individual-level treatment effect heterogeneity, the unconfoundedness assumption theoretically permits identification and consistent estimation of the conditional average treatment effect (CATE) for all possible values of the set of covariates $X$ used in adjusting for selection bias. One way to think about these covariates is that they are ex-ante predictors of an individual's potential outcomes with and without treatment, and hence are highly correlated with the treatment participation decision as well. Unconfoundedness states that the econometrician observes all relevant predictors so that conditional on $X$, the treatment takeup decision is no longer statistically related to the potential outcomes.\footnote{This condition was formalized by RR83; since then, unconfoundedness (or “selection on observables” or “conditional independence”) has become one of the standard paradigms for modeling selection effects. See, e.g., IW09 for further discussion.} Nevertheless, in many situations the individual deciding on treatment participation is likely to have access to private signals about their potential outcomes. Relying on the unconfoundedness assumption amounts to hoping that a set of publicly observed characteristics can still proxy for the information content of these signals. Therefore, the unconfoundedness assumption is more plausible in applications in which $X$ is a rich, detailed set of covariates, i.e., the dimension of $X$ is high.
While CATE as a function of $X$ provides a detailed characterization of treatment effect heterogeneity across observable subpopulations, this information is very hard to analyze and convey if $X$ is high-dimensional. Of course, one could examine slices of this function along some component(s) $X_1$ of $X$ while holding the other components $X_{-1}$ of $X$ constant. Nevertheless, how CATE varies as a function of $X_1$ will generally depend on the level at which $X_{-1}$ is held constant, requiring the examination of (infinitely) many different slices. For this reason, instead of holding the variables in $X_{-1}$ constant, AHL15 (henceforth AHL) suggest integrating them out with respect to the conditional distribution of $X_{-1}$ given $X_1$ or, in practice, a smoothed estimate of this distribution. This gives rise to a reduced dimensional CATE function that is easier to present and interpret.\footnote{If all covariates $X$ are integrated out, one obtains an estimator of ATE as in H98 or HIR03.}
In this paper we propose two-step estimators of the reduced dimensional CATE function where in the first step the required high-dimensional nuisance regressions are conducted by machine learning methods designed specifically to handle such problems, while the second integration step is implemented by a traditional local linear nonparametric regression.\footnote{This step assumes that $X_1$ is a continuous variable, which is the technically challenging and interesting case.} We derive the statistical properties of two variants of the estimator. In the first case, the first step (nuisance function estimation) and the second step (local linear regression) are both implemented over the full sample of available observations. In the second case, the available sample is split into parts, and the first step is implemented in one subsample while the second step is done in the complement sample. The roles of the subsamples are then rotated and the results are averaged. This is the “cross-fitting” approach to machine-learning-aided causal inference advocated by CC18. The first approach is used by BCFH13 in estimating unconditional treatment effects.
In proposing and studying these estimators, we contribute to two recent strands of the econometrics literature. First, we advance the currently available flexible methods for the estimation of reduced dimensional CATE functions due to AHL and LOW16 (henceforth LOW). Second, we make technical contributions to the recent literature that employs machine learning methods in tackling the prediction component of causal inference problems; see, e.g., BCH14b,BCH14a, BCFH13, CC17, CC18. Taking a broader perspective, our paper is also related to a large statistics literature on regular estimation and the use of orthogonal (doubly robust) moment conditions.
Regarding the first set of papers, AHL use an inverse probability weighted conditional moment of the data to identify CATE. They consider both kernel-based and parametric estimation of the propensity score in the first step, and derive the asymptotic distribution of the estimated CATE function evaluated at a fixed point $x_1$ in the support of $X_1$. LOW advance these results in two respects: their estimator is based on a Neyman-orthogonal moment condition and they also provide a method for uniform inference about the CATE function as a whole (rather than point by point). While LOW only use parametric models to estimate the nuisance functions involved in the moment condition, orthogonality lends their estimator a “double robustness” property: either the model for the propensity score or the models for the conditional means of the potential outcomes are allowed to be misspecified (but not both).
The CATE estimator proposed here is based on the same orthogonal moment condition as in LOW, but the required nuisance functions are estimated by machine learning methods, which allow for data-driven flexible functional forms as well as a (very) high-dimensional set of covariates. Neyman orthogonality is crucial in ensuring that the proposed CATE estimators are robust to the regularization bias inherent in the first stage, making post-selection inference possible. As the asymptotic theory is derived from high-level assumptions, there are a number of applicable first-stage estimation methods in practice, such as a random forest or $\ell_1$-penalized lasso or post-lasso. In this paper we use lasso estimation as the leading example.
In light of the discussion of the unconfoundedness assumption above, replacing the parametric estimators in LOW with machine learning methods greatly enhances the applicability and empirical relevance of flexible CATE estimation. At the same time, the asymptotic theory remains tractable: we provide methods for pointwise as well as uniform inference about the CATE function under both the full sample and sample-splitting implementation schemes. The uniform methods utilize the multiplier bootstrap, while pointwise inference can be based either on the bootstrap or the analytic results.
Turning to the literature on machine learning in treatment effect estimation, we build primarily on BCFH13 for the full-sample method and CC17,CC18 for the split-sample method, while providing the necessary extension of the theory to account for the use of local linear regression in second step. In these papers the parameter of interest is identified by the restriction that the unconditional expectation of a “score function” evaluated at the true parameter value (and the true nuisance functions) is zero. By contrast, the identifying restriction in our case is that the conditional expectation of the same score function is zero. Hence, our estimation procedure does not simply consist of substituting in the estimated nuisance functions and setting the sample average score to zero; instead, the score function will enter a local linear regression with kernel weight $\mathcal{K}((X_{1i}-x_1)/h)/h$ on each observation, where $h$ denotes a smoothing parameter (bandwidth).
The key high-level assumptions we employ in deriving our asymptotic results involve bounding the $L_\infty$ norm of the difference between the true and estimated nuisance functions, and the $L_2$ norm of the same difference multiplied by the kernel. The rates at which these error bounds are required to converge to zero are closely linked to the rate at which the bandwidth sequence converges to zero. From a purely technical standpoint, incorporating the bandwidth conditions into the high-level norm bounds in the full-sample as well as the cross-fitting case is a central contribution of the paper. Similarly to AHL and LOW, the resulting convergence rate of the CATE estimators is $\sqrt{Nh^d}$, where $N$ is the sample size and $d=dim(X_1)$.
In addition to the error bounds, the full-sample estimator also requires controlling the complexity (entropy) of the function space in which the nuisance functions take values. In the case of lasso estimation, this can be accomplished by restricting how fast the number of covariates and the sparsity indices associated with the nuisance functions are allowed to increase with the sample size. These conditions are more stringent than in the case of estimating ATE.
There are several papers in the broader statistics literature that have considered estimation problems related to ours R04, L13, Luedtke2016a, Luedtke2016b, NW17, L19.\footnote{We thank Edward Kennedy and an anonymous referee for these references.} NW17, in particular, estimate the full-dimensional CATE function in a data-rich environment using penalized regression, and establish the quasi-oracle error bounds for their estimator. While we also use a high-dimensional set of covariates and machine learning methods to deal with selection into treatment, the ultimate parameter of interest, being a function of a low-dimensional subset of the covariates, is then targeted by a traditional nonparametric estimator. We also complement NW17 by establishing both pointwise and uniform inference procedures. In a related paper, ZL19 consider the local constant estimation of CATE in the high-dimensional setting but only provide pointwise asymptotic results.
Another closely related paper, CS19, proposes an approach to CATE estimation that also includes a dimension-reduction step. There are, however, substantial technical differences between their paper and ours. First, the traditional nonparametric estimator used by CS19 in the second stage is series regression rather than local linear regression. Second, they only consider the cross-fitting approach and do not address the problem of estimating both the nuisance functions and the target function on the full sample. Third, we also provide a reasonably detailed discussion of the primitive conditions under which lasso estimation fulfills the high-level conditions posited in the paper, while CS19 restrict attention to high-level analysis.
Finally, our estimation method based on doubly robust moments is tied to the classic literature on regular estimation and semiparametric efficiency B83, P90, B93, N94, V00. As mentioned above, in a parametric setting, estimation of CATE based on doubly robust moments is consistent as long as either the treatment assignment process or the outcome processes is correctly specified. If both processes are nonparametrically estimated, the method can achieve a faster convergence rate than the nuisance estimators employed. The use of doubly robust methods for causal inference has also been considered by RR95, H98, VR03, HIR03, V06, f07, t07, VR11, BCFH13, F15, K17, R17, WA18, and suz19, among others.
In addition to providing theoretical results, we study and illustrate our methods through Monte Carlo simulations. The proposed estimators perform well in terms of bias, MSE, and coverage rates. In general, we find that the cross-fitting estimator has somewhat better finite sample properties than the full sample estimator, and thus we suggest using the cross-fitting estimator in empirical studies with reasonably large sample sizes.
Our application uses vital statistics data from North Carolina to estimate the effect of a (first-time) mother's smoking during pregnancy on the baby's birth weight as a function of the mother's age. Despite a number of previous analyses, the application is well worth revisiting with the help of machine learning methods, as there are a large number of covariates describing the mother's characteristics and events during pregnancy, and the specification of the propensity score is known to have a substantial impact on the results (see AHL, Section 4.2). Our results provide some corroborating evidence that the negative effect of smoking on birth weight becomes more detrimental with age. This pattern is less prevalent than some of the results reported in AHL but stronger than that found by LOW.
The rest of the paper proceeds as follows. In Section 2 we describe the formal setup and the estimators. Section 3 states and discusses the assumptions underlying the first-order asymptotic theory and provides the main results. Section 4 describes how to conduct uniform inference using the multiplier bootstrap. The application is presented in Section 5, while Section 6 concludes. An online supplement contains additional empirical studies, the Monte Carlo exercise as well as detailed proofs of the theoretical results.
Population units are characterized by a random vector $(D, Y(1),Y(0), X)$, where $D\in\{0,1\}$ indicates the receipt of a binary treatment, $Y(1)$ and $Y(0)$ are the potential outcomes with and without the treatment, respectively, and $X$ is a vector of pre-treatment covariates. The observed variables are given by the vector $W=(D, Y, X)$, where $Y=DY(1)+(1-D)Y(0)$. The distribution of $(D, Y(1),Y(0), X)$, and hence $W$, is induced by an underlying probability measure $\mathbb{P}$; parameter values computed under $\mathbb{P}$ will be denoted by the subscript “0” and represent the true values of these parameters. The expectation operator corresponding to $\mathbb{P}$ is denoted by $\mathbb{E}$, but we also use the linear functional notation $\mathbb{P} f:=\int f(w)d\mathbb{P}=\mathbb{E}[f(W)]$.
In order to accommodate high-dimensional data, we follow the conceptual considerations in F15 and treat the DGP (the measure $\mathbb{P}$) as dependent on the sample size $N$, allowing, in particular, the dimension of $X$ to grow with $N$.\footnote{This implies that the nuisance functions $\mu_0(j,X)$, $j=0,1$ and $\pi_0(X)$, to be defined below, may generally depend on $N$ as well.} This has two practical interpretations. First, the number of raw controls can already be comparable to the sample size or, second, $X$ may be composed of a large dictionary of sieve bases derived from a fixed dimensional vector $X^*$ through suitable transformations (e.g., powers and interactions). Thus, the high dimensionality of $X$ can also stem from the desire to provide a flexible approximation to the required nuisance functions. We explicitly allow for the use of lasso-type methods in the first stage that select a smaller subset of terms from the dictionary to approximate these functions.
To ease the already heavy notational burden in the paper, the dependence of the DGP on the sample size is left implicit throughout, but is of course accounted for in the theoretical analysis. Most arguments in the paper are based on concentration inequalities, which are non-asymptotic in nature. In our Assumption (ref) below, we also take into account the fact that, as the dimension of $X$ grows, the complexity of the first-stage estimator will generally diverge, which can affect the rate of convergence of our second-stage estimator. Furthermore, we establish the uniform inference results using the multiplier bootstrap based on the strong approximation theory developed by CCK14, which does not require the existence of an asymptotic distribution.
Given a $d$-dimensional subvector $X_1\subset X$ composed of continuous variables, the reduced dimensional CATE function is defined as \[ \tau_0(x_1)=CATE(x_1)=\mathbb{E}[Y(1)-Y(0)|X_1=x_1].\footnote{The most relevant case in practice is $d=1$ or perhaps $d=2$, otherwise the motivating properties of the reduced dimensional CATE function (interpretabilty and presentability) are lost. As we will see below, under a fourth-moment condition on $Y$, the general theory requires $d\le 3$. There are no restrictions on $d$ for bounded outcomes.} \]
The identification of $\tau_0(x_1)$ from the joint distribution of $W$ is facilitated by the unconfoundedness assumption along with some technical conditions:
Assumption (ref)(i) is the standard unconfoundedness condition. Although we are interested in CATE for a low-dimensional subset $X_1$ of the covariates, we still use the full vector of $X$ to address selection into treatment. Allowing for $X$ to be high-dimensional makes it more plausible to have conditional independence between the potential outcomes and the treatment indicator. Assumption (ref)(ii) is a usual sufficient condition for the estimation of standard errors. Assumption (ref)(iii) is the overlapping support condition commonly assumed in the literature. We also need it to establish that our CATE$(x_1)$ estimator converges at the usual nonparametric rate.
Let $\mu_0(j,x)=\mathbb{E}[Y|X=x, D=j]$, $j=0,1$. It follows immediately from Assumption (ref) that $E[Y(j)|X_1=x_1]=E[\mu_0(j,X)|X_1=x_1]$, and hence $\tau_0(x_1)$ is identified as $ \tau_0(x_1)=\mathbb{E}[\mu_0(1,X)-\mu_0(0,X)|X_1=x_1]. $
We now state a less obvious but more robust result based on a Neyman-orthogonal moment condition. Given any probability measure satisfying Assumption (ref), let $\tau(\cdot)$, $\mu(1,\cdot)$, $\mu(0,\cdot)$, $\pi(\cdot)$ denote the functions corresponding to $\tau_0(\cdot)$, $\mu_0(1,\cdot)$, $\mu_0(0,\cdot)$, $\pi_0(\cdot)$, respectively. Let $\eta=(\pi(\cdot),\mu(1,\cdot),\mu(0,\cdot))$ represent the infinite dimensional nuisance parameters needed to identify CATE, and define \[ \psi(W;\eta)=\frac{D(Y-\mu(1,X))}{\pi(X)}+\mu(1,X)-\frac{(1-D)(Y-\mu(0,X))}{1-\pi(X)}-\mu(0,X). \] The following theorem gives a moment condition that is (at least approximately) satisfied at $(\tau_0, \eta)$ even when $\eta$ deviates from $\eta_0$.
\paragraph{Remarks:}
The following assumption describes the properties and use of the sample data:
We now propose two versions of the CATE estimator, depending on whether the first-stage approximation to $\eta_0$ and the second-stage local linear regression targeting $\tau_0$ take place over the same sample or not.
We use the local linear smoother studied by f92,fg92,f93 to estimate $\tau_0(x_1)$ in the second stage. While it is possible to extend out results to local polynomial estimators with an extra degree of smoothness, we focus on the linear case for simplicity.\footnote{An early version of the paper considered kernel-based (local constant) nonparametric regression in the second stage. The results are available upon request.}
Our second-stage estimator is related to the partial mean estimator studied by n94b and lee18. However, lee18 and our paper are distinct in two important ways. First, the parameters of interest, and thus the estimators, are different. We are interested in CATE$(x_1)$ when the treatment variable is binary, while lee18 considers a model with a continuous treatment. Second, as explained above, our analysis is compatible with the use of high-dimensional data. For the full-sample first-stage estimation, we allow the complexity of our first-stage estimator to increase with the dimensionality of the data and investigate its impact on the rate of convergence. For the split-sample first-stage estimation, we show that the impact of the increasing complexity is eliminated due to the independence between the observations used in the first- and second-stage estimations.
When $X_1$ is discrete and takes the values $x_{1,1},\cdots,x_{1,M}$, the function CATE($x_{1,m}$), $m = 1,\ldots,M$ can be interpreted as the average treatment effect for the subpopulation $X_1 = x_{1,m}$. In this case one can restrict the sample to observations with $X_1 = x_{1,m}$, and directly apply the full-sample or cross-fitting estimation methods developed in BCFH13 and CC17.
In this section we first provide the fundamental asymptotic results for our CATE estimators which form the basis of the uniform inference procedures to be given in Section (ref). To this end, we state and discuss several assumptions. Let $\mathcal{X}_1\subset \mathbb{R}^d$ denote the support of $X_1$ and let $\overline{\mathcal{X}}_1$ be the subset of $\mathcal{X}_1$ over which $\tau_0(x_1)$ is to be estimated. In addition, let $f(x_1)$ denote the p.d.f.\ of $X_1$.
For the most part, Assumption (ref) is a collection of standard regularity conditions used in the nonparametric treatment effect estimation literature. The functions $f(x_1)$, $\mu_0(0,x_1)$, and $\mu_0(1,x_1)$ are required to be sufficiently smooth over $\overline{\mathcal{X}}_1$, the density of $X_1$ must be bounded away from zero over the same set, and the kernel function $\mathcal{K}$ must obey some mild restrictions, satisfied by usual choices of $\kappa$ such as the Gaussian or the Epanechnikov kernel (in the simulations and the empirical study we use the former). Of course, Assumption (ref)(ii) also implies that we restrict attention to the technically more interesting case in which the distribution of $X_1$ is continuous, which means that one cannot simply use sample splitting to estimate CATE at various points in the support of $X_1$.
The conditions imposed on the bandwidth in Assumption (ref)(iv) are motivated as follows. The restriction $H>1/(4+d)$ means that $h$ converges to zero faster than the MSE-optimal bandwidth choice; this undersmoothing condition is needed to ensure that the bias from the second-stage kernel regression is asymptotically negligible. In addition, we require $H<(1-2/q)/d$ to be able to use a Gaussian approximation as in CCK14. If the outcome variable is bounded, one can set $q=\infty$ in Assumption (ref)(ii) so that $H<1/d$ as in CCK14. If one only assumes $q=4$, then the convergence rate must satisfy $H \in (1/(4+d),1/2d)$. For this interval to be nonempty, $d$ can be at most 3, which is consistent with Assumption 1 in LOW. In principal, it would also be possible to use the optimal bandwidth, i.e., $H = 1/(4+d)$, and conduct bias correction as in c18, while accounting for the impact of the estimated bias correction term on the standard error. This approach is, however, beyond the scope of the present paper. A key conceptual difference between our setup and c18 is that in our case the dependent variable is not directly observed, but is rather constructed based on first-stage nuisance estimators. Finally, Assumption (ref)(v) is a standard bound for the Taylor remainder, commonly assumed in local linear regression theory. See, for example, lr07.
We now state high-level conditions that specify the convergence rates required of the first-stage nuisance function estimators. The stated rates are linked to the bandwidth sequence $h$ used in the second-stage regressions. More specifically, we make the following assumption about the full-sample first-stage estimator $\hat\eta(I)$.
\paragraph{Remarks:}
The corresponding assumption about the cross-fitting (split-sample) estimator is as follows.
\paragraph{Remarks:}
Theorem (ref) provides the linear (Bahadur) representations of the nonparametric estimators $\hat{\tau}(x_1)$ and $\check \tau(x_1)$ with uniform control of the remainder terms. It serves as a building block for both pointwise and uniform inference about $\tau_0(x_1)$.\footnote{In the online supplement, we provide the linear (Bahadur) representations of the nonparametric estimators $\hat{\beta}(x_1)$ and $\check \beta(x_1)$ with uniform control of the remainder terms, which can be of independent interest.} Starting with the former, we define \[ \sigma^2_N(x_1)=h^d Var\left(\frac{1}{h^df(x_1)}(\psi(W,\eta_0) - \tau_0(x_1))\mathcal{K}_h\left(X_1-x_1\right)\right), \] and suppose that $\sigma^2_N(x_1)$ satisfies:
Then Theorem (ref), together with Lyapunov's CLT, implies
for any fixed $x_1\in\overline{\mathcal{X}}_1$. One can estimate the variance $\sigma^2_N(x_1)$ as \[ \hat\sigma^2_N(x_1)=\frac{1}{Nh^d\hat f^2(x_1; I)}\sum_{i=1}^n\big(\psi(W_i,\hat\eta(I)) - \hat\tau(x_1)\big)^2\mathcal{K}^2_h\left(X_{1i}-x_1\right), \] and we will show that
Of course, this means that inference in practice can proceed based on ((ref)) with $\hat\sigma_N(x_1)$ replacing $\sigma_N(x_1)$. Furthermore, result ((ref)) remains valid if one uses the estimator $\check \tau(x_1)$ in place of $\hat\tau(x_1)$; in this case $\sigma^2_N(x_1)$ can be estimated as
As can be seen from the proof, the $o_p(1)$ term actually vanishes polynomially in $N$.
While the high-level assumptions stated in Section (ref) can accommodate multiple machine learning procedures for estimating $\eta_0$, here we describe the first stage using lasso estimation as a leading example. We now discuss some primitive conditions under which lasso estimation of $\eta_0$ will satisfy Assumptions (ref) and (ref). Specifically, let $b(X)=(b_1(X),\ldots,b_p(X))$ be a dictionary of control terms based on $X$, where $p$ is potentially larger than the sample size $N$ and can grow with $N$.\footnote{To be fully consistent with the general notation, it would be more precise to denote the dictionary as $X=b(X^*) = (b_1(X^*), \ldots, b_p(X^*)) $; see the discussion in the second paragraph of Section 2. We opt for simplicity at a small cost in notational consistency.} Typically, $b(X)$ consists of $X$, and powers and interactions of the components of $X$. The lasso approximates the nuisance functions $\eta_0$ with linear combinations of the components $b_i(X)$; in particular, for $p$-vectors $\beta$, $\alpha$ and $\theta$, set
where $\Lambda(\cdot)$ is the logistic c.d.f. A primitive condition that justifies using the lasso is approximate sparsity. Intuitively, this means that it is possible to make the approximation errors $r_\alpha$, $r_\beta$, $r_\theta$ small with just a small number of approximating terms, i.e., with $\alpha$, $\beta$ and $\theta$ having only a handful of non-zero components.\footnote{The linear index structure and approximate sparsity are specific to the lasso; other machine learning methods provide different types of approximations which do not necessarily rely on sparsity.} The coefficients $\alpha$, $\beta$ and $\theta$ are estimated by penalized least squares or maximum likelihood, where a penalty is imposed for any non-zero component.
In the lasso computations, we set the tuning parameter to be $2c\sqrt N\Phi^{-1}(1-0.1/(\log(N)2p))$ and $c\sqrt N\Phi^{-1}(1-0.1/(\log(N)4p))$ for the conditional mean and propensity score functions estimation, respectively, following BCH14a and BCFH13.\footnote{The constant values are usually chosen as $c=1.1$. Also in practice, one could use cross-validations to choose the tuning parameter here.}
To formalize the idea that the dimension $p$ of $b(X)$ is comparable with or larger than the sample size, we let $p=p_N$ be a function of $N$ and allow $p_N$ to grow to infinity as $N$ increases, possibly (much) faster than $N$. For example, one could set $p_N=O(N^\lambda)$ for any $\lambda>0$, but even $\log(p_N)=O(N^\lambda)$ is allowed if $\lambda$ is not too large. The linear approximation errors to the components of $\eta_0$, defined in display ((ref)), will be controlled by the sparsity index $s=s_N$, a nondecreasing sequence of positive numbers potentially converging to infinity with $N$. Also needing control is the upper bound on the components of $b(X)$; to this end, let $\zeta=\zeta_N=\max_{1 \le j\le p_N}\|b_j(X)\|_{\mathbb{P},\infty}$, and note that $\zeta$ (weakly) increases as $p_N\rightarrow\infty$. The following assumption formalizes the notion of approximate sparsity.
Part (i) of Assumption (ref) states that the number of nonzero coefficients in the $b(X)$-based linear approximations to $\eta_0$ is at most $s$. Part (ii) requires that the approximation errors associated with these linear combinations asymptotically vanish both in $L_2$ and $L_\infty$ norm. This generally requires $s\rightarrow\infty$, but $s$ needs to stay small relative to $N$ in the sense that $s^2\zeta^2\log(p)/N\rightarrow 0$.
Given Assumption (ref) and additional regularity conditions, results by BCFH13 imply that conditions ((ref)) and ((ref)) hold with
where the last inequality holds because $||X_1 - x_1||_2\mathcal{K}_h(X_1-x_1) \lesssim h$. Furthermore, BCFH13 also establish (ref) with $\delta_{4N} = s$ and $A_N = p$ for the following function classes:
where $\ell_N$ is some slowly diverging sequence, e.g., $\ell_N=\log(\log(N))$ and $M>0$. (As $\pi_0(\cdot)$, $\mu_0(1,\cdot)$, and $\mu_0(0,\cdot)$ are uniformly bounded, $\mathcal{G}_{N}^{(0)}$, $\mathcal{G}_{N}^{(1)}$, $\mathcal{G}_{N}^{\pi}$ have bounded envelope functions.)
Given these results, Assumption (ref) with first-stage lasso estimation boils down to the following conditions:
These conditions all hold if $\frac{s_\mu s_\pi\log^2(p \vee N) \log(N)}{Nh^d} = o(1)$ and $\zeta^2 (s_\mu +s_\pi)^2\log(p \vee N)\log(N) = o(N^{1-2/q})$. For example, if $q=4$, $p=O(N^\lambda)$, $\lambda>0$, and $\zeta=O(N^{1/4})$, then $\max(s_\mu,s_\pi)=o(\sqrt{Nh^d})$ is essentially sufficient for Assumption (ref), ignoring logarithmic factors of $N$.
By contrast, Assumption (ref) holds under substantially weaker sparsity conditions. Given the rates in ((ref)), the l.h.s.\ of ((ref)) is at most of order $O\big(\sqrt{s_\pi}\sqrt{s_\mu}\log(p)/N\big)$, as $\mathcal{K}_h$ is a bounded function. Hence, Assumption (ref) essentially reduces to $\sqrt{s_\pi}\sqrt{s_\mu}\log(p)/(Nh^d)=o\big((\log(N)Nh^d)^{-1/2}\big)$. Again, setting $p=O(N^\lambda)$, $\lambda>0$, and ignoring the logged factors of $N$ gives $s_\pi s_\mu=o(Nh^d)$. This condition is of course satisfied if $s_\pi=s_\mu=o(\sqrt{Nh^d})$, but there can be tradeoffs between the two sparsity indexes. For example, if $s_\pi=O(1)$, i.e., the propensity score essentially obeys a finite dimensional model linear in parameters, then $s_\mu=o(Nh^d)$ is possible, i.e., $\mu_0(j,\cdot)$ can be a function that is substantially harder to approximate. Given Remark 5 after Assumption (ref), we can see that our sparsity conditions for Assumption (ref) are essentially equivalent to those in CS19. On the other hand, LOW16 is based on parametric first-stage estimators with the dimension of the regressors fixed. Therefore, they do not need sparsity conditions (though one could regard the parametric assumption as an extreme form of sparsity).
Other types of lasso methods such as the group lasso by F15 and the penalized local least squares and maximum likelihood methods by suz19 can also be used. One can verify the rate restrictions in a manner similar to the above.
Turning to uniform inference, one option is to construct uniform confidence bands analytically similarly to LOW. We provide an alternative method based on the multiplier bootstrap. Our multiplier bootstrap procedure is computationally efficient and takes the nuisance function estimators from the first stage as given and only recomputes the nonparametric regression estimator(s) from the second stage. This step simply involves a random rescaling of the terms in the sums ((ref)) and ((ref)). As lasso estimation is usually time consuming, our procedure is less costly to implement than, say, a standard nonparametric bootstrap requiring new samples from the original data and recomputing the whole estimator.
To describe the procedure formally, we make the following assumption.
Assumption (ref) is standard for multiplier bootstrap inference. For example, a normal random variable with unit mean and standard deviation satisfies this assumption. The bootstrap is implemented as follows:
The following theorem is the bootstrap version of Theorem (ref), and it forms the basis of our inference procedure.
Theorem (ref) justifies the validity of the multiplier bootstrap in implying that $\sqrt{Nh^d}(\hat{\tau}^b(x_1)-\hat{\tau}_0(x_1))$ converges in distribution to the limiting distribution of $\sqrt{Nh^d}(\hat{\tau}(x_1)-{\tau}_0(x_1))$ conditional on the sample path (data) with probability 1. Therefore, if Assumption (ref) also holds, then, conditional on data,
The same statements of course hold true if $\check\tau^b(x_1)$ and $\check\tau(x_1)$ replaces $\hat\tau^b(x_1)$ and $\hat\tau(x_1)$, respectively.\footnote{In the online supplement, we also shows similar results regarding $\hat\beta^b(x_1)$ and $\check\beta^b(x_1)$ that might be of separate interest.} In addition to pointwise inference, the uniform control of the error term $R_N^b(\cdot)$ in Theorem (ref) makes it possible to employ the multiplier bootstrap for uniform inference. For the rest of the paper, we focus on the inference of $\tau_0(x_1)$. The uniform inference of $\beta_0(x_1)$ can be implemented in the same manner. We propose the following algorithm.
\paragraph{Uniform Confidence Band Implementation Procedure}
The following theorem formally states the asymptotic validity of the confidence regions proposed above.
\noindentRemarks:
We discuss bandwidth choice in practice. To obtaining our theoretical results, we require undersmoothing to eliminate bias asymptotically. When $d=1$ as in the simulations, we suggest setting $h_N=\hat h \times N ^{1/5} \times N^{-2/7}$, where $\hat h=1.06\cdot \hat{\sigma}_{x_1} N^{-1/5}$ and $\hat{\sigma}_{x_1}$ is the estimated standard deviation of $X_1$. The formula for $\hat h$ corresponds to the rule-of-thumb bandwidth with a Gaussian kernel suggested by Silverman1986.\footnote{When $d=2$ or $3$, we suggest setting for $j=1,\ldots,d$, $h_{jN}=\hat h_j \times N ^{1/(4+d)} \times N^{-2/(4+3d)}$ and $\hat h_j=1.06\cdot \hat{\sigma}_{x_{1j}} N^{-1/(4+d)}$ and $\hat{\sigma}_{x_{1j}}$ is the estimator of the standard deviation of the $j$-th element of $X_1$.} The bandwidth selection is done with the entire sample even for the cross-fitting method. Also the same selection method is employed in the empirical application.
We investigate the finite sample properties of the proposed high-dimensional CATE estimators and the inference procedure outlined above using Monte Carlo experiments. Because of the space constraint, this material can be found in the online supplement.
In this section, we employ the proposed high-dimensional CATE estimators to analyze the average effect of maternal smoking on birth weight while allowing for virtually unrestricted treatment effect heterogeneity conditional on the mother's age. Birth weight has been associated with health and human capital development throughout life (BDS07, AC11), and maternal smoking is considered to be the most important preventable cause of low birth weight K87. In recent studies, AHL and LOW both explored this causal relationship using the CATE approach, and found different degrees of heterogeneity by age. Using observations from 3,754 white mothers in Pennsylvania, LOW found that the CATE of smoking is decreasing from 17 to around 29 years of age, but they differ from AHL in that the contrast between young and 30-year-old mothers is still not large.\footnote{As the smoking effect is negative, “decreasing” means that the detrimental effects of smoking become stronger with age.}
Our study improves on these previous investigations by considering a much larger pool of covariates and explicitly incorporating a variable selection mechanism into the estimation. This initial pool consists of a vector $X$ of raw covariates as well as technical regressors (powers and interactions) to account for the fact that the functional form of $\pi_0$ and $\mu_0$ is unknown. By contrast, AHL assume that a low dimensional parametric model (known up to its coefficients) is correctly specified for $\pi_0$, while LOW assume that either $\pi_0$ or $\mu_0$ obeys such a model. While we still assume that $\pi_0$ and $\mu_0$ are sparse functions, we let a data-driven procedure (lasso) select the most relevant regressors.
We start with the same data set as AHL, composed of vital statistics collected by the North Carolina State Center Health Services, and extract the records of first-time mothers\footnote{The motivation for focusing on first-time mothers is discussed in AHL. In effect, the restricted sample enables more credible identification of the causal effect, as there cannot be uncaptured feedback from the previous birth experience to the current one.} between 1988 and 2002. The variables include whether the mother smokes (the treatment dummy), the baby's birth weight (the main outcome variable, measured in grams), the parents' socio-economic information, such as age, education, income, race, etc., as well as the mothers' medical and health records. The dataset includes 45 raw covariates and 591,547 observations in total. Table (ref) summarizes the most important pre-treatment covariates in the data set.\footnote{We drop some covariates from the analysis for various reasons. For example, the mother's weight gain during pregnancy is arguably not a pre-treatment variable, and the Kessner index of prenatal care is basically a function of the number of prenatal visits and the timing of the first visit.}
In this section we estimate the CATE of maternal smoking on the baby's birth weight with mother's age as the conditioning variable. Following AHL and LOW, we estimate CATE separately for black mothers and white mothers. We only report the estimation results for white mothers in this section; the results for the black mothers can be found in the online supplement. The dependent variable $Y$ is the baby's birth weight measured in grams. The treatment dummy $D$ takes on the value 1 if the mother smokes and 0 otherwise. We start from the set of variables displayed in Table (ref), and construct an even larger dictionary $b(X)$ by adding polynomial terms to account for the unknown form of the nuisance functions in a flexible way. Specifically, we include, up to degree 3, the powers and interaction terms of key dummy variables and continuous and integer covariates. We then end up with 792 covariates in total.
With such a large set of covariates, it is not clear which variables are important in estimating the CATE function. The true set of variables which belong to the estimating equations is assumed to be sparse, as discussed in the previous sections. We hence apply the lasso method in BCFH13 to estimate propensity score ($\pi_0$) and conditional mean function ($\mu_0$). We then compute the robust score function $\psi$ for each observation $i$, and run a local linear regression of $\psi_i$ on mother's age evaluated at numerous grid points in the interval [15, 36] (years of age). We use the cross-fitting variant of the estimator, i.e., the nuisance function estimation and the kernel regression take place in different subsamples, and then these roles are rotated. In the empirical study, we use the same $K$ ($=4$) as in the simulations. Granted that the theoretical property of the proposed $K$-fold cross-fitting estimator is the same as the full-sample estimator in large samples, we recommend using sample-splitting estimator with $K=4$ or 5 following the suggestion of CC18. We refer to the resulting point estimates as HDCATE (HD stands for “high-dimensional”).
AHL include the mother's age, education, month of first prenatal visit (=10 if prenatal care is foregone), number of prenatal visits, and indicators for the baby's gender, the mother's marital status, whether or not the father's age is missing, gestational diabetes, hypertension, amniocentesis, ultrasound exams, previous (terminated) pregnancies, and alcohol use as the confounding factors. The variables selected by our first-step estimation are similar to those used in AHL, with some notable differences.\footnote{Given that we use the cross-fitting method, there are $K=4$ first-stage estimates, and each has its own variable selection. The reported set of variables selected in the first stage is the union of the selected variables in the four split-sample first stages.} In the propensity score function, we also select father's age, and father's education, besides the ones used in AHL, but not gestational diabetes and amniocentesis. In the conditional mean function, we have father's education, and the rest overlap with that of AHL.
The HDCATE estimates are displayed in Figures (ref), (ref) and (ref), along with 90%, 95% and 99% confidence bands, respectively. For a given confidence level, we compute two types of intervals. “HDCATE CB” is the proposed uniform confidence band computed according to the algorithm given in Section 4. “PW CB” is a pointwise confidence band, given for purposes of comparison, where the critical value $\widehat{C}^{\text{2-sided}}_\alpha$ is replaced by the corresponding value from the standard normal distribution (e.g., 1.96 for $\alpha=5\%$), and “LOW CB” is the uniform confidence band by that of LOW. The constant function labeled “ATE” represents the estimated average treatment effect across all ages.
Figures (ref)--(ref) show that maternal smoking has a negative effect on birth weight at all ages (the upper bounds of the confidence bands are negative), and the average effect is likely to become more negative with age. For example, the point estimates show that for teenage mothers of age 18 or younger the negative effect of smoking is, on average, less than 180 grams in absolute value. For mothers around age 24, the same effect is $-220$ grams, and it approaches $-250$ above 35 years of age.\footnote{The non-monotonicities in the point estimate between ages 25 and 35 could be due to undersmoothing and the quickly declining number of first-time mothers toward the top of this age range.} Thus, there is substantial variation in the estimated average treatment effect by age. A potential explanation is that older mothers are likely to have smoked for a longer period, and the detrimental effects of smoking are cumulative (the smoking dummy does not control for duration or intensity of smoking). The figures also shed some light on the gains from using the proposed method compared with the original study of AHL. The CATE function changes the shape and location of the treatment effect estimates, especially for the younger mothers, and shows that the estimated treatment effects are always significantly negative using our method.\footnote{ Granted, the choice of bandwidth is different for the two studies, which affects the shape of the estimated heterogeneous treatment curve to some extent, but it is not the key reason for the different results. If we were to use the same bandwidth choice as in AHL, we would still observe the difference, as we mentioned in the main text.} Another important difference is that in this study we provide a valid uniform confidence band.
Examining the confidence bands qualifies the analysis of the point estimate in important ways. In Figure (ref), the lower bound of the 99% uniform confidence band (dashed line) attains its maximum at around 16 years of age, and the value of this maximum lies just below the minimum of the upper bound attained at around age 24. Thus, it is possible to fit a constant function (at about $-185$ grams) inside the uniform confidence bands. Nevertheless, if one is less conservative and uses the 95% or 90% uniform confidence bands displayed in Figures (ref) and (ref), respectively, then it is no longer possible to do so. Thus, there is fairly compelling (statistically significant) evidence that the smoking effect becomes more negative at least between the ages of 16 and 24. Based on the pointwise confidence band, there is some evidence of further decline in HDCATE at higher ages but it is possible to fit constant functions even within the 90% uniform confidence bands over the interval $[25,35]$. (Again note that these bands become rather wide at higher ages due to the relatively small number of observations.) The LOW confidence band is visibly wider than our confidence band, which is consistent with our simulation results.
We advance the literature on the estimation of the reduced dimensional CATE function by proposing that the nuisance functions necessary for identification be estimated by flexible machine learning methods, followed by a traditional local linear regression. The asymptotic theory we develop builds on previous work by BCFH13 and CC18. Nevertheless, the theory requires non-trivial modifications to accommodate local linear regression in the second stage. Moreover, CATE is a functional parameter, and our results can be used to conduct uniform inference through a bootstrap procedure. In line with CC18, we also advocate using the cross-fitting approach to estimate the nuisance functions and conduct the second-stage regression.
Using the proposed methods, we revisited the problem of estimating the average effect of smoking during pregnancy on birth weight as a function of the mother's age. Our results fall in between AHL and LOW in the sense that we do find age-related heterogeneity (unlike LOW), but it is less marked than in the former study. In particular, there is evidence that the negative effect of smoking becomes somewhat more pronounced with age.