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.
86,987 characters · 31 sections · 31 citation commands
Machine Learning Estimation of Heterogeneous Causal Effects: Empirical Monte Carlo Evidence
\onehalfspacing
\doublespacing
Economists and many other professionals are interested in causal effects of policies or interventions. This has triggered substantial advances in microeconometrics and statistics in understanding the identification and estimation of different average causal effects in the recent decades \cite<see, e.g.,>[and references therein]{Imbens2009,Athey2017}. However, in most applications it is also interesting to look beyond the average effects in order to understand how the causal effects vary with observable characteristics. For example, finding those individuals who benefit most from active labor market policies, promotion campaigns or medical treatments is important for the efficient allocation of public and private resources.
In recent years, methods for the systematic estimation of heterogeneous causal effects have been developed in different research disciplines. Those methods adapt standard machine learning methods to flexibly estimate heterogeneity along a potentially large number of covariates. The suggested estimators use regression trees Su2009,Athey2016, Random Forests Wager2017,Athey2017a, the least absolute shrinkage and selection operator (Lasso) Qian2011,Tian2014,Chen2017, support vector machines Imai2013, boosting Powers2018, neural nets johansson2016learning,Shalit2016EstimatingAlgorithms,Schwab2018PerfectNetworks or Bayesian machine learning Hill2011,Taddy2016.\footnote{\citeA{Hastie2009} introduce the underlying machine learning algorithms. \citeA{Athey2018TheEconomics} and \citeA{belloni2014high} provide an overview how those methods might be used in the estimation of average causal effects and other parameters of interest.} Recently, the first applied studies using these methods appeared in economics \cite<e.g.,>{Bertrand2017,Davis2017,Knaus2017,Andini2018,Ascarza2018,Strittmatter2018WhatEvaluation}.
In contrast to the rather mature literature about the estimation of average causal effects, the literature on the estimation of effect heterogeneity is still lacking guidance for practitioners about which methods are well suited for their intended applications. Theoretical asymptotic approximations are currently either not available, incomplete, or they are based on non-overlapping assumptions preventing comparisons of estimators. Furthermore, the available information about the finite sample performance is of limited use to practitioners. Most comparisons are based on data generating processes (DGPs) that are very unrealistic for real applications. One exception is \citeA{Wendling2018} who base their simulation study on data from medical records. However, they focus in their study on the special case of binary outcomes and data structures that are unusual in economics.
In this study, we categorize major approaches from different fields. We distinguish between generic approaches that can be combined with a variety of different off-the-shelf machine learning estimators and estimator specific approaches that modify an existing method. The generic approaches are combined with the machine learning estimators Random Forest and Lasso. This leads to eleven different causal machine learning estimators under investigation. As opposed to standard simulation methods that rely on a synthetic DGP, we investigate the finite sample performance of these estimators in an Empirical Monte Carlo (EMCS) approach \cite<e.g.,>{Huber2013,lechner2013sensitivity}. An EMCS informs the DGPs as much as possible by real data and reduces synthetic components in the DGP to a minimum. We consider six different specifications of the heterogeneous causal effects, two different sample sizes, and DGPs with and without selection into treatment.\footnote{We focus on point estimation and leave the investigation of inference methods for further research.}
Our contribution to the aforementioned literature is three-fold: First, we provide a comprehensive comparison of different estimators and DGPs. Second, we consider the finite sample properties of causal machine learning estimators for effect heterogeneity under DGPs that are arguably realistic at least in some fields of economics. Third, this is the first simulation study that considers also different aggregation levels of the heterogeneous effects. In particular, we consider an intermediate aggregation level between the most individualized causal effects and the average population effect. Such intermediate aggregation levels are important as feasible action rules for practitioners.
Our findings suggest that no causal machine learning estimator is superior for all DGPs and aggregation levels. However, four estimators show a relatively good performance in all settings: Random Forests combined with a doubly robust outcome modification \cite<based on>{Chernozhukov2018}, Causal Forest with local centering Athey2017a, Lasso combined with a covariate modification and efficiency augmentation Tian2014, and Lasso with R-learning Nie2017. All those methods use multiple estimation steps to account for the selection into treatment and the outcome process. Several other estimators may be suitable in specific empirical settings but their performance is unstable across different DGPs. Lasso estimators tend to be more unstable than Random Forests, which frequently prevents them from achieving a normal distribution.
In the next section, we introduce the notation and the estimation targets. In Section (ref), we describe and categorize causal machine learning approaches to estimate heterogeneous causal effects. In Section (ref), we explain the implementation of the estimators. In Section (ref), we discuss the EMCS approach. In Section (ref), we provide our simulation results. The final Section concludes and hints at some topics for future research. Appendices (ref)-(ref) contain supplementary statistics and results. Appendix (ref) describes the data that are used for the EMCS. Appendices (ref) and (ref) provide details about the DGPs and the implementation of estimators, respectively. Finally, Appendix (ref) shows and discusses the full simulation results of all DGPs. We provide code that implements the estimators under investigation in the R package \href{https://github.com/MCKnaus/CATEs}{CATEs} on GitHub.
We describe the parameters of interest using Rubin’s Rubin1974 potential outcome framework. The dummy variable $D_i$ indicates a binary treatment, e.g. participation in a training program. Let $Y_i^1$ denote the outcome (e.g., employment) if individual $i$ ($i=1, ..., N$) receives the treatment ($D_i=1$). Correspondingly, $Y_i^0$ denotes the outcome if individual $i$ does not receive the treatment ($D_i=0$). Each individual can either receive the treatment or not. This means that only one of the two potential outcomes ($Y^d_i$) is observable:
Thus, the individual treatment effect (ITE) $\xi_i = Y_i^1 - Y_i^0$ of $D_i$ on $Y_i$ is never observed. However, the identification of expectations of $\xi_i$ may be possible under plausible assumptions. For example, the identification of the average treatment effect (ATE), $\tau = E[\xi_i]$, or the average treatment effect on the treated (ATET), $\theta = E[\xi_i \mid D_i=1]$, are standard in microeconometrics \cite<see, e.g.,>{Imbens2009}.
The focus of this study is on conditional average treatment effects (CATEs). CATEs take the expectations of $\xi_i$ conditional on exogenous pre-treatment covariates.\footnote{Covariates are also called features or predictors in parts of the machine learning literature.} We call the finest conditioning level that uses all available covariates $X_i$ individualized average treatment effect (IATE),
where $\mu^d(x) = E[Y_i^d \mid X_i=x]$ denotes the conditional expectation of the unobserved potential outcomes. IATEs provide an approximation of ITEs for the set of covariates that are at the disposal of the researcher in a specific application. However, researchers may additionally be interested in intermediate aggregation levels that are coarser than IATEs but finer than ATEs. Especially groups based on a smaller set of pre-defined characteristics, $G_i$, may be of interest if the estimated IATEs need to be summarized for the research community, communicated to practitioners, or acted upon.\footnote{For example, if interest is in gender differences, $G_i \in \lbrace female, male \rbrace$.} We call the effects defined on this aggregation level group average treatment effects (GATE),
The identification of any aggregation level of individual treatment effects in observational studies is complicated by non-random treatment assignment. However, identification of the IATE and any coarser aggregation level is still possible if the observable covariates $X_i$ contain all confounders.\footnote{$X_i$ represents the union of confounders and heterogeneity variables for notational convenience. In principle, they may be completely, partly or non-overlapping \cite<see, e.g.,>{Knaus2017}.} These are covariates that jointly affect the treatment probability and the potential outcomes. Although there are alternative ways to identify the various effects, here we focus on the case where all the confounders are contained in the data available to the researcher. This means that we operate throughout the paper under the following assumptions.
Assumption 1: (Conditional independence): $Y_i^1,Y_i^0 \mathrel{\text{\scalebox{1.07}{$\perp\mkern-10mu\perp$}}} D_i \mid X_i=x$, for all $x$ in the support of $X_i$.
Assumption 2: (Common support): $0 < P[D_i = 1 \mid X_i =x] = p(x) < 1$, for all $x$ in the support of $X_i$.
Assumption 3: (Exogeneity of covariates): $X_i^1 = X_i^0$.\footnote{The potential confounders $X_i^d$ are defined equivalently to potential outcomes.}
Assumption 4: (Stable Unit Treatment Value Assumption, SUTVA): $Y_i=Y_i^1 D_i + Y_i^0 (1-D_i)$.
Assumption 1 states that the potential outcomes are independent of the treatment conditional on the confounding covariates. According to Assumption 2, the conditional treatment probability (often called propensity score) is bounded away from zero and one. Assumption 3 requires that the covariates are not affected by the treatment. Assumption 4 excludes spillover effects between treated and non-treated. Under Assumptions 1-4,
and thus IATEs, GATEs and ATE are identified from observable data. We denote the conditional expectations of the outcomes in one treatment arm by $\mu(d,x) = E[Y_i \mid X_i = x , D_i = d]$, the conditional expectation of the outcome as $\mu(x) = E[Y_i\mid X_i =x]$, and the conditional treatment probability by $p(x) = P[D_i = 1 \mid X_i =x]$.
Equation (ref) shows that the fundamental task is to estimate the difference of two conditional expectations. However, we never observe the differences at the individual level and have to estimate them in two different subpopulations. Thus, the estimation of IATEs is a non-standard machine learning problem. In this section, we present different approaches to target the estimation of IATEs. We distinguish between generic approaches and one estimator specific approach. Generic approaches split the causal estimation problem into several standard prediction problems and may be combined with a large variety of supervised machine learning estimators. On the other hand, Causal Forest \citeA{Athey2017a} is a modification of a specific machine learning estimator to move the target from the estimation of outcomes to the estimation of IATEs.
A straightforward generic approach follows directly from Equation (ref). Conditional mean regressions takes the difference of conditional expectations that are estimated in the two samples of treated and non-treated separately using off-the-shelf machine learning methods to estimate the conditional outcome means $\hat{\mu}(d,x)$:\footnote{This approach is also referred to as T-Learner Kunzel2017,Nie2017 or Q-Learning Qian2011.}
Conditional mean regressions are straightforward to implement. Any supervised machine learning methods for conditional mean estimation may be used. However, their usual target is to minimize the mean squared error (MSE) in two separate prediction problems and they are not tailored to estimate IATEs.\footnote{For an intuition why this is not optimal: Biases that for the same value of $x$ go in the same direction are less harmful than if they go in opposite directions. However, this cannot be accounted for by separate MSEs that are not directly linked (for a new Causal Forest estimator that takes up this theme directly, see \citeA{Lechner2018PenalizedEffects}).} This suggests that they may be outperformed by more specialized methods for this causal problem. Three generic multi-step approaches that target IATE estimation are presented in the following and a framework to summarize them is provided.
\citeA{Abadie2005SemiparametricEstimators} introduces the idea of modifying the outcome to estimate conditional average treatment effects on the treated in studies based on difference-in differences. For IATEs in the experimental and observational setting, the idea is formulated by \citeA{Signorovitch2007} and \citeA{Zhang2012}, respectively. The latter discuss two modifications of the outcome, which we summarize as modified outcome methods (MOM). The first is based on inverse probability weighting (IPW) \cite<e.g.,>{Horvitz1952,Hirano2003} where the modified outcome is
The second is based on the doubly robust (DR) estimator of \citeA{Robins1995SemiparametricData},
The crucial insight here is that $\tau(x) = E[Y_{i,IPW}^* \mid X_i = x] = E[Y_{i,DR}^* \mid X_i = x]$. This means that a regression with one of these modified outcomes and covariates $X_i$ can be used to obtain estimates of IATEs, $\hat{\tau}_{IPW}(x)$ or $\hat{\tau}_{DR}(x)$. In practice, the researcher has no access to the true parameters $p(x)$ and $\mu(d,x)$, the so-called nuisance parameters. The conditional expectations need to be approximated in a first step and plugged into Equations (ref) and (ref). Any suitable prediction method can be used to estimate the nuisance parameters as well as the IATEs.
The asymptotic properties of $E[Y_{i,DR}^*]$ as estimator for ATEs are well understood Belloni2014InferenceControls,Farrell2015,Belloni2017,Chernozhukov2017,Chernozhukov2018. Furthermore, \citeA{Abrevaya2015EstimatingEffects} and \citeA{Lee2017} analyze the properties of estimating $\tau(z) = E[Y_{i,IPW}^* \mid Z_i = z]$ and $\tau(z) = E[Y_{i,DR}^* \mid Z_i = z]$ for a low-dimensional subset of covariates ($Z_i$), respectively. However, both do not consider machine learning to estimate the nuisance parameters. We are currently not aware of theoretical results for the case where nuisance parameters and IATEs are estimated with machine learning.
Simulation evidence of \citeA{Powers2018} suggests that estimators based on $Y_{i,IPW}^*$ may exhibit high variance due to potentially extreme values of the propensity score in the denominator. Estimators based on $Y_{i,DR}^*$ might be more stable, because of the double-robustness property, but this is unexplored until now.
\citeA{Tian2014} introduce the modified covariate method (MCM) for experiments and \citeA{Chen2017} extend it to observational studies. They show that we can estimate IATEs by solving the objective function
where $T_i = 2D_i-1 \in \lbrace -1,1 \rbrace$. The name MCM results from the practice to replace the non-parametric function of the IATE with a linear working model, $\tau(x) = x \beta$. This enables us to rewrite the minimization problem (ref) as
where $X_i^{MCM}=T_i/2X_i$ are the modified covariates. The estimated IATEs are then obtained by $\hat{\tau}_{MCM}(x) = x\hat{\beta}_{MCM}$. The nuisance parameter $p(x)$ needs to be estimated in a first step using any suitable method.
In principle, rewriting (ref) as
allows to apply any machine learning estimator that is able to solve weighted minimization problems. However, we are not aware of any study that notices and pursues this possibility.
MCM does not require to specify any model of the so-called main effects $\mu(x)$ or $\mu(d,x)$. However, \citeA{Tian2014} describe that an estimate of $\mu(x)$ might be useful to increase efficiency. The efficiency augmented version replaces the outcome $Y_i$ in Equations (ref) to (ref) by the residuals $Y_i - \mu(X_i)$. Thus, MCM with efficiency augmentation (EA) requires additionally to estimate the nuisance parameter $\mu(x)$ in the first step. \citeA{Tian2014} show that MCM with a linear working model provides the best linear predictor of the potentially non-linear $\tau(x)$. However, we are not aware of any further theoretical analyses of the statistical properties of this approach.
\citeA{Nie2017} propose R-learning that is based on the partially linear model of \citeA{Robinson1988Root-n-consistentRegression}. It is equivalent to MCM with EA for 50:50 randomization but solves otherwise a different minimization problem to estimate IATEs:\footnote{\citeA{Murphy2003} and \citeA{Robins2004} develop a similar method for optimal dynamic treatment regimes called A-Learning. It is adapted for the binary treatment case by \citeA{Chen2017}. We focus on R-Learning because it explicitly applies machine learning methods in all estimation steps.}
Like for MCM, current implementations consider a linear working model for the IATE Nie2017,Zhao2017 and solve
where $X_i^{RL}=(D_i - p(X_i))X_i$ can be considered as an alternative way to modify covariates. The estimated IATEs are then obtained by $\hat{\tau}_{RL}(x) = x\hat{\beta}_{RL}$. Similar to MCM, (ref) can be rewritten as
and solved with any suitable method after estimating the nuisance parameters in a first step \cite<see also>{Schuler2018AEffects}.
\citeA{Nie2017} show that R-learning can perform as good as an oracle estimator that knows the real nuisance parameters in the special case of solving (ref) with penalized kernel regression. This result requires that the estimators of the nuisance parameters need to be fourth root consistent in the semiparametric case.
One goal of this paper is to structure approaches that estimate IATEs coming from different literatures. One way to put the approaches above on more common ground is by noting that they can all be considered as solving a weighted minimization problem with modified outcomes:
Table (ref) summarizes the weights, $w_i$, and outcome modifications, $Y_i^*$, underlying the different approaches. This common representation is helpful to see and understand the differences of the methods. The two MOM methods require no additional weighting because $\tau(x) = E[Y_{i,IPW}^* \mid X_i = x] = E[Y_{i,DR}^* \mid X_i = x]$. For MCM, $\tau(x) = E[2 T_i Y_i \mid X_i = x] = E[2 T_i (Y_i - \mu(X_i)) \mid X_i = x]$ in the special case of 50:50 randomization of the treatment ($p(x)=0.5$). Any deviating assignment mechanism requires reweighting with IPW weights to control for the deviation from 50:50 randomization. The modified outcome of R-learning is equivalent to MCM with efficiency augmentation if $p(x)=0.5$. However, while the MCM modified outcome does not change with other propensity scores and preserves the interpretation as mean comparison under 50:50 randomization, $E[(Y_i - \mu(X_i)) / (D_i - p(X_i)) \mid X_i = x]$ is lacking such an intuitive interpretation.
Another strand of literature modifies machine learning algorithms based on regression trees breiman1984 to estimate IATEs. We focus on Causal Forest that is a special case of the Generalized Random Forest of \citeA{Athey2017a} as the most recent estimator in this line of research.\footnote{Previous works concerned with estimating IATEs by modifying tree-based methods are \citeA{Su2009}, \citeA{Athey2016} and \citeA{Wager2017}.} Causal Forests build on the idea of Random Forests Breiman2001 and boil down to taking the difference of two weighted means in our case of binary treatments,
where the weights $w_i^d(x)$ define an adaptive local neighbourhood around the covariate value of interest, $x$. These weights are obtained from the tailored splitting procedure that we describe in Section (ref). \citeA{Athey2017a} show that this estimator can be consistent and asymptotically normal for a fixed covariate space.
Section (ref) describes generic approaches to split the estimation of IATEs into several prediction problems. The machine learning literature offers a large variety of potential methods such that the investigation of every possible combination of approaches and machine learning methods is not feasible given our restrictions on computational costs. Thus, this study is restricted to two machine learning methods for the implementation of the prediction problems. They are chosen to be representative for more general approaches. First, we consider Random Forests Breiman2001 that serve as a popular representative for methods that attempt local approximations of conditional mean functions. Second, we consider Lasso Tibshirani1996 as a method that attempts global approximation of conditional mean functions.\footnote{By 'local' we mean that each point of the conditional mean function is approximated by the (weighted) average of neighbouring observations. By 'global' we mean the attempt to approximate the conditional mean by a flexible functional fitted to all data simultaneously.} Both methods are increasingly popular in econometrics and are used in methodological contributions as well as in applications.
We consider the combinations of the generic approaches in Section (ref) with Random Forest and Lasso.\footnote{We consider only 'pure' combinations where all estimation steps are conducted with one of the two machine learning methods and neglect the possibilities to estimate, e.g., the nuisance parameters via Random Forests and the IATEs via Lasso. In principle, this is possible but not pursued due to computational constraints. For the same reason, we do not pursue ensemble methods that combine different estimators for the nuisance parameters or for IATE estimation \cite<see, e.g., >{Rolling2018CombiningEffects,Schuler2018AEffects}.} Following \citeA{Chernozhukov2018}, we apply cross-fitting to all approaches that require the estimation of nuisance parameters. This means that the nuisance parameters and the IATEs are estimated in different samples to avoid overfitting.
Table (ref) summarizes all estimators under investigation. We are currently not aware of a Random Forest implementation that supports weighted minimization. Thus, we implement MCM and R-learning only with Lasso. We add further an infeasible benchmark estimator that has access to the true ITEs and uses them as outcome in a standard prediction problem.
Like all machine learning methods, Random Forests and Lasso involve a variety of choices in the implementation. The following sections briefly explain Random Forests and Lasso, present the details of the implementation and explain the use of cross-fitting. The resulting estimators target the estimation of IATEs. Additionally, we consider the methods to estimate GATEs and ATE as a computationally cheap by-product of our analysis. To this end, we average the estimated IATEs, $\hat{\tau}(x)$, to GATEs by $\hat{\tau}(g) = N_g^{-1} \sum_{i=1}^{N} \mathds{1}[G_i = g] \hat{\tau}(X_i)$, with $N_g = \mathds{1}[G_i = g]$, and to the ATE by $\hat{\tau} = N^{-1}\sum_{i=1}^{N} \hat{\tau}(X_i)$.\footnote{Specialized estimators for GATE and ATE might outperform these aggregate estimators but this is beyond the scope of this paper \cite<see, e.g.,>{Lee2017,Chernozhukov2018}.}
The building block of Random Forests for conditional mean estimation are regression trees breiman1984. Regression trees recursively partition the sample along covariates to minimize MSE of the outcome. This leads to the tree structure and the means in the final leaves are used as predictions \cite<for an introduction, see>{Hastie2009}. However, regression trees are unstable and exhibit a high variance. The Random Forest of \citeA{Breiman2001} addresses this issue by combining a large number of decorrelated regression trees. The decorrelation is achieved by growing each tree on a random subsample (generated either by bootstrapping or subsampling) and by randomly choosing the covariates for each split decision.
The standard regression and probability forests split the trees to minimize the MSE of the observed outcomes. Those trees can then be used to form predictions for a realization of $X_i$. These predictions are formed as a weighted average of the observed outcomes where the weights are larger, the more often the observed outcome shares a final leave with the realization of $X_i$. These kind of forests are required for the conditional mean regression and the modified outcome approaches.
The Causal Forest of \citeA{Athey2017a} follows a similar structure. However, instead of splitting the sample according to observed outcomes, Causal Forests split the samples along the gradient of the mean difference with the pseudo outcomes
where $\bar{D}_P$ and $\bar{Y}_P$ are averages of treatment indicator and outcome, $\hat{\beta}_P$ is the mean difference and $Var_P(D_i)$ is the variance of the treatment in the parent node. This splitting rule is tailored to maximize heterogeneity and produces splits that are used to calculate the weights for the weighted mean comparison in Equation (ref).
We also consider the Causal Forest with local centering. This means that $D_i$ and $Y_i$ in Equation (ref) are replaced by $D_i - p(X_i)$ and $Y_i - \mu(X_i)$, respectively. The nuisance parameters are again estimated in a first step and this partialling out should remove the confounding at a high level before building the Causal Forest. \citeA{Athey2017a} show that local centering can improve the performance of Causal Forests substantially in the presence of confounding.
We implement the regression forests for the conditional mean regression and the modified outcomes as well as the Causal Forests using the R package grf Athey2017a. We provide the forests with 105 baseline covariates for the prediction and set the number of variables that are considered at each split to 70. The minimum leaf size is set to one and one forests consists of 1,000 trees. We build honest trees such that the building of the tree and the estimation of the parameters are conducted in separate samples Athey2016. To this end, we split the sample randomly for each tree in three parts: 25% are used to build the tree, 25% to calculate the predictions, and 50% are left out.
The Lasso is a shrinkage estimator and can be considered as an OLS estimator with a penalty on the sum of the absolute coefficients. The standard least squares Lasso solves the following minimization problem
where $w_i$ are weights, $p$ is the number of covariates and $\lambda$ is a tuning parameter to be optimally chosen.\footnote{For the estimation of the propensity score, we use the equivalent logistic regression.} We obtain the standard OLS coefficients if the penalty term is equal to zero and we have at least as many observations as covariates. For a positive penalty term, at least some coefficients are shrunken towards zero to satisfy the constraint. The Lasso serves as a variable selector because some coefficients are set to zero if the penalty term is sufficiently increased. By incrementing the penalty term to a sufficiently large number, eventually all coefficients besides the constant are zero. The idea of this procedure is to shrink those variables with little or no predictive power to zero and to use the remaining shrunken coefficients for prediction. The degree of shrinkage should be chosen to balance the bias-variance trade-off and is the crucial tuning parameter of the Lasso.
We apply the R package glmnet to produce the predictions in the different approaches Friedman2010. We provide the estimator a set of 1,749 potential covariates including second order interactions and up-to fourth order polynomials.\footnote{We exclude binary variables that represent less than 1% of the observations. Furthermore, we keep only one variable of variable combinations that show correlations of larger magnitude than $\pm0.99$ to speed up computation.} The tuning parameter is selected via 10-fold cross-validation.
Some approaches require the estimation of the nuisance parameters in a first step. We follow \citeA{Chernozhukov2018} and apply cross-fitting to remove bias due to overfitting that is induced if nuisance and main parameters are estimated using the same observations. We implement a 50:50 version of their DML1 procedure. This means that we split the sample in two parts of the same size. In the first half, we estimate models for the nuisance parameters. We take these models to predict the nuisance parameters in the second half. These predicted nuisance parameters are then used in the estimation of the IATEs, $\hat{\tau}_1(x)$. We reverse the role of the two halves to obtain $\hat{\tau}_2(x)$. The estimates of the IATE are then calculated as $\hat{\tau}(x) = 1/2 (\hat{\tau}_1(x) + \hat{\tau}_2(x))$.
Table (ref) above lists all estimators that we consider in this study. This list comprises not all alternatives, because we are not able to consider all estimators that have been proposed or would be possible combinations of the generic approaches and existing machine learning methods. We do not consider methods that are tailored for experimental studies \cite<e.g.,>{Imai2013,Grimmer2017,McFowland2018}. Furthermore, restrictions in computation power force us to commit to approaches where we can leverage synergies in the implementation as illustrated in Figure (ref) of Appendix (ref). Thus, we are not able to consider the X-learner of \citeA{Kunzel2017}, the three conditional outcome difference methods proposed by \citeA{Powers2018}, Orthogonal Random Forests Oprescu2018, Penalized Causal Forests Lechner2018PenalizedEffects, methods based on neural nets \cite<e.g.,>{johansson2016learning,Shalit2016EstimatingAlgorithms,Schwab2018PerfectNetworks}, Bayesian approaches like those based on Bayesian additive regression trees (BART) Hill2011,Hahn2017 or Bayesian forests Taddy2016, and potentially other approaches that we are currently not aware of.
Further, the generic approaches discussed in Section (ref) could be implemented using different machine learning algorithms like Boosting, Elastic Nets, Neural Nets, Ridge or any other supervised machine learning algorithm that minimizes the required loss functions \cite<see for an overview>{Hastie2009}.
The simulation study of \citeA{Wendling2018} is close in spirit to our approach, in the sense that their and our DGP relies as much as possible on real data. They compare eight conditional outcome difference estimators for binary outcomes, i.e., they focus on probability models. Their four DGPs are based on the covariates and the observational treatment assignment of four medical datasets. Thus, the IATEs and the resulting binary potential outcomes are the only components that need to be specified. The outcomes are simulated based on predictions of $\mu(0,x)$ and $\mu(1,x)$ from logistic neural networks \cite<for more details, see>{Wendling2018}.\footnote{\citeA{Nie2017} use a similar EMCS for binary outcomes to assess the performance of different implementations of R-learning. However, they do not estimate the IATE but specify it to depend on two covariates.} This is a realistic approach in the medical context. However, it removes two important features from the true outcome generating process. First, the projection of the outcome on observable covariates removes the impact of unobservable variables. Second, the true error structure is lost by imposing a logistic error term. Our EMCS aims to preserve these features of the data at least for the non-treatment outcome. \citeA{Wendling2018} find that conditional mean regressions \cite<implemented with BART, see>{Hill2011} and causal boosting Powers2018 perform consistently well, while causal MARS Powers2018 and Causal Forests Athey2017a are found to be competitive for complex IATE but perform poorly if the variance of the IATE is relatively low. We implement conditional mean regressions and Causal Forests, but omit causal boosting and MARS because of computational restrictions.
Similar to \citeA{Wendling2018} for the medical context, our study aims to approximate a real application in economic policy evaluation as close as possible. The idea of an EMCS is introduced by \citeA{Huber2013} and \citeA{lechner2013sensitivity}. It aims to take as many components of the DGP as possible from real data. We build this EMCS on 96,298 observations of Swiss administrative social security data that is already used in previous evaluation studies behncke2010caseworker,behncke2010unemployed,Huber2017. In particular, the EMCS mimics an evaluation of job search programs as in \citeA{Knaus2017}. The interest is in the heterogeneous effects of such a program on employment over the 33 months after the program start.\footnote{Appendix (ref) provides more details about the outcomes and the rest of the dataset.}
Before we describe and motivate our EMCS approach, we list the general steps to evaluate estimators for IATEs, GATEs and ATEs in Table (ref). We leave out a validation sample (10,000 randomly drawn observations) to compare the estimated IATEs against the true ITEs, while previous EMCS compare in-sample estimates to true values of a known IATE. This modification is intended to focus on the out-of-sample predictive power of the estimated causal effects. The advantage of this procedure is that we can specify the ITEs as ground truth without knowing the IATE as we describe below.
After removing the 10,000 observations of our validation sample, the remaining 78,844 observations form our 'population' from which we draw random subsamples of size 1,000 and 4,000 for estimation. We replicate this 2,000 times for the smaller and 500 times for the larger samples. The precision of the estimators and the computational costs increase with the sample size. Thus, we reduce the number of replications when we increase the sample size to restrict the latter. In case of $\sqrt{N}$-convergence, this will keep the simulation error approximately constant. Table (ref) below shows the variants of the DGP for different $N_s$, $R$, $p(x)$ and $\xi_i$. Before, we explain the specification of the two latter functions.
The 'population' propensity score is estimated in the full sample with 7,454 treated and 88,844 controls. After this estimation step, all treated are removed from the sample. The specification of the propensity score is taken from \citeA{Huber2017} and estimated using a standard logistic regression. We manipulate the constant to create a 50:50 split into treated and non-treated in the simulated samples.\footnote{We remove the 342 observations with a modified propensity score below 5% and above 95%. We deviate at this point from the real dataset and make the problem better behaved than in reality in terms of common support \cite<see discussion in, e.g.,>{lechner2017practical}. We leave the investigation of performance in the presence of unbalanced ratios and insufficient common support for future studies and focus here on a relatively nice setting to start with.} Appendix (ref) provides the details of the specification of the original propensity score and the distribution of the modified propensity score.
We are not able to observe the ITEs or any of its aggregates in a real world dataset. Therefore, we either need to estimate or to specify them. We choose the latter because estimation might favor similar estimators under investigation. Thus, our goal is to create a challenging synthetic ITE that uses components from real data. In observational studies, the estimators must be able to disentangle selection bias and effect heterogeneity. We make it hard for the estimators by using the 'population' propensity score $p^{HLM}(x)$ directly to calculate the ITEs. To this end, the propensity score is normalized and put into a sine function
where $\varepsilon_i$ is random noise. This highly non-linear function of the propensity score is standardized to have mean zero and variance one before it is scaled by the parameter $\alpha$:
where $\bar{\omega}$ is the mean of $\omega(x)$ and $SD(\omega(x))$ is its standard deviation. Finally, we force the ITEs to respect two features of our outcome variable. This means that they are rounded to the next integer and that they must respect that $\hat{Y}_i^1$ falls between zero and 33.\footnote{The histogram of the (observed) $Y_i^0$ is provided in Figure (ref). Given the censored and integer nature of the outcome, we considered also using Poisson Lasso to estimate the outcome regressions. However, the computation time compared to least squared Lasso is substantially longer, while the predictive performance is very similar for our outcomes. Thus, we chose the least squares version for the simulations.} Thus, the final ITEs take the form
where $\lfloor \cdot \rceil$ indicates that we round to the nearest integer.
$\xi(x,y^0)$ is highly non-linear and complicated due to the logistic function, the sine function and the rounding. Additionally, enforcement of the censoring makes it dependent on $Y_i^0$ that is taken directly from the data and thus depends on the covariates in an unknown fashion. Thus, we know the true ITEs but we do not know the functional form of the true IATE.
The true ITEs depend on the observables $X_i$ and additionally on some unobservables through $Y_i^0$.\footnote{These unobservables do not invalidate the CIA in our simulated samples, as $Y_i^0$ and thus the unobservables are not part of the population propensity score. The alternative to ensure a valid CIA in an EMCS is to keep the true treatment allocation structure and to specify the potential outcomes as function of the observables. This is the approach of \citeA{Wendling2018} that is discussed in (ref).} This means that the estimators approximate $\xi(x,y^0)$ using observables and produce estimates $\hat{\tau}(x)$ of $\tau(x)$. The goal of the EMCS is to figure out which estimators approximate the ITE comparatively well in this arguably realistic setting. The relative performance of the estimators translate then directly into the ability to approximate the unknown IATEs because estimators that minimize the MSE of the ITE also minimize the MSE of the IATE \cite<see, e.g.,>{Kunzel2017}. Note that the aggregation of IATEs to GATEs and ATE in step 6 of Table (ref) can be considered as true values because the influence of $Y_i^0$ is averaged out for them asymptotically. This implies that the MSE of GATEs and ATEs would be approximately zero if $\hat{\tau}(x) = \tau(x)$, while the MSE of ITEs might still be positive in this case.
We consider three different values of $\alpha$ in Equation (ref) to vary the size of the ITEs: $\alpha = 0$ (ITE0), $\alpha = 2$ (ITE1) and $\alpha = 8$ (ITE2). Additionally, we create one specification without random noise and one with error term in Equation (ref), $\varepsilon_i = 0$ and $\varepsilon_i \sim 1 - Poisson(1)$, respectively. Table (ref) reports basic descriptive statistics of the resulting potential outcomes, ITEs, and GATEs. ITE0 without random noise creates a benchmark scenario that is most likely to be informative about which estimators are prone to confuse effect heterogeneity with selectivity. ITE1 leads to a scenario with moderate variance of the resulting ITEs. Their standard deviation amounts to about 14% of the non-treatment outcome. ITE2 produces bigger ITEs with a standard deviation of about 6, which is roughly 50% of the standard deviation of the non-treatment outcome. Thus, they should be less difficult to detect.
\doublespacing
The ITEs without and with random noise are created to be similar in their first two moments. However, they differ substantially in their variation that can be explained by observables. The influence of $Y_i^0$ on the ITEs is substantial because between 27% and 44% of the observations are censored for the non-zero ITEs. This explains why ITE1 and ITE2 without random noise are not deterministic either and therefore not perfectly predictable by $X_i$. Still, the out-of-sample $R^2$ of Random Forest and Lasso predictive regressions shown in Table (ref) document that we can explain between about 50% and 70% of the ITEs with our covariates. With random noise, this explainable part decreases to close to zero for ITE0 and 6.3% for ITE1. We consider the latter to be a more realistic scenario because the individual component is expected to be relatively large.\footnote{For example, \citeA{Djebbari2008HeterogeneousPROGRESA} provide evidence that the ITEs in their applications show only little systematic variation.} Thus, we select ITE1 and ITE2 with random noise as our baseline DGPs additionally to the benchmark scenario ITE0 without random noise.
The first column of Table (ref) shows that we specify the mean of the ITEs, the ATE, close to zero.\footnote{Appendix (ref) shows in detail how the ITEs and potential outcomes are distributed, how the ITEs are related to the propensity score and $Y_i^0$, as well as an interpretation of the simulated selection behavior of caseworkers. Note that the lower standard deviations of $Y_i^1$ compared to $Y_i^0$ result from the censoring that moves mass away from the bounds (see Figures (ref) and (ref)).} Appendix (ref) describes how we aggregate the ITEs into 64 groups with sizes between 32 and 420 observations to specify the true GATEs.
In summary, we consider six different scenarios defined by different choices for the scale of the ITEs and the random noise variables. Additionally to the DGP with selection into the treatment, we consider also the case of an experiment with 50:50 random assignment. These twelve different DGPs are considered for the sample sizes of 1,000 and 4,000 observations leading to a total number of 24 different settings. Table (ref) summarizes all parameter settings in which the eleven estimators are compared.
We consider three major performance measures: mean squared error (MSE), absolute bias ($|Bias|$) and standard deviation ($SD$) for the prediction of each observation $v$ in the validation sample:\footnote{The formulas are written for the ITE. The same measures are used for GATE and ATE.}
Most simulation studies are interested in only few parameters such that the performance measure for each parameter can be reported and interpreted. However, in our case we have 10,000 parameters such that we need to summarize the performance over the whole validation sample by taking the averages $\overline{MSE}$, $|\overline{Bias}|$ and $\overline{SD}$.\footnote{For example, $\overline{MSE} = N_v^{-1} \sum_{v=1}^{N_v} MSE_v$.} Additionally, we apply the Jarque-Bera test (JB) to the distribution of predictions for each observation $v$ in the validation sample and report the fraction of observations for which normality is rejected at the 5% confidence level.\footnote{Appendix (ref) discusses and provides also alternative performance measures.}
Table (ref) shows the main performance measures for the three baseline DGPs.\footnote{The full tables with more performance measures are provided in Tables (ref) for ITE0, (ref) for ITE1, and (ref) for ITE2.} At first, we compare estimators within similar approaches to identify the competitive versions. Afterwards, we compare the competitive versions over all approaches to identify those estimators that show an overall good performance and provide a general comparison of Random Forest and Lasso based methods.
\doublespacing
The Random Forest version of conditional mean regressions clearly outperforms the Lasso version in terms of mean MSE. The differences are particularly striking in the smaller sample estimation of ITE0 where the mean MSE of the Lasso version is more than three times larger compared to the Random Forest version. The substantially worse performance of the Lasso version is consistently observed over all baseline DGPs and sample sizes. This is is mostly driven by a substantially lower mean SD of Random Forest based conditional mean regressions that is thus the dominant choice within the two considered versions of conditional mean regressions.
The ranking of the MOM estimators depends on the sample size. Table (ref) shows that Random Forests are superior to the Lasso versions in the smaller samples. Especially the DR modification with Random Forest performs well due to relatively low mean SD. In contrast, the Lasso equivalent is by far the worst estimator in the smaller samples. It has up to twice as large mean SD compared to the next worst estimator and shows consequently a very high mean MSE. One potential reason is that only the DR estimators require the estimation of $\mu(d,x)$ as a nuisance parameter. These predictions are then based on only 250 observations when using cross-fitting, while $\mu(x)$ and $p(x)$ are based on 500 observations. The instability of the Lasso as outcome predictor in small samples seems to spillover to the IATE estimation.\footnote{\citeA{Chernozhukov2018} observe a similar problem of global approximations for the estimation of average effects. See also \citeA{Waernbaum2017} for conditions under which a poor approximation of the outcome leads to worse performance of DR estimators compared to IPW. Similarly, \citeA{Kang2007} demonstrate that double robust methods can perform poorly if both nuisance parameters are misspecified.} The results for the larger sample size indicate that the poor performance is a small sample issue. The DR modification with Lasso outperforms the other versions of MOM in ITE0 and ITE1 and is also close to its Random Forest equivalent for ITE2. The good performance is mainly driven by relatively low mean SDs.
As expected from the results of \citeA{Powers2018}, the IPW modification has relatively high mean SD and is therefore not competitive. This is despite the fact that our DGP does not lead to extreme propensity scores and creates thus a relatively favorable setting for IPW. Therefore, the DR modification seems to be in general the dominant choice as long as the Lasso version is not used in small samples.
The results of the three estimators with modified covariates are similar to the results for the MOM. The MCM is clearly outperformed by its efficiency augmented version and R-learning that both use the outcome regression additionally to the propensity score as a nuisance parameter. For MCM, the efficiency augmentation more than halves the mean SD in all baseline DGPs. Thus, the additional computational effort is fruitful when using MCM. For all DGPs, efficiency augmented MCM and R-learning perform very similar along all dimensions. This finding is in line with the synthetic simulation in Appendix D of \citeA{Chen2017} who find also very similar results for those two.
The Causal Forest is specialized to maximize heterogeneity in experimental settings but it is not build to explicitly account for selection. Thus, it is prone to choose splits that do not sufficiently remove selection bias. However, Causal Forests with local centering address this problem by partialling out the selection effects in a first step. They are specialized to maximize effect heterogeneity and to account for selection bias. Consequently, they uniformly perform better than Causal Forests. This is driven by a relatively low mean absolute bias, but higher mean SD partly offsets this advantage. However, the differences between the two Causal Forest versions are moderate but the version with local centering is the dominant choice if the goal is to minimize mean MSE. However, the improvement comes at the cost of estimating additionally two nuisance parameters before estimating the Causal Forest.\footnote{Together with the conditional mean regression based on Random Forests, the Causal Forest is thus attractive if computation time is a concern. Appendix (ref) shows that both require very similar computation time and are the fastest Random Forest based estimators under consideration.}
The results in Table (ref) show that no estimator is uniformly superior for all sample sizes and DGPs. However, we can categorize the estimators into those that show a relatively good performance in all settings, the volatile ones with outstanding performance only in particular settings, and those that are never competitive.
The first category comprises Random Forest MOM DR, MCM with efficiency augmentation, R-learning and Causal Forest with local centering. These four estimators are in a similar range over all DGPs and sample sizes and belong consistently to the five best estimators. Thus, they seem to be reasonable choices to estimate IATEs. Causal Forest with local centering is the only one of those four that shows never the best performance in terms of mean MSE. This is driven by a larger mean SD that works against the very competitive mean absolute bias. The feature that unifies all four best performing estimators is that they use propensity score and outcome regressions as nuisance parameters in the estimation process.
The MOM DR with Lasso belongs to the second category because it is very competitive for larger samples but the worst choice in smaller samples. Thus, it remains a risky choice for applications because the critical sample size for good performance may depend on the particular dataset.
Finally, conditional mean regressions, MOM IPW with Lasso and Causal Forest should not be considered in settings like ours if minimizing MSE has a high priority. However, if computational constraints are binding, conditional mean regressions with Random Forests and Causal Forests can be attractive options.
A direct comparison of Random Forest and Lasso is possible for conditional mean regressions and MOM. For the smaller sample size, Random Forest clearly outperforms the Lasso based versions. This is driven by the substantially lower mean SD of Random Forest based estimators. The reason is that the global approximations of Lasso are rather instable for small samples. This instability is reduced for larger samples and the Lasso based MOM performs better than the Random Forest equivalents for ITE0 and ITE1.
This dependence on the sample size is not observed for the Lasso specific estimators MCM with efficiency augmentation and R-learning. Both show competitive performance regardless of the sample size. This is particularly surprising given the highly non-linear ITEs. However, all Lasso based methods are far from being normally distributed. For at least 30% of the validation observations, the JB test rejects normality. For many estimators it is even rejected for all observations in the validation sample, while we would expect only a fraction of 5% to be rejected under normality. Columns 9 in the Tables of Appendix (ref) show that this is due to excess kurtosis, which indicates that the Lasso based methods are prone to produce outliers. It is mitigated for the sample size of 4,000 but still the JB test is rejected for a large majority. This reflects the theoretical results of \citeA{Leeb2005ModelFiction,Leeb2008SparseEstimator} that shrinkage estimators like Lasso exhibit non-normal finite sample distributions.
In contrast, all Random Forest based estimators appear to be approximately normally distributed. This is also reflected by a mean skewness close to zero and a kurtosis close to three. Decently performing Random Forest based estimators might be therefore preferable to slightly better performing Lasso based estimators. The former produce less outliers and seem therefore more reliable and robust in empirical applications as well as more amendable to statistical inference.
Table (ref) shows the main performance measures of GATE estimation for the three baseline DGPs.\footnote{The full tables with more performance measures are provided in Tables (ref) for ITE0, (ref) for ITE1, and (ref) for ITE2. The results for all DGPs are provided in Appendix (ref).} We observe similar patterns as for the IATE estimation and the categorization of estimators in Section (ref) remains by and large the same. The four estimators that show a consistently good performance for IATEs are also good choices for the estimation of GATEs.
\doublespacing
For GATE estimation, we observe a new candidate with outstanding performance in a particular setting. MCM performs remarkably well for ITE2 showing the second lowest mean MSE. The mean absolute bias of MCM is already competitive for the estimation of IATEs of ITE2 in Table (ref). However, the mean SD is more than twice as large compared to the best estimators, which prevents a competitive performance in terms of mean MSE. The averaging of these noisy but relatively unbiased estimators seems to produce a competitive estimator for the higher aggregation level. Still, MCM performs poorly for the other DGPs.
The averaging improves also the performance of locally centered Causal Forests relative to its uncentered version. The results for the estimation of IATEs show that the advantage in terms of mean absolute bias is partly offset by a higher variability. The aggregation step reduces this difference such that the lower bias translates also into a substantially lower mean MSE.
The aggregation leads further to a substantial reduction in the excess kurtosis of all Lasso estimators (see Tables of Appendix (ref)). However, the JB test is still rejected for most observations. Note that we observe for all estimators a substantial amount of bias although the influence of $Y_i^0$ and the irreducible noise is averaged out to a large extent. This indicates that the estimators are not able to completely remove the selection bias, which is particularly problematic if we would be interested in statistical inference. The results for ITE0 without noise and with random assignment in Appendix (ref) provide evidence in this direction.
The results for the estimation of ATEs in Appendix (ref) are mostly in line with those for GATEs. Again, MCM is highly competitive and provides the best performing estimators for ITE2. Also the benefits of averaging the locally centered Causal Forest are observed. The bias is halved compared to their uncentered version while both show similar SDs.
The skewness and kurtosis show that the ATE estimators are mostly normally distributed with mean skewness close to zero and mean kurtosis close to three. The obvious exception is MOM DR with Lasso for which also averaging the IATEs does not mitigate the bad performance due to extreme outliers.
Finally, we note that the comparison of the mean MSE for the two sample sizes indicates that GATEs and ATEs estimators show a substantially faster convergence rate compared to the respective IATE estimators. This indicates that the additional averaging of noisily estimated IATEs results in faster convergence and the ATEs may be estimable with close to parametric rates. However, we do not overemphasize this finding as it is only based on two sample sizes.
The discussions in the previous section focus on the results of the three baseline DGPs. This section summarizes the major insights from the alternative DGPs. Their results are provided in Appendices (ref) to (ref) where we also discuss details and peculiarities of the specific DGPs and aggregation levels. In general, the four estimators that are identified as the best performing for the baseline DGPs belong also to the best performing ones in the alternative DGPs.
For the estimation of IATEs, we observe new candidates that are only successful in particular DGPs with selection into treatment. For example, Random Forest MOM IPW performs very good for ITE2 without noise. However, these and other peculiarities discussed in the Appendix hold only for either mean MSE or median MSE, while the four best performing candidates are usually competitive in both measures. Additionally, we assess whether our findings stay robust if we ignore the natural bounds of our outcome variable when creating the DGP. The results in Appendix (ref) show that this is the case in a DGP that allows treated outcomes outside the natural bounds of the original outcome when defining the ITEs.
We also consider all previously discussed DGPs with random treatment assignment. This means that the estimation problem becomes easier because selection bias is no longer a concern. By and large the results are in line with the respective results for the DGPs with selection into treatment, especially the conclusions about the best performing estimators are not changed. As expected, the mean MSEs for the DGP with random assignment are lower for most of the estimators and thus closer to the infeasible benchmark. This is always driven by a lower mean absolute bias while the mean SDs are very similar to the equivalent DGPs with selection. This suggests that the methods are not able to completely remove the selection bias.
Two other differences to the baseline DGPs are noteworthy. First, MOM DR with Lasso shows competitive performance already in small samples. This indicates that the bad performance is related to large errors made in the outcome and the selection equation in small samples, which is in line with the simulation evidence of \citeA{Kang2007} for DR ATE estimators. Second, Causal Forest and its locally centered version show a nearly identical performance. This illustrates that local centering is only beneficial when there is selection into treatment.
The results for GATE and ATE estimators confirm the observation in the baseline DGPs that IATE estimators with low bias but high variance can provide competitive estimators if they are averaged to higher aggregation levels. In particular, averaging MCM IATEs shows in many alternative DGPs a relatively good performance. However, MCM performs worst in some other DGPs in a non-systematic way. Thus, the results show that noise can be averaged out for these higher aggregation levels but there is no guarantee for this. The estimators that are already successful for the IATEs remain the most reliable choices for GATEs and ATEs.
In general, we find that the differences between the estimators become smaller, the more the IATEs are aggregated. Especially, the SDs become more similar by averaging IATEs such that the differences between the estimators are mainly driven by bias.
This is the first comprehensive simulation study in economics that investigates the finite sample performance of a large number of different causal machine learning estimators. We rely on arguably realistic DGPs that have potentially more external validity than the mostly synthetic DGPs considered so far in the limited simulation literature for these estimators. We consider DGPs with and without selection into treatment. Our main goal is to estimate individualized average treatment effects. Additionally, we report the performance of estimators that aggregate individualized average treatment effects to an intermediate and the population level.
We do not find any single causal machine learning estimator that consistently outperforms all other estimators. However, we do find a group of four estimators that show competitive performance for all DGPs. This group includes the Causal Forest with local centering, Random Forest based MOM DR, MCM with efficiency augmentation, and R-learning. These estimators explicitly use both, the outcome and the treatment equations in a multiple step procedure. The estimators that use the Lasso have heavy tails in the smaller samples.
The best performing estimators for the individualized average treatment effects produce also the most reliable estimators for higher aggregation levels. However, in some settings also noisily estimated individualized average treatment effects with low bias produce competitive estimators for higher aggregates because the noise is averaged out.
Despite relying as much as possible on arguably realistic DGPs, the external validity of every simulation study is uncertain. Future research will show if our findings hold in other empirical settings. Furthermore, it may be possible to improve the performance of each method with more tailored implementations. Finally, we focus in this study on the finite sample performance of point estimates and leave the investigation of inference procedures to future research.