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.
81,474 characters · 19 sections · 100 citation commands
\thispagestyle{empty}
\setcounter{page}{1} \onehalfspacing \allowdisplaybreaks
In observational studies treatment effects estimation is a delicate task. Unlike experimental designs such settings heavily rely on what is called a conditional independence or exogeneity assumption. Credible identification of causal effects requires that the biasing effect of confounding variables is purged out by controlling for them. For example, evaluating active labour market programs typically aims at estimating the causal effect of a certain policy like a training program for unemployed on employment outcomes. However, in situations that were not explicitly designed for scientific evaluation of a certain policy the major challenge arises from the fact that policy assignment can in general not be regarded as random. In our example unemployed that are more educated might self select into training programs if participation is voluntary and not controlling for potential confounders can heavily bias the effect one wants to measure. Thus, the selection of the `right' covariates is at the heart of observational studies econometrics.\\ Recent advances in the causal inference literature cope with problems where the dimension of the covariate space is large -- potentially much larger than the number of observations (Belloni_Chernozhukov_Hansen_2014, Farrell_2015, Kennedy_2016, Athey_Imbens_Wager_2018, Chernozhukov_Chetverikov_Demirer_Duflo_Hansen_Newey_2017). This issue may arise in two situations that should be distinguished. First of all, administrative and other datasets make a huge amount of potential confounders available to researchers. Second, including covariates not only in levels but also as higher order polynomials and interactions might help to enhance credible causal identification in observational studies. While the first situation might be seen as a symptom of the age of `Big Data', the second situation basically represents the well-known case of global nonparametric approximation. Being capable to cope with the latter situation without being exposed to the `curse of dimensionality' would be desirable since it would ease functional form dependence in causal inference. In both cases this brings up the challenge that in settings where the dimension of the covariate space is large compared to the number of observations, the researcher has to select the covariates that should enter the model. Typically, applied researchers select their empirical model based on some theoretical prediction, data availability, or in an ad-hoc fashion just having some feeling that a variable should be included. Estimation methods like Lasso labelled as supervised `machine learning' are able to reduce the dimension of the model and find the variables that predict another variable best. Clearly, introducing a machine learning step brings up new issues when thinking about inference. At a first glance, it might look like a rather philosophical question if the uncertainty of an estimator is influenced by the decision of the researcher which model to choose. Actually the potential additional uncertainty coming from model selection can only be correctly estimated if the selection step is performed explicitly -- that is within the controllable statistical framework of machine learning -- and not ad-hoc.
Formally the problem under investigation extends Rubin_1974's (Rubin_1974) causal framework. The often cited `fundamental problem of causal inference' (Holland_1986) is that in the case of binary treatments one observes only one of the two potential realizations. An unemployed either decides to participate in a training program or not. We never observe the counterfactual situation. Let $D_i\in\{0,1\}$ denote the treatment status of observation $i=1,...,n$ and let $Y_i$ denote the observed outcome of interest for that observation. Then the observed outcome can be written as a function of potential outcomes
In an observational study the triple $(Y_i,D_i,X_i)$ is observed where $X_i$ is a vector of potential confounders. The setting of interest involves an untestable conditional independence assumption
Thus, knowing the confounder space $X$ solves the selection into treatment problem discussed above. Statistics of interest in typical applications are given by
Only under the necessary condition that the CIA holds, the difference in the potential outcomes for the whole population (ATE) or the population that has received treatment (ATET) can be identified and the measured statistics may be interpreted as reflecting causal effects. Several estimator classes ranging from the classical parametric outcome regression to inverse probability weighting (IPW) (Hirano_Imbens_Ridder_2003), K-nearest-neighbour matching (Abadie_Imbens_2006), propensity score matching (Heckman_Ichimura_Todd_1998, Dehejia_Wahba_2002) and optimal weighting approaches (notably Graham_Campos_Egel_2012, Hainmueller_2012 and Zubizarreta_2015) have been proposed in the literature. While the traditional outcome regression potentially suffers from functional form misspecification, the other approaches balance the outcome distributions in the treatment and control group conditional on the covariates. In practice the latter approaches enable semiparametric estimation of treatment effects. A standard strategy is to estimate weights parametrically and then identify the treatment effect of interest nonparametrically using the reweighted outcomes. For example Rosenbaum_Rubin_1983 show that conditioning on the so called propensity score defined as $e(x_i)=Pr(D_i=1|X_i=x_i)$ is as good as conditioning on the covariates directly. Angrist_Pischke_2009 among others find that indeed all econometric estimators can be regarded as reweighting observed outcomes with weights being a function of the covariates included. Not relying on a functional form in the nonparametric step is more honest in the sense that it prevents the researcher to draw conclusions from extrapolation in areas where she just does not have any information on. In other words, overlap of the covariate distributions is required which can be expressed as
For good surveys on the literature with low-dimensional covariate spaces see Imbens_Wooldridge_2009 and Athey_Imbens_2017.\\ In the context of high-dimensional data the additional challenge arises that $X\in\mathbb{R}^p$ and $p>>n$. Thus, we allow for a case that was implicitly ruled out for standard parametric or nonparametric estimation techniques. The pioneering work of Belloni_Chernozhukov_Hansen_2014, Farrell_2015 and Belloni_Chernozhukov_FernandezVal_Hansen_2017 uses the semiparametric efficient influence function theory behind the original doubly robust estimation techniques of Robins_Rotnitzky_1995 and Hahn_1998 for the integration of machine learning methods to the framework of treatment effects identification. Since both the outcome and the treatment equation are predicted using machine learning (for an overview of the different algorithms see Hastie_Tibshirani_Friedman_2009), they call their estimation technique double machine learning. De facto they use the classical IPW estimator that is augmented with some outcome regression terms (AIPW). Athey_Imbens_Wager_2018 develop what they call approximate residual balancing (ARB) under a different regime of assumptions. In contrast to AIPW their estimator does not rely on propensity score estimation, but uses suggestions from the optimal weighting literature to weight the debiasing residual from the outcome equation term.
The contribution of this article is twofold. In a first step we outline the statistical concepts necessary to compare different reweighting schemes in the context of high-dimensional covariate spaces. Particularly, we review the new literature on high-dimensional treatment effects estimation from the perspective of traditional semiparametric theory. Understanding the latter as the basis for the machine learning based estimators enables to link them with popular methods like matching. Guided by asymptotical considerations we focus on the efficient influence function approach to treatment effects estimation. This motivates the discussion of recent contributions of Belloni_Chernozhukov_Hansen_2014 and Athey_Imbens_Wager_2018. Further we also present a modified version of the radius matching estimator of Lechner_Miquel_Wunsch_2011 that is claimed to have good statistical properties while being able to address some finite sample concerns of the efficient influence function approach. However, unlike ARB it is suitable for non-linear machine learning estimators. The second step is to validate our basic claims via various simulation experiments. We show that alternative reweighting schemes like the one we propose perform well in finite sample.
The rest of the article is organized as follows. In the two following sections we will motivate a semiparametric estimation structure that is suitable for the incorporation of machine learning estimators. Section (ref) studies a naive approach to combine treatment effects estimation with Lasso prediction. The failure of this approach will motivate the so called doubly robust methods discussed in section (ref). We will theoretically review the different features and assumptions behind AIPW, ARB and our matching based estimator. A second part in section (ref) then uses various Monte Carlo experiments to investigate the properties of the different estimators in finite sample. The last section concludes.
Given previous considerations at a first glance it seems to be a good idea to apply a machine learner directly to a standard treatment effects estimator. To make this point a bit more concrete, assume that a normalized IPW estimator (see Busso_DiNardo_McCrary_2014 for a discussion) is used to estimate
and $n_1$ and $n_0$ are the number of treated and control observations. In the high-dimensional case when $p>>n$ standard methods to estimate the propensity score will be infeasible. However, it can in principle be estimated using any supervised machine learning technique. Lasso always optimizes a $L_1$-regularized version of the standard loss function. For the case of logistic regression the problem can be written as
When the tuning parameter $\lambda$ is zero, the problem transforms to the standard logistic regression optimization. In settings where $p\rightarrow n$ the in-sample fit of the unpenalized model can become arbitrary large by just increasing the number of covariates. Thus, for any $\lambda>0$ the goal is not to optimize the in-sample fit but to minimize the out-of-sample mean squared error (MSE). Thus, Lasso instead of fitting the model optimally in sample aims at achieving good out-of-sample predictions by avoiding overfitting. This means that including additional covariates in the model is penalized in the minimization problem. In fact, the $L_1$-norm will shrink some of the coefficients to zero. It follows that with Lasso one achieves the dimensionality reduction necessary to fit a propensity score model in the case when $p>>n$. When a high $p$ is achieved by series expansion Lasso approximates global nonparametric methods. In some sense machine learning methods overcome the curse of dimensionality introduced by nonparametric methods. However, it is important to notice that these methods should not be confused with nonparametric approaches. Rather by shrinking some parameters to zero Lasso might be seen as a data-driven compromise between fully flexible nonparametric and rich parametric methods. In principle, it is therefore very appealing to fit the propensity score model in a regularized form since this approach represents a trade off between functional form dependence and the dimensionality of the problem.\\ It follows from this discussion that the choice of the tuning parameter is crucial for the performance of the Lasso. Theoretical results are typically derived under a sparsity assumption which means that only few covariates $s$ really matter in the model. Additionally there are usually some assumptions that guarantee that the correlation between feutures is not too strong. For example Bickel_Ritov_Tsybakov_2009 and Belloni_Chernozhukov_2011 show that for $p>n$ under such a set of assumptions the empirical $L_2$ prediction norm\footnote{In the linear model implied by the Lasso the prediction norm upper bounds the prediction error. This result follows from the triangle inequality and suggests the same convergence rates for the prediction error. For a discussion see Belloni_Chernozhukov_2011.} is bounded by
Despite the asymptotic appeal of such well-defined convergence rates, practitioners typically use $k$-fold cross-validation as it mimics the idea of optimizing out-of-sample MSE. In a recent paper Chetverikov_Liao_Chernozhukov_2017 show that under Gaussian errors the prediction norm is bounded by
The authors also show that by relaxing the normality assumption lower convergence rates may be obtained. It follows that asymptotically the price to pay for not relying on theory based tuning parameters are worse convergence rates for the Lasso, though we acknowledge that in finite samples the more data-driven cross-validation based approaches might be less sensitive against violations of the main Lasso assumptions. Alternatively, practitioners also choose the tuning parameter by taking a value that is one standard error to the right of the minimal cross-validation criterion and thus select a smaller than the cross-validation optimal model. Originally the method was proposed for regression trees by Breiman_Friedman_Olshen_Stone_1984 without giving any theoretical justification. Buehlmann_vanderGeer_2011 show that the true model is nested within the chosen model if the tuning parameter is determined via cross-validation. This result immediately implies that when having variable selection and not prediction in mind a larger penalty term might be more appropriate. These issues are discussed in more detail in section (ref) of this analysis.\\ Despite the very appealing nature of the Lasso, we will develop two arguments why a procedure that uses the machine learning prediction of the propensity score solely in a second stage is not a valid estimation strategy.\\
Figure (ref) shows the distribution of the ATE for IPW with a logistic Lasso model for predicting the propensity score. The distribution was generated with a Monte Carlo experiment using a design where both the outcome and the treatment equation are determined by an approximately sparse sequence of covariates for the high-dimensional case $n=2000$ and $p=2000$ (design 1; see section (ref) for details). The distribution is recentered at zero such that an unbiased estimator would have mean zero. Obviously, the estimator is heavily biased and probably non-normal in both specifications. If the relation between the covariates and the outcome equation gets stronger and $R^2_Y$ is increased from 0.2 to 0.8 the bias gets stronger since now covariates left out by the Lasso when estimating the propensity score potentially have a stronger effect on the outcome. Therefore, the need to include them in the model gets more urgent and the failure when not doing so gets heavier. The result complements the results in Belloni_Chernozhukov_Hansen_2014 for the case of treatment effects estimation based on the propensity score.\\ To conclude, the naive idea of just applying machine learning tools to econometrics turns out to be a delicate approach when aiming at inference and causal interpretation. Broadly speaking, the reason for this is that for valid inference typically both the outcome and the treatment equation have to be used -- a feature so called doubly robust estimators incorporate. We proceed by reviewing their properties in the next section.
The beauty of the so called doubly robust treatment effects estimators is that they can be directly derived from semiparametric theory. In fact, all other classes of treatment effects estimators can then be regarded as `only' adjusting their weighting and bias-adjustment terms to some perceived finite sample needs.\\ In the following we will continue with a short review of semiparametric theory that will prove to be useful to understand causal treatment effect estimation using machine learning techniques. In particular, considering the very general question how an asymptotically optimal estimator should look requires the notion of influence functions or curves. Let $\psi$ be a regular asymptotically linear estimator for estimating the effect of $D$ on $Y$ in a parametric outcome regression model. Then its influence function $\phi$ for the random realization $Z_i=(Y_i,D_i,X_i)$ is determined by
Knowing the influence function of the estimator will enable the researcher to infer the asymptotic properties. Consistency implies that $E[\phi(Z)]=0$ and the asymptotic variance will then be given by $E[\phi(Z)\phi(Z)^T]$. There are well-known results for parametric estimators. For example maximum likelihood estimation (MLE) implies that the Hessian converges to the product of scores, the Cram\'{e}r-Rao lower bound is achieved (see for example Wooldridge_2010) and the influence function can be estimated using
The problem for semiparametric treatment effects estimators, however, is slightly more sophisticated for two reasons.
Conceptually, semiparametric theory builds on the separation of the two problems (i) and (ii). We proceed by first considering a problem where the nuisance parameter $\eta$ is of finite dimension in section (ref) and investigate the asymptotic effects of using the estimated nuisance parameter as an input in the identification step. Linking this to the analysis in (ref) that translates the problem to infinite-dimensional nuisance settings requires that we formulate it very generally. Unlike for standard outcome regression problems we replace the concepts of Euclidean with Hilbert space vector geometry. Roughly speaking, in contrast to Euclidean spaces a Hilbert space $\mathcal{H}$ also captures the case when vectors that span this space have elements that are random variables itself. Most importantly the smallest distance between a Hilbert space vector $h$ and a linear subspace of the Hilbert space $\mathcal{U}$ is orthogonal to any other element of that subspace $\mathcal{U}$. The unique point in $\mathcal{U}$ that is closest to $h$ is then defined as the projection $\Pi(h|\mathcal{U})$ of $h$ onto the subspace. For a more rigorous treatment of these concepts we refer to Tsiatis_2006 and vanderLaan_Robins_2003.
Suppose the estimation problem generalizes to $\psi=\{\theta^T, \eta^T\}^T$. Specifically let $\theta$ be the statistic of interest like for example the ATE and $\eta$ the first stage nuisance parameter like for example the propensity score. Since we are interested in the asymptotic behaviour of $\theta$, our goal is to derive the moments of its influence function. Like in the standard M-estimation setting, scores might be particularly helpful. Following Tsiatis_2006 define the scores of the log-likelihood with respect to the two parameters as
Now suppose both score functions span spaces which are denoted as tangent spaces $\mathcal{T}_\theta\subset \mathcal{H}$ and $\mathcal{T}_\eta\subset \mathcal{H}$ where the latter is denoted the nuisance tangent space. A fundamental result of semiparametric statistics is then that all influence functions reside in the orthogonal complement of the nuisance tangent space, i.e. $\phi\in\mathcal{T}_\eta^\perp$. In fact an influence function of $\theta$ is orthogonal to the score function of $\eta$ and $E[\phi(Z)s_\eta(Z,\psi_0)]=0$ (for details see Tsiatis_2006).\\ Since the variance of the estimator is the second moment of its influence function, the influence function with the lowest second moment is of particular interest. The efficient influence function is given by the projection of any influence function on the whole tangent space $\phi_{eff}=\Pi(\phi|\mathcal{T})$. Thus, the efficient influence function incorporates all the information coming from the first order conditions of the parameters. More precisely, the efficient influence function turns out to be a function of the efficient score
Notice that this rather abstract formulation is well-known in parametric M-estimation. In contrast to the case without a first-stage estimation step of the nuisance parameter, the asymptotically optimal estimator is now defined in terms of the efficient score. The projection captures the effect of the nuisance parameter estimation step on the score of the parameter of interest. The efficient score in this sense represents a quantity where this effect is purged out. Thus, if estimating the nuisance parameter has no effect on the score of the parameter of interest, i.e. $\Pi(s_\theta|\mathcal{T}_\eta)=0$ then the standard MLE result without nuisance parameter applies.
The results ((ref)) and ((ref)) can now be generalized for the case where the nuisance parameter is of infinite dimension. Fortunately, Tsiatis_2006 formulates the problem such that the results depend on projections on subspaces of the Hilbert space. Therefore, introducing nuisance parameters of infinite dimension does not require a fundamentally new formulation of the problem. Instead one can use so called parametric submodels as methodological devices to derive the efficient influence function with infinite nuisance. In particular one can imagine a semiparametric model $\mathcal{P}$ as a set of densities such that
A submodel of $\mathcal{P}$ with finite-dimensional nuisance such that the true data-generating density is contained in the set of these submodels is then defined as the parametric submodel. More formally let the set of possible distributions within a submodel be given by $\{P_\epsilon : \epsilon\in \mathbb{R}\}$. If $\epsilon=0$ the parametric submodel equals the true distribution of the semiparametric model (Kennedy_2016). The parametric submodel is a methodological trick that allows to approximate the semiparametric model arbitrarily close but eases the mathematical concepts required drastically. In particular, an estimator for $\theta$ is a valid estimator for a semiparametric model if it is also a valid estimator for every parametric submodel of that semiparametric model. If every $\hat{\theta}$ in a semiparametric model is contained in the set of semiparametric submodels, this also has to apply for every influence function of every estimator for $\theta$. It follows that every influence function of $\theta$ is contained in the set of influence functions of the parametric submodels. This reasoning also carries over to efficiency considerations. If every influence function of the semiparametric model is contained in the influence functions of the parametric submodels, the minimal variance of a semiparametric model that can be achieved is the supremum over all parametric efficient influence function variances.\\ Efficient influence functions can again be derived by constructing the nuisance tangent space. For that purpose one can first of all define score functions for the parametric submodels that span tangent spaces. The closure of these tangent spaces is then the tangent space of the semiparametric model. Thus, we obtain $\mathcal{T}_\eta$ again and the efficient influence function can again be derived by using its efficient score representation.
We will now make this last point a bit more concrete by showing how it can be used to check if a proposed influence function is an efficient influence function. A fundamental result that bridges the gap between the rather abstract concepts discussed above and concrete estimation techniques is given by Newey_1990, Bickel_Klaassen_Ritov_Wellner_1998 and Newey_1994. They show that every influence function of a regular and asymptotically linear estimator has to obey
Thus, one needs to know the pathwise derivative of the estimator of interest and the score of the parametric submodel to check if a proposed function is indeed an influence function of the parameter of interest. The efficient influence function for the binary semiparametric ATE estimation problem is given by
Hahn_1998 indeed shows that equation ((ref)) is an influence function for the binary treatment effects estimation problem (see also Kennedy_2016). Since ((ref)) is an element of the tangent space, results from section (ref) suggest that the proposed form is also the efficient influence function.\\ Obviously, the approach presented here is unsatisfactory in the sense that the form of the efficient influence function has to be guessed and can then be checked to obey the properties of an efficient influence function. The theory for average treatment effects estimation is well-established. However, in other contexts developing estimators directly from asymptotic theory might be especially appealing. For the case of treatment effects estimation Tsiatis_2006 uses the fact that influence functions for a known nuisance are much easier to derive. He develops a theoretical framework how these full-data results can be transformed to yield results when the nuisance has to be estimated.\\ Note that due to the linear form of ((ref)), the efficient influence function for the ATE equals its efficient score. Recently, Chernozhukov_Chetverikov_Demirer_Duflo_Hansen_Newey_2017 use the work of Neyman_1979 to obtain orthogonal scores for various estimation problems. In different frameworks (ref) can be derived by adjusting the original score with the projection term. For the case of ATE estimation this is equivalent to adjusting the original score as suggested by Newey_1994.
A natural estimator using the moment condition $E[\phi_{eff}]=0$ implied by ((ref)) is
where we advice to use standardized weights (see Busso_DiNardo_McCrary_2014) as in ((ref)). Originally termed as doubly robust estimation techniques this class of estimator goes back to papers by Robins_Rotnitzky_Zhao_1994 and Robins_Rotnitzky_1995. The first terms are the same as for the IPW estimator while the latter terms are weighted representations of the outcome regressions for the treated and the nontreated sample. Therefore, ((ref)) is also labelled as Augmented IPW (AIPW). Rearranging the estimator as in appendix (ref) reveals that the estimator incorporates terms like
In order to achieve consistency, they have to converge to zero in expectation. Thus, either the propensity score model or the outcome regression model or both have to be correctly specified to achieve consistency. It follows that the estimator is protected against misspecification if at least one model is correct. This feature lead Scharfstein_Rotnitzky_Robins_1999 to call such an estimator doubly robust.\\ Particularly, double robustness is a result that stems from the fact of combining two moment conditions multiplicatively. Actually this is exactly what follows from the definition of the efficient influence function as the projection of any influence function on the whole tangent space. Incorporating all information in the case of treatment effects estimation means to model the selection into treatment as well as the outcome regression step as nuisance parameters.\\ In the context of machine learning estimators the feature of efficient scores relying on these double moment conditions is particularly helpful. In section (ref) the failure of combining IPW estimation with machine learning estimation of the propensity score was described. In contrast to the efficient influence function motivated AIPW that takes into account all nuisances, IPW solely relies on the propensity score model. The corresponding terms to ((ref)) are given by
Now using the Lasso convergence rates as in ((ref)) or ((ref)) indicates that such terms will not exhibit $\sqrt{n}$-convergence. In particular, ((ref)) diverges if it is scaled up due to the fact that $s\log p\rightarrow\infty$. Hence, it is the decreased Lasso convergence rate that leads to the failure of IPW in this setting. On a more profound level this explains the bad behaviour of IPW with Lasso predicted propensity scores as shown in figure (ref). Belloni_Chernozhukov_Hansen_2014 therefore propose to use AIPW where the multiplicative structure of the moment term ((ref)) protects against the divergent behaviour of a single Lasso prediction. In particular, using the Lasso convergence rate ((ref)) from the plug-in tuning parameters and rescaling term ((ref)) results in a term of order $\frac{s\log p}{\sqrt{n}}$. For consistency this term has to vanish which implies a sparsity condition $\frac{s^2\log^2p}{n}\rightarrow 0$. Similarly using the cross-validation convergence rate ((ref)) implies the more restrictive sparsity condition $\frac{s^2\log^2 p\log^\frac{7}{2}(pn)}{n}\rightarrow 0$. Investigating the properties of their estimator more deeply, Chernozhukov_Chetverikov_Demirer_Duflo_Hansen_Newey_2017 propose different ways to incorporate sample splitting in the estimation approach. Instead of using the same sample for prediction and second stage estimation, they suggest to partition the sample in $k$ different subsamples such that the first stage nuisance prediction step is performed on $k-1$ subsamples and the less sophisticated second stage ATE estimation step on the remaining subsample. Iterating over the different left out possibilities and taking the average of all $k$ estimators finally results in reduced sparsity conditions.\\ A potential problem with AIPW, independent of how the nuisance functions are estimated, is its high sensitivity to misspecification of the propensity score model. In a controversial paper Kang_Schafer_2007 show in a simulation study that slight misspecification of both nuisance parameters may lead to serious bias in estimating the potential (or missing) outcome. Most importantly this bias can be bigger for AIPW than for estimators that only use the outcome or the treatment model in certain situations. Their paper highlights a major drawback of the estimator for practical empirical research. Even though double robustness is a very desirable theoretical property, it rests on the restrictive assumption of at least one nuisance parameter to be fully correctly specified. In the context of machine learning estimators dropping below a minimum $o\left(n^{-\frac{1}{4}}\right)$ rate might be counterbalanced by a better convergence rate for the other nuisance parameter if sample splitting is applied (for details see Chernozhukov_Chetverikov_Demirer_Duflo_Hansen_Newey_2017).\\ A second potential problem of AIPW is that by using the inverse prediction of the propensity score directly, very low and very large values of the propensity score lead to extreme weights. Although appropriate from an asymptotic perspective, slight misspecification for such predictions can increase the bias drastically (see also Cao_Tsiatis_Davidian_2009, Tan_2010 and Vermeulen_Vansteelandt_2015).\\ General forms of misspecification might be more likely when using very sophisticated prediction methods like machine learning that tend to be sensitive to violations of their major identifying assumptions. In the next two sections we will therefore review alternative suggestions that are less closely related to semiparametric theory but more concerned with the finite sample performance of the estimator.
Approximate Residual Balancing (ARB) proposed by Athey_Imbens_Wager_2018 follows from the same estimation structure using the efficient score function. However, instead of weighting the observed outcomes and the outcome model nuisance parameters with the inverse of the propensity score, weights are optimized in order to minimize MSE. The optimal weighting framework most notably goes back to ideas of Graham_Campos_Egel_2012, Hainmueller_2012 and Zubizarreta_2015. The intuition behind this direction of research is that a `good' treatment effects estimator should estimate the potential outcomes by reweighting the observed outcomes in a way that the covariate distributions in both treatment samples are identical. Given ((ref)), if both samples are identical in terms of observable characteristics then calculating the mean difference between the reweighted observed outcomes is a consistent estimator because the samples are then as if they would come from an experiment. However, by balancing the two covariate distributions nearly perfectly the weights across observations may become very volatile. This directly maps into a higher variance of the estimator (Zubizarreta_2015). Thus, the main contribution of this literature is that by controlling the weighting process directly within an optimization that maximizes balancing subject to some second moment upper bound one can achieve good performance in terms of MSE.\\ Consequentely, instead of predicting the propensity score using a machine learner, Athey_Imbens_Wager_2018 solely rely on outcome model predictions and then use the structure as in ((ref)) with bias-variance trade-off optimized weights. For ATE their estimator can be written as
Crucially, the estimator assumes linearity of the outcome regression. This is necessary for two reasons. First of all, the conditional expectations in ((ref)) can then be estimated using a linear machine learner like the Lasso that guarantees certain convergence rates for a sparse model. Secondly, in order to control the bias that arises from Lasso predictions, one needs the estimator to achieve approximate balancing. This requires that terms like $\Vert\sum_{i=1}^NX_i-\sum_{i=1}^ND_i\hat{W}_i^{ARB}X_i\Vert_{\infty}$ need to be bounded and converge at rates fast enough to counterbalance the decreased convergence rate of the Lasso. By controlling for the weighting scheme directly, such rates can be achieved while still keeping the second moments of the weights from getting too large.\\ Thus, the price to pay for asymptotically optimal MSE weights is the linearity assumption that kicks in twice. Without linearity neither the Lasso rate is sufficient nor is the weighting scheme optimal. We therefore expect ARB to drastically fail in settings where the outcome process is highly non-linear and tree-based methods or neural nets are more appropriate estimators. Strictly speaking the ARB might therefore not even be interpreted as a semiparametric estimator but rather as a flexible linear model with a bias adjustment motivated by the efficient score structure.
The previous discussion suggests that there is potential to reweight the observed outcomes in a more sophisticated way than just using the inverse of the estimated propensity score. Essentially we build our considerations on the treatment effects estimation literature where instead of labelling the outcome equation nuisance parameters augmentation terms also takes into account the outcome process as a bias adjustment.\\ Most notably Abadie_Imbens_2006 for K-NN matching find that the vanishing bias of their estimator is a function of the outcome model nuisance parameters. Their key finding is that this bias term can converge to zero at a rate potentially slower than $N^\frac{1}{2}$. Based on this observation Abadie_Imbens_2011 develop a bias adjusted K-NN matching estimator where the outcome nuisances are estimated using a global nonparametric method and are then equally weighted using the $K$ observations used for estimating the potential outcome. Similar to AIPW estimation their estimator is doubly robust but incorporates the desirable properties of matching. Generally the idea is that instead of basing the reweighting for a single potential outcome prediction on only one particular propensity score prediction, one should use more predictions close by. The distance is defined in terms of the covariates of the other observations. This feature gives matching the very intuitive interpretation of estimating treatment effects using the comparison of similar observations that differ in their treatment status.\\ Another class of estimators incorporates the same principal ideas but instead of using the covariates directly smooths around an estimated propensity score. A first example of this kind of propensity score or Kernel matching was introduced by Heckman_Ichimura_Todd_1998. Lechner_Miquel_Wunsch_2011 build on this idea and use a triangular Kernel with data-adaptive bandwidth choice. In reference to Dehejia_Wahba_1998 they call their estimator radius matching (RM). Following the findings in Abadie_Imbens_2006 they also implement a bias adjustment term that, however, only uses the estimated propensity scores as predictors for the outcome nuisance parameters. Finally Abadie_Imbens_2016 derive the asymptotic properties of propensity score matching taking into account the first stage estimation of the propensity score. Unlike the estimator of Lechner_Miquel_Wunsch_2011, their estimator uses a fixed number of matches that are equally weighted. From the perspective of nonparametric estimation the authors establish a result for propensity score matching using a uniform Kernel.\\ The big advantage of all of these estimators is that by smoothing around the best match, they are more robust towards misspecification of a single observation propensity score prediction. Especially for propensity scores close to zero IPW is very prone to large biases since minor misspecification of the treatment effects model lead to large bias weights. Also the critique of Kang_Schafer_2007 should only apply to a lesser extent since neither the treatment nor the outcome model have to be completely correctly specified but only the local average around one of the nuisance predictions. Formal results have been shown by Heckman_Ichimura_Todd_1998 for parametric and nonparametric estimation of the propensity score. More generally, advances in the field of nonparametric regression with generated data also suggest that the second step Kernel smoothing improves convergence rates (see Mammen_Rothe_Schienle_2012). However, to the extent of the author's knowledge there is no result so far that describes the smoothing properties of a Kernel when faced with a machine learning generated input.\\ Machine learning estimators of the nuisance parameters can be seen as a further advancement in robustifying treatment effects estimators against misspecification, though, none of the estimators considered here are directly suitable for combining with machine learning approaches. Therefore we propose an Augmented Radius Matching (ARM) estimator that like RM uses a triangular Kernel with data-adaptive bandwidth choice to smooth Lasso predicted propensity scores around the best match to reweight the observed outcomes. The same procedure is used to weight the outcome nuisance parameters also predicted by Lasso. Drawing on the structure of ((ref)) for ATE we now use the weights
To account for outliers the bandwidth $h$ is implicitly chosen such that some multiple of the 90$\%$ quantile and not the maximum of the distribution of the propensity score distance defines a radius determining the maximum distance up to which observations are taken into account. In this sense the weighting strategy represents a sophisticated Kernel matching approach with data-adaptive bandwidth. For some more concrete specifications of the bandwidth choice we refer to the simulation study in Huber_Lechner_Wunsch_2013 and just follow the recommended parameter choices further investigated in Huber_Lechner_Steinmayr_2015.\\ We conjecture that by combining the extended double robustness property from the smoothing behaviour with the extended robustness towards functional form misspecification in the first step nuisance predictions ARM performs well in finite sample. Additionally like K-NN matching with bias adjustment it should also exhibit the double robustness property necessary to cope with the decreased convergence rates of machine learning estimators. Unlike ARB the estimator is not fine-tuned to minimize MSE in finite sample which may be a disadvantage. However, by still relying on a propensity score weighting, it does not have to make specific assumptions about the linearity of the outcome model and therefore in principle would allow to use non-linear machine learning estimators like random forests or neural nets. Seen like this our estimator represents a compromise between the neat structure of AIPW and weighting concerns that are likely to become more prominent when highly complex methods as machine learning techniques are used.
The variable selection property of Lasso enables to construct an additional class of estimators that builds upon the double selection suggestions in Belloni_Chernozhukov_Hansen_2014. Instead of using the union of selected covariates from Lasso regression on outcome and treatment to select the relevant covariates for the outcome model, one may use this set to refit the propensity score model (see also Athey_Imbens_Wager_2018).
The validity of this approach follows the same reasoning as the double selection procedure for the outcome equation. Since for consistency it does not matter if one controls for all relevant confounders in the propensity score or the outcome model, the estimator should be able to cope with the concerns of single selection discussed in section (ref). However, it is unlikely that the estimator achieves semiparametric efficiency because the propensity score is not estimated nonparametrically. In particular, by including the covariates from the outcome model the refitted parametric model does not approximate a global nonparametric method. An additional disadvantage of this estimator is that it does not allow to use any other machine learning estimator than the Lasso since it explicitly depends on the variable selection feature. Though, from an economic perspective this can also be interpreted as a strength of the estimator because it allows the researcher to gain some intuition what channels should be controlled for when trying to identify causal effects.
To assess the finite sample behaviour of the estimators discussed, we use a specific Monte Carlo set-up. In general the data-generating process (DGP) employed is loosely inspired by designs in Busso_DiNardo_McCrary_2014 for low-dimensional and Athey_Imbens_Wager_2018 and Belloni_Chernozhukov_Hansen_2014 for high-dimensional settings. In particular, we model the selection into treatment with an index model $D_i=I(D_i^*>0)$ where $D_i^*$ follows $D_i^*=\beta_d X_i+c_dv$ where $i=1,...,n$ with $n=2000$ and the dimension of the covariate space is $p=2000$. Similarly the outcome is modelled using $Y_i=D_i\theta+D_iX_i\beta_{g1}+(1-D_i)X_i\beta_{g0}+c_yu_i$.\\ We set the error terms as $
\sim N(0,I)$ and determine the feature matrix as jointly normal such that $X\sim N(0,\Sigma)$. For most specifications a Toeplitz structure of the variance-covariance matrix is assumed. In particular every element of $\Sigma$ obeys $\sigma_{jk}=q^{\left|j-k\right|}$ where $j$ and $k$ are row and column indices. Notice that the ATE can be explicitly derived from the model as
since the first moment of $X$ is zero. However, this does not hold for the ATET except for the case when $\beta_{g1}=\beta_{g0}$ which we exclude such that we obtain effect heterogeneity.\\ Finally, $c_d$ and $c_y$ are scaling factors that are chosen in order to achieve a certain $R^2$ for the two processes. While there is a closed form solution for the treatment process, the parameter value for the outcome process has to be simulated.\\ With the Monte Carlo study we want to approximate four fundamental questions concerning the finite sample performance of the estimators.
Most prominently the relative within design performance of the different estimators shall be compared. Table (ref) depicts the different specifications used.\\
Before considering the performance of the different estimators, one should clarify the question under which machine learning estimator the comparison is made. In fact there are as many Lasso estimators as there are tuning parameter choices. Therefore the three different bandwidth choices discussed in section (ref) are compared in the ideal setting of design 1. The treatment and the outcome process are modelled with $\beta_{j,d}=\beta_{j,g0}=\left(\frac{1}{j}\right)^2$ and $\beta_{j,g1}=\left(\frac{1}{s_yj}\right)^2$ where the index $j$ denotes the column of the covariate in the feature matrix and $s_y\neq 1$ is an arbitrarily chosen parameter to make the effects heterogeneous. Under this DGP the prediction problem is approximately sparse and we expect Lasso based methods to work well. Indeed this is what figure (ref) depicts. Under their best bandwidth choice all methods show good performance. While for the double selection (DSRM, DSIPW) approaches the 1-standard-error rule seems to be dominant, the nuisance prediction approaches (AIPW, ARM, ARB) seem to work best with either the theoretically justified or the cross-validation choice. Intuitively, by fitting a smaller model the 1-standard-error rule Lasso trades a bit more bias against hopefully less variance out of sample. Since in principle the union of both covariates sets would not be necessary to achieve consistency, variable selection with the union is prone to including too much variables in the model and a higher penalty is optimal. Thus, by explicitly relying on the model selection instead of the prediction feature, cross-validation minimization cannot be optimal if it is optimal for prediction. For all further specifications it turns out that the pattern prevails with the exception that the theoretically justified tuning parameter leads to more and more arbitrary results as the Lasso assumptions get violated. Since the plug-in tuning parameter in no setting strongly dominates the cross-validation based methods, we proceed by comparing the performance of DSRM and DSIPW using the 1-standard-error rule and AIPW, ARM and ARB using the minimized cross-validation criterion.
Given these a priori considerations, table (ref) shows that methods which take into account both the treatment and the outcome process perform comparatively well. As already shown graphically in the motivation part and in line with theoretical claims, methods based on simply predicting the propensity score (we also use cross-validation minimization) and using this as the only nuisance are heavily biased. A second observation is that alternative weighting schemes in the context of doubly robust estimators perform better than IPW based methods. We suspect the enhanced smoothing properties of the estimators to be more suitable for the lower machine learning convergence rates. Among the estimators based on the efficient score structure ARB exhibits the best performance in terms of RMSE and ARM has the lowest bias.\\ Figure (ref) shows that the standardized draws from the estimators studied are pretty close to 1000 draws from the standard normal distribution. This indicates that despite potential differences in asymptotic efficiency all estimators are nicely behaved and it is likely that they converge towards a normal distribution.
In designs 2-9 we consider two different forms of misspecification. Instead of including unobserved effects in the model, we introduce misspecification as coming from the deteriorating quality of nuisance parameter prediction. For the Lasso this may be achieved by violating the sparsity (designs 2-7) or the restricted eigenvalue assumption (designs 8 and 9). Before examining the latter by analysing the effects of clusters in the feature variance-covariance matrix, the processes are specified more `densely'.\\ For both the treatment and the outcome process a moderately dense process is characterized by coefficients of the form $\beta_j=\left(\frac{1}{j+10}\right)$ and a dense process by $\beta_j=\sqrt{\frac{1}{j}}$. We notice that the dense case represents an upper bound for the degree of misspecification since for such a DGP the asymptotic variance of the estimator becomes unbounded. Moving from the neat design 1 to these misspecified cases for both processes simultaneously shows that the performance of the estimator first decreases slowly before RMSEs explode. For the moderately misspecified design 2 a similar pattern as in design 1 emerges such that IPW-based estimators are again dominated by the alternative double type estimators. Design 2 represents a more realistic situation where not all assumptions necessary to achieve good predictions are specified but estimators are only moderately biased. In such a setting ARM and ARB perform particularly well and dominate AIPW. The enhanced performance of ARM stems from two sources. Most importantly the smoothing property of the estimator allows a higher degree of misspecification because predicted counterfactual outcomes are based on more than one propensity score prediction. Indeed figure (ref) shows that RM not using any adjustment term is probably inconsistent but performs way better than IPW. Furthermore this property also improves the double robustness mechanism of the estimator. Under both DGPs the correlations between the augmentation term and the propensity score based estimator are stronger for ARM. Thus, whenever RM exhibits positive bias the augmentation term corrects for this by adding a negative term and vice versa. Again the smoothed propensity score allows to give higher weights to the outcome nuisances in areas of the support where it is actually needed most.\\ As the results for design 4 and 5 indicate a potential problem with ARM is its sensitivity regarding the choice of the bandwidth. While it beats AIPW in terms of bias, it is worse in terms of RMSE. A higher bandwidth that trades a bit of bias against a decrease in variance might therefore more appropriate in such settings. Moreover, double selection based estimators exhibit an increasing number of draws where they cannot be estimated at all. Since choosing between the different covariates becomes increasingly difficult for the Lasso, the probability that the union of the two sets consists of very many variables becomes bigger such that the refitted propensity score model can only hardly be estimated as in designs 3,6 and 7.\\ Having either of the two processes to be strongly misspecified (designs 6 and 7), forces the estimator to heavily rely on only one of the processes. When the propensity score is very hard to predict, the matching based models become increasingly worse. Reweighting the outcomes and the predicted outcome nuisances with weights that depend on smoothed distances between propensity scores makes these weights increasingly useless. ARB is particularly well-suited for design 6 since it does not depend on the treatment model and shows decent performance. For design 7 AIPW, ARM and ARB perform similarly. This is in line with theoretical considerations because AIPW and ARM mostly depend on the propensity score weighting and ARB should not exhibit any relative losses as long as the outcome model is linear.\\ Besides the violation of the sparsity assumption, nuisance prediction may also suffer from correlation clusters in the feature variance-covariance matrix. This case may be even more relevant in practice since often high dimensions in covariate spaces are generated by either including classes of variables that actually capture the same economic channels or by just generating interactions and polynomials out of an existing dataset. While especially in the latter case one may believe in the sparsity of the model, variables will exhibit clustered correlation structures by construction. Simulations designs 8 and 9 were generated with stochastic correlation matrices using the method by Hardin_Garcia_Golan_2013. Since design 9 puts the cluster structure on neighbouring covariates with similar importance for outcome and treatment, not much additional bias compared to design 1 is expected. However, the standard errors increase because the nuisance prediction correlation between Monte Carlo draws increases.\\ Design 8 puts the cluster structure randomly such that potentially very different covariates in terms of effects become strongly correlated reflecting the dimensionality scaling discussed before. In general a larger bias compared to design 9 is observed since Lasso now becomes partly unable to discriminate covariates with different effects while standard errors increase less drastic. AIPW, ARM and ARB perform similarly but dominate the double selection based estimators. Again ARM and ARB are significantly less biased than AIPW. However, the differences are not as drastic as in the sparsity based misspecification.
We now consider the interesting DGP that assigns treatment completely at random. Hence, we model the case of an experiment in which the researcher is actually unaware of the fact that treatment is assigned randomly and therefore estimates the treatment effect using one of the discussed selection on observables methods. In particular we model $D^*=v$ such that given this knowledge no covariate is necessary to consistently estimate the treatment effect. In such a setting all estimators with the exception of the radius matching based estimators exhibit considerable bias. The reason is that by smoothing over regions with comparable propensity scores, radius matching should weight the different observations nearly equally, while IPW based methods where the weight one observation receives only depends on the propensity score prediction for this particular observation fail to equalize weights over the sample. Therefore for settings where $R^2_D$ is close to zero, smoothing methods have a built-in mechanism that enables them to be near the simple unconditional mean difference estimator which would be the appropriate estimator in such settings. The double selection estimators in general also exhibit comparatively low bias as they include unnecessary variables that however do not hurt when estimating ATE. However, unsurprisingly all selection on observables estimators fail in comparison to the simple difference in means estimator.\\ The discussion highlights a major problem of plugging in machine learning predictions in the efficient score or its radius matching based differences as not only unnecessarily many covariates enter the estimator but they heavily distort the reweighting scheme.
In this paper we analysed and reviewed the doubly robust treatment effects estimator structure for average treatment effects estimation. In the context of machine learning this particular estimation structure turns out to be suitable for incorporating machine learning by combining less then $\sqrt{n}$-convergence rates for both the treatment and the outcome model. However, instead of relying on IPW weights that follow from purely asymptotic arguments we have argued for more sophisticated weighting schemes that are more robust to misspecification -- a problem potentially very prevalent for sophisticated machine learning estimators. Within our simulation design it turned out that a weighting scheme that is based on Kernel matching performs well in finite sample. In contrast to the other AIPW alternatives considered the proposed estimator should also be feasible for nonlinear specifications of the nuisance parameter by still relying on propensity score estimation.\\ For the sake of brevity, many aspects of machine learning approaches for treatment effects estimation were not considered in this study. While the goal of our Monte Carlo designs was to examine some theoretically implied properties, it would be particularly interesting to investigate the performance of the different weighting schemes in less artificial settings. As a further practical concern the weak stability of model selection based estimators highlights the need for variable importance measures for other supervised machine learning estimators. Although from a statistical point of view only the predictive performance for the nuisance function estimators matters, economists are nevertheless interested in developing a feeling for the different driving forces behind the models.\\ Further any rigorous asymptotic analysis of Kernel weighting schemes would require the exact convergence rates of nonparametric methods with machine learning generated inputs. Additional research in this field seems to be very promising. \printbibliography