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.
237,077 characters · 33 sections · 274 citation commands
Meta-Learners for Estimation of Causal Effects: Finite Sample Cross-Fit Performance
\thispagestyle{empty}
\setcounter{page}{1}
In recent years there has been a growing interest in the estimation of causal effects using machine learning algorithms, particularly in the field of economics Athey2018. The newly emerging synthesis of machine learning methods with causal inference has a large potential for a more \ul{comprehensive} estimation of causal effects Lechner2019. On the one hand, it enables a more flexible estimation of average effects which are of main interest in microeconometrics Imbens2009. On the other hand, it advances the estimation beyond the average effects and allows for a systematic analysis of effect heterogeneity Imbens2017. Both of these aspects contribute to a better description of the causal mechanisms and thus to a possibly more efficient treatment allocation \parencites{Zhao2012}{Kitagawa2018}{Athey2021}{Nie2021}. Hence, applied empirical researchers can greatly benefit from the usage of machine learning methods ranging from evaluation of public policies and business decisions to designing \ul{personalized interventions} \parencites{Andini2018}{Bansak2018}.
Machine learning estimators as such are, however, primarily designed \ul{for} prediction problems and thus cannot be used \ul{directly} for causal inference. Therefore, new approaches for the estimation of causal parameters using machine learning emerged Athey2019a. In particular, the development of the so-called meta-learners have received considerable attention \parencites(see e.g.)(){Kunzel2019}{Kennedy2020}[or][]{Nie2021a}. Meta-learners decompose the causal problem into separate prediction problems that can be solved by standard machine learning algorithms and subsequently combined to estimate the causal parameters of interest Curth2021. Such an approach is advantageous for several reasons. First, the meta-learners do not modify the objective function of the machine learning methods but rather combine their predictions in order to estimate the causal effect. This enables to directly leverage the superior prediction power of machine learning estimators. Second, the meta-learners are generic algorithms refraining from a specific usage of any particular machine learning method. This allows to apply any suitable supervised learning method for the particular prediction problem at hand. Third, the meta-learners are attractive due to the ease of implementation using standard statistical software. This permits researchers to apply the meta-learners without any potential restrictions due to limited availability in software packages and enables tailored implementation for particular types of data. Despite the attractive features of the meta-learners, there is little guidance for applied empirical researchers on how to choose from a variety of the meta-learners proposed in the literature, with lack of unifying simulation evidence for an assessment of the performance of the meta-learners in applied empirical settings.
The complexity of the meta-learners \ul{for the estimation of the causal parameters} varies widely and often hinges on \ul{a prior} estimation of the conditional means of the outcome and the treatment, \ul{also referred to as the nuisance functions, that are not of primary interest} Chernozhukov2018c. Due to the machine learning, \ul{or in general flexible} estimation of such nuisance functions the meta-learners are prone to the overfitting bias, i.e. own observation bias, \ul{which stems from fitting the data too well such that the prediction performance gets compromised} Hastie2009. \ul{This bias then pollutes the estimation of the causal parameters, when the machine learning estimators of the nuisance functions are directly plugged into the estimation of the causal effect using the same data} Newey2018. Therefore, sample-splitting has been proposed in the literature to reduce the overfitting bias by using one part of the sample for estimation of the nuisance functions and the other part of the sample for estimation of the causal effect. \ul{In order to regain the full sample size efficiency of the estimator cross-fitting repeats the estimation by swapping the samples and averaging the estimated causal effects} Chernozhukov2018c. \ul{However, the usage of sample-splitting and cross-fitting is not well understood in practice and the specific definitions of meta-learners differ substantially in their implementation of these procedures.} Despite the ambiguous definitions, there is a lack of simulation evidence concerned with the usage of sample-splitting and cross-fitting within the meta-learning framework and thus limited guidance for or against specific implementations. Moreover, there appears to be limited knowledge about how the asymptotic arguments translate into finite sample properties of the meta-learners.
In this paper, we address both of the above issues and study the finite sample properties of the \ul{machine learning based} meta-learners for the estimation of causal effects based on the specific implementations using the full-sample, sample-splitting and cross-fitting procedures for varying sample sizes. We focus on \ul{evaluating the estimation} of heterogeneous treatment effects \ul{as these provide the most detailed description of the underlying causal mechanisms and thus allow for a better assessment of the individualized impacts of an intervention}. For this purpose, we review the most widely used meta-learning algorithms together with their \ul{assumptions} with respect to sample-splitting and cross-fitting and identify their strengths and weaknesses. We conduct both synthetic and empirically grounded semi-synthetic simulations comparing the performance of the meta-learners in various settings featuring unequal treatment shares, non-linear functional forms and large-dimensional feature sets. \ul{Importantly}, within the simulations we explicitly study the convergence performance of the meta-learners based on growing sample sizes \ul{up to $32'000$ observations.} \ul{Furthermore}, we derive practical recommendations on the choice of specific meta-learners and the respective estimation procedures for applied empirical work.
The results of our simulation experiments reveal that the choice of the estimation procedure has a \ul{large} impact on the performance of the \ul{machine learning based} meta-learners in finite samples. For sufficiently large samples we provide evidence for the theoretical arguments of bias reduction via sample-splitting and cross-fitting, while for smaller samples we observe adverse effects of these procedures \ul{when using machine learning}. The results show that, \ul{if computation time is not a constraint,} cross-fitting is always preferable to sample-splitting as it keeps the bias low, while successfully reducing the variance of the estimators even in small samples. \ul{However}, the results imply heterogeneous impacts of the estimation procedures on the performance of the different meta-learners. \ul{We detect meta-learners for which the performance is stable regardless of the estimation procedure and for which the performance is more sensitive to the choice thereof. This holds not only with regard to bias and variance, but also with regard to convergence rates.} \ul{Beyond the impacts of the estimation procedures, the results reveal clear patterns for the performance of particular meta-learners based on data characteristics. As such, we identify meta-learners suitable for empirical settings with highly unbalanced treatment shares, irrespective of the sample size, as well as meta-learners that are unstable in small samples but have superior performance once larger samples become available. Finally, we set apart meta-learners with undesirable statistical properties that should be avoided for estimation of causal effects.}
This paper contributes to the causal machine learning literature in several ways. First, we provide unifying simulation evidence of meta-learning algorithms for the estimation of heterogeneous causal effects in large-dimensional and highly non-linear settings based on synthetic and semi-synthetic simulations. Second, we explicitly study the meta-learners under the full-sample, sample-splitting and cross-fitting implementations, respectively and thereby provide evidence on the contrast between the asymptotic arguments and finite sample properties. Third, we empirically investigate the convergence \ul{performance} of the meta-learners by repeating the simulation experiments with growing sample sizes. Finally, we derive relevant practical recommendations for applied empirical work which are based on the particular observable data characteristics.
\ul{The rest of} this paper is organized as follows. We briefly discuss the related literature in Section (ref). Section (ref) introduces the notation, the parameters of interest and their identification. Section (ref) reviews the considered meta-learners and the estimation procedures. Section (ref) describes the synthetic as well as semi-synthetic simulations and presents the corresponding results. The main findings of the study are discussed in Section (ref). Section (ref) concludes. Further details including descriptive statistics, an exhaustive summary of the main and supplementary results as well as a computation time analysis are provided in Appendices (ref), (ref) and (ref), respectively.
Within the fast developing causal machine learning literature, many different meta-learning algorithms of different complexities have been proposed. Following the naming convention of meta-learners based on the capital letters, the less complex meta-learners include the S-learner and the T-learner \parencites[see][and]{Lo2002}[][for early ideas of these approaches]{Hansotia2002}[and][for applications in the marketing domain]{Kane2014}\footnote{In the marketing literature the so-called uplift modelling is concerned with the same target of estimating causal effects and developed in parallel to the classical econometric and causal machine learning literature. For a comprehensive overview of the literature on uplift modelling see Devriendt2018. Gutierrez2016 and Zhang2022 provide unified surveys of uplift modelling harmonized with the econometric and causal machine learning literature, respectively.} which besides the treatment effect function do not require estimation of any additional nuisance functions. However, the most prominent and widely used meta-learners in the literature consist of the more complex X-learner Kunzel2019, the DR-learner Kennedy2020, and the R-learner Nie2021a, which all require prior estimation and combination of several nuisance functions, such as the conditional means of the outcome and the treatment, to estimate the causal effect. We define the meta-learners formally and discuss them in more detail in Section (ref). In this paper, we focus on the above-listed meta-learners in order to provide a contrast between less and more complex algorithms for the estimation of causal effects. Moreover, these particular meta-learners have been extensively studied theoretically as well as applied in various empirical settings, including economics \parencites{Knaus2020b}{Jacob2021}{Sallin2021}{Valente2022}, public policy \parencites{Kristjanpoller2021}{Shah2021}, marketing \parencites{Gubela2020}{Gubela2021}, medicine \parencites{Lu2018}{Duan2019} or sports Goller2021b. Some further examples of meta-learners proposed in the literature consist of the U-learner and Y-learner Stadie2018, or more recently the IF-learner Curth2020 and RA-learner Curth2021, which are, however, beyond the scope of this paper.
Besides the meta-learning framework, \ul{there has been also a substantial development of specific causal estimators based on direct modifications of particular machine learning algorithms. Especially, the tree-based estimators have been studied extensively in this respect. These include Causal Trees} Athey2016 \ul{as well as} Causal Boosting Powers2018 \ul{and Causal Forests} Wager2018 \ul{with the extensions of the Modified Causal Forests} Lechner2019 \ul{and the Generalized Random Forests} Athey2019. These methods are based on the underlying predictive algorithms of Regression Trees Breiman1984, Boosted Trees Friedman2001 and Random Forests Breiman2001, respectively. Furthemore, Bayesian versions of Regression Trees\ul{, the so-called BART} Chipman1998 have been adapted for estimation of causal effects as well \parencites{Hill2011}{Taddy2016}{Hahn2020}. Besides the estimators based on recursive partitioning, important causal adjustments have been applied in respect to regularization based estimators such as the Lasso \parencites{Qian2011}{Belloni2013a}{Tian2014} or Lasso-augmented Support Vector Machines Imai2013. Additionally, further machine learning algorithms such as the Nearest Neighbours Fan2018 or Neural Networks \parencites{Johansson2016}{Shalit2017}{Schwab2018}{Shi2019} have been transformed towards causal inference as well. For a comprehensive overview of many of these estimators, we refer the interested reader to Athey2019a, Kreif2019 and Jacob2021, or to Baiardi2021 for numerous empirical examples. While many of these methods have well-established theoretical properties, they restrict the researcher in the choice of the machine learning method. In this paper, although we focus on the machine learning estimation of causal effects, we refrain from an analysis of these methods due to major conceptual differences to the meta-learning framework and the lack of comparability in terms of the usage of sample-splitting and cross-fitting procedures.
In general, the literature on the finite sample properties of causal machine learning estimators under unified framework seems to be rather scarce. One of the exceptions in the econometric literature\footnote{Wendling2018 conduct similar empirically grounded simulation study in medical context, while McConnell2019 use synthetic data.} is Knaus2021 who study a wide range of estimators for heterogeneous as well as (group) average treatment effects, including direct estimators as well as some meta-learners in an Empirical Monte Carlo Study as developed in Huber2013 and Lechner2013. Knaus2021 find no estimator to perform uniformly best, but notice that estimators which model both the outcome as well as the treatment process are substantially more robust throughout all data generating processes considered, which can be observed in our simulations as well. Among the meta-learners considered, the DR-learner and the R-learner perform especially well in terms of the root mean squared error. Moreover, using the Random Forest as a base learner turns out to be more stable with better statistical properties in contrast to using the Lasso, particularly in smaller samples, which also motivates the usage of the Random Forest in our simulations. \ul{However, although} both meta-learners are implemented with cross-fitting, an explicit consideration of different sample-splitting or cross-fitting schemes is missing. A recent work by Naghi2021 also focuses on a variety of causal machine learning estimators, inclusive of selected meta-learning algorithms, in an empirically grounded simulation design. Naghi2021 find Bayesian estimators, including the DR-learner with BART as a base learner, to perform best in estimating the heterogeneous treatment effects. In addition to Knaus2021, the simulations of Naghi2021 provide results on statistical inference in terms of coverage rates and length of the confidence intervals for the estimated causal effects. Nevertheless, albeit Naghi2021, similarly as Knaus2021, implement the relevant estimators using cross-fitting, a devoted assessment of this procedure is omitted. Curth2021 \ul{focus directly on meta-learning algorithms for estimation of heterogeneous treatment effects, but refrain from studying sample-splitting and cross-fitting procedures and rely fully on the full-sample estimation.} In this regard, Zivich2021 study the performance of treatment effect estimators based on cross-fitting, including some meta-learners as well. Similarly to Knaus2021 they find the DR-learner with an ensemble machine learning base learners together with cross-fitting to perform the best among all considered estimators, both in comparison to \ul{cases without} cross-fitting and to parametric base learners. However, Zivich2021 study \ul{exclusively} the estimation of average effects without examining convergence performance of the estimators, considering only a single sample size of $3'000$ observations.\footnote{Numerous other recent studies in the epidemiology literature focus on the estimation of average effects using the DR-learner based on cross-fitting, typically in small samples. For details see Balzer2021, Meng2021, Naimi2021, Zhong2021 and ConzueloRodriguez2022.} Recently, Jacob2020 focuses on the estimation of heterogeneous treatment effects under various cross-fitting schemes for selected meta-learning algorithms. Also, in this simulation study the DR-learner together with the R-learner achieve consistently the best results. Jacob2020 stresses the \ul{heterogeneous impacts} of the particular sample-splitting and cross-fitting procedures on each meta-learner, which is documented in our simulations as well. \ul{Nevertheless, even though considering varying sample sizes within the simulation experiments, the considered sample sizes are limited to $2'000$ observations.} Overall, none of the above studies focuses directly on the convergence performance of the meta-learners under various estimation procedures which still remains an open question. To the best of our knowledge, this is the first paper that \ul{empirically} studies the convergence properties of the meta-learners under full-sample, sample-splitting and cross-fitting implementations with growing sample sizes up to several thousands of observations, \ul{reaching $32'000$ in our simulations.}
In order to describe the effects of interest and their corresponding identification assumptions we rely on the potential outcome framework Rubin1974. We assume a population \(\mathcal{P}\) \ul{from which a realization of} \(N\) i.i.d. random variables is given consisting of a random sample \(\{Y_i(1), Y_i(0),W_i,X_i\} \sim \mathcal{P} \). Here, we consider a binary treatment variable \(W_i\) that is equal to \(1\) for the treated group and equal to \(0\) for the control group, respectively. According to the treatment status we define the potential outcome \(Y_i(1)\) under treatment for the case when \(W_i=1\) and correspondingly the potential outcome \(Y_i(0)\) under control for \(W_i=0\). Additionally, we define a \(p\)-dimensional vector of exogenous pre-treatment covariates such that \(X_i \in \mathbb{R}^p\). Given this definition we can characterize the Individual Treatment Effect (ITE) as follows: \[\xi_i=Y_i(1)-Y_i(0).\] However, the fundamental problem of causal inference is that we never observe both potential outcomes at the same time Holland1986. \ul{Hence, the observed outcomes are defined according to the observational rule as} \(Y_i=Y_i(W_i)\). \ul{The observed data then consists of the triple} \(\{Y_i,W_i,X_i\}_{1\leq i\leq N}\). Nevertheless, it is still possible to identify the expectation of \(\xi_i\) under additional assumptions \parencites(compare)(){Rubin1974}[or][]{Imbens2015}. Thus, we shift the effect of interest towards the Conditional Average Treatment Effect (CATE) which takes the expectation of \(\xi_i\), conditional on covariates \(X_i\) and is given as: \[\tau(x)=\mathbb{E}\big[\xi_i \mid X_i=x\big]=\mathbb{E}\big[Y_i(1)-Y_i(0) \mid X_i=x\big]=\mu_1(x)-\mu_0(x)\] where \(\mu_1(x)=\mathbb{E}[Y_i(1) \mid X_i=x]\) and \(\mu_0(x)=\mathbb{E}[Y_i(0) \mid X_i=x]\) are the response functions \ul{for potential outcomes} under treatment and under control, \ul{respectively}. \ul{In this paper we always refer to the CATE with conditioning on all observed exogeneous covariates and thus focusing on the finest level of heterogeneity} Knaus2021.\footnote{In general, the term CATE describes conditional average treatment effects on various aggregation levels Lechner2019. \ul{In our case, the CATE corresponds to the Individualized Average Treatment Effect (IATE).} Additionally, researchers and especially policy makers might be interested in a \ul{low-dimensional heterogeneity} level based on some pre-specified heterogeneity covariates of interest, which are referred to as the Group Average Treatment Effects (GATEs). Such effects are, however, beyond the scope of our study and the interested reader is referred to Zimmert2019, Jacob2019 and Semenova2021 for a theoretical analysis and to Knaus2021 for simulation based results or to Cockx2019, Knaus2020a, Hodler2020 and Goller2021 for empirical applications estimating policy relevant GATEs.} Kunzel2019 point out that the best estimator for \(\tau(x)\) is also the best estimator for \(\xi_i\) in terms of the mean squared error (MSE).
In order to identify the effects of interest, we need a set of identification assumptions. We operate under the selection-on-observables strategy\footnote{For estimation of heterogeneous treatment effects under different identification strategies see e.g. Athey2019, BargagliStoffi2020 and Biewen2021 for the case of instrumental variables and Zimmert2018, SantAnna2020 and Wager2021 for the case of difference-in-differences.} Imbens2015 and assume that we observe all relevant confounders, i.e. all covariates \(X_i\) that jointly influence both the treatment \(W_i\) and the potential outcomes, \(Y_i(0)\) and \(Y_i(1)\). We state the following identification assumptions:
According to Assumption (ref), also referred to as the conditional ignorability or unconfoudedness assumption, we assume that the potential outcomes are independent of the treatment assignment once conditioned on the covariates, i.e. we assume that there are no hidden confounders. Assumption (ref), also known as the overlap assumption, states that the conditional treatment probability is bounded away from \(0\) and \(1\) and thus it is \ul{possible} to observe treated as well as control units for each realization of \(X_i=x\). Assumption (ref) is known as the stable unit treatment value assumption and indicates that the observed treatment value for a unit is independent of the treatment exposure for other units, which rules out any general equilibrium or spillover effects between treated and controls. Lastly, Assumption (ref) specifies that the covariates are not influenced by the treatment.\footnote{\ul{Analogously to the definition of potential outcomes, we denote potential covariates under control and under treatment as $X_i(0)$ and $X_i(1)$, respectively.}} Under these assumptions it follows that
and thus the CATE can be nonparametrically identified from observable data Hurwicz1950.
\ul{In the machine learning literature meta-learning represents algorithms that exploit knowledge about learning to improve the algorithm's performance, as generally defined by} Vilalta2002. \ul{These include various algorithms that learn to solve new task from prior learning experience, i.e.} learning to learn \parencites{Schmidhuber1987}{Thrun1998}, \ul{algorithms that learn from multiple related tasks, i.e. multi-task learning} Caruana1997, \ul{or algorithms that learn from multiple models solving identical task, i.e.} ensemble learning Dietterich2000.\footnote{\ul{For a recent survey on meta-learning, see} Vanschoren2019.} \ul{Recently, the meta-learning framework has been adopted within the causal machine learning literature for learning causal effects from multiple prediction models} Kunzel2019, \ul{which could be termed accordingly as causal learning.}
At a high level the meta-learners for estimation of heterogeneous causal effects are two-step algorithms. In the first step they define regression functions, in the causal machine learning literature often denoted as the nuisance functions \parencites{Chernozhukov2018c}{Kennedy2020}, which can be estimated by any \ul{supervised learning method fulfilling suitable regularity conditions}, i.e. the base learner.\footnote{The base learners can be in general any set of black-box methods as long as they are consistent estimators of the nuisance functions and sufficiently minimize the prediction error in terms of the MSE \parencites{Kunzel2019}{Kennedy2020}{Nie2021a}. However, in order to provide explicit error bounds and achieve specific convergence rates on the estimation of the causal effect, the base-learners must fulfill a set of regularity conditions such as smoothness or sparsity Alaa2018. We discuss the particular conditions required for each meta-learner in Section (ref).} In the second step they use the estimated nuisance functions to construct an estimator for the causal effect, \ul{i.e. the meta-learner}. Various meta-learners then differ in the definitions of the nuisance functions and their subsequent usage to obtain the final estimator for the causal effects. Depending on the algorithm complexity, some meta-learners require estimation of only one single model whereas others require estimation of multiple models. This raises the question of data usage within the estimation procedure and thus the possible need for sample-splitting and cross-fitting, respectively.\footnote{\ul{Recently, related issue of data usage of the meta-learning algorithms with respect to splitting into training and validation set for the learning to learn domain has been discussed by} Bai2020 and Saunshi2021.}
In general, the nuisance functions are defined as conditional expectations of various types. The most common types are the propensity score function and the response function. First, the propensity score is defined as the conditional \ul{probability} of a binary treatment \(W_i\) given the covariates \(X_i\) as follows: \[e(x)=\mathbb{P}\big[W_i=1 \mid X_i=x\big].\] In the causal inference literature the propensity score plays a central role Rosenbaum1983 in many matching and reweighting methods to balance the distributions of treated and controls \parencites(see)(){Hahn1998}[and][among others]{Huber2013}. Second, the response function is broadly defined as the conditional expectation of an outcome variable \(Y_i\) given a conditioning set of explanatory variables. The particular definitions of the response function then differ in the specification of the conditioning set and the subset of the data used. For the meta-learners studied in this paper, the following definitions of the response function are of interest:
where Equation (ref) defines the response function with conditioning on both the covariates \(X_i\) as well as the treatment indicator \(W_i\), while \(\mu(x,1)\) and \(\mu(x,0)\) describe the response functions with conditioning on the covariates \(X_i\) in the \ul{subpopulation} under treatment \(W_i=1\) and under control \(W_i=0\), accordingly. Similarly, Equation (ref) defines the response function with conditioning only on covariates. The meta-learners then use selected nuisance functions together with the available data as inputs for the estimation of the CATE function which can be generally denoted as follows: \[\tau(x)=\zeta\big(W_i,X_i,Y_i,e(x),\mu(x,w),\mu(x)\big)\] \ul{where $\zeta(\cdot)$ is a function of the respective inputs, which} is detailed for each particular meta-learning algorithm in Section (ref). The problem arises when estimating the nuisance functions using flexible machine learning methods as these are prone to the overfitting bias, i.e. the `own observation bias'. \ul{The overfitting bias emerges when the in-sample data is fitted too well such that the out-of-sample performance is compromised }Hastie2009. \ul{Hence, a single observation $i$ can have a large influence on the predictions for covariates $X_i$ as pointed out by} Athey2019a. Chernozhukov2018c and Newey2018 thus propose sample-splitting procedures that allow for elimination of such overfitting biases.\footnote{Original ideas of using sample-splitting procedures to eliminate own observation bias stem from the literature on density estimation going back to Bickel1982, Bickel1988 and Powell1989 among others.}
Theoretical arguments express the need for sample-splitting when the causal estimator involves several estimation steps such as the estimation of nuisance functions. Within the meta-learning framework the nuisance functions are typically highly complex and potentially high-dimensional functions estimated by supervised machine learning methods such as penalized regression, tree-based methods, neural networks, etc. Using the same data sample for machine learning estimation of the nuisance function as well as for estimation of the causal effect leads to overfitting which induces a bias in the CATE estimator. On a high level, the bias of the CATE estimator can \ul{generally be} decomposed into an estimation error of learning the CATE function itself, and the estimation error in learning the nuisance functions, encompassing the overfitting bias Kennedy2020.\footnote{\ul{The machine learning estimation of the nuisance functions might additionally induce so-called regularization bias. This is due to the effective MSE minimization that restricts the variance of the estimator, but introduces a bias, which might not decay fast enough. For a discussion of the regularization bias in treatment effect estimation, see} Chernozhukov2018c.} Chernozhukov2018c show that \ul{for the ATE estimation} the \ul{overfitting bias} can be controlled by using sample-splitting, \ul{while} Kennedy2020 and Nie2021a \ul{extend this concept for the CATE estimation}. In that case one part of the sample is used to estimate the nuisance functions and the other part is used to estimate the causal effect.\footnote{Sample-splitting procedures are frequently used in causal machine learning literature including Double Machine Learning Chernozhukov2018c, Causal Forests \parencites{Wager2018}{Lechner2019} or the here-discussed meta-learners \parencites{Kennedy2020}{Nie2021a}.} As a result, \ul{the bias term stemming from overfitting} can be shown to be bounded and to converge to zero \ul{sufficiently fast}. Building upon this result, Newey2018 propose \ul{an advanced} sample-splitting scheme called double sample-splitting. In this case, not only the nuisance functions are estimated together on a separate part of the sample but each single nuisance function is estimated on an own separate part of the sample. In practice, the training data is split into \(M+1\) equally sized parts, with \(M\) being the number of nuisance functions to estimate and the remaining part of the data serves for estimation of the causal effect. Newey2018 show that under the double sample-splitting the bias term converges to zero at a faster rate \ul{compared to standard sample-splitting where all nuisances are estimated on the same sample}.\footnote{\ul{The intuition for this result comes from the observation that for estimators using multiple nuisance functions, such as the doubly robust estimators as e.g. the herein discussed DR-learner, the estimation error involves a product of the biases from the estimation of the} \(M\) \ul{nuisance functions. This induces additional nonlinearity bias if all} \(M\) \ul{nuisance functions are estimated using the same data, which gets effectively removed by using separate samples for estimation of each of the} \(M\) functions. For more details see Newey2018 and Kennedy2020.} The double sample-splitting procedure has also been recently implemented by Kennedy2020 in the context of the DR-learner.
In general, the overfitting bias could also be controlled for by restricting the complexity of the nuisance functions which would, however, prevent high-dimensional settings as well as usage of a variety of machine learning estimators or ensembles of those.\footnote{For results in the context of the Lasso estimation under sparsity see Belloni2017.} Hence, the advantage of using sample-splitting is to allow for a \ul{high degree of} complexity of the nuisance functions estimated by a wide class of machine learning estimators Kennedy2020.
It follows that, \ul{theoretically}, sample-splitting prevents overfitting and thus reduces the bias in the final causal estimator \parencites{Chernozhukov2018c}{Wager2018}. At the same time, however, the variance of the estimator increases as less data is effectively used \ul{for estimation}. Cross-fitting Chernozhukov2018c and respectively double \ul{cross-fitting} Newey2018 have been proposed in the literature in order to reduce the variance loss induced by sample-splitting. \ul{In this procedure}, the roles of the data parts get switched such that each part has been used for both the estimation of nuisances as well as the causal effect estimation. The final CATE estimator is then an average of the separate effect estimators produced. This \ul{method} can be further extended to use more than \(M+1\) splits denoted as \(K\)-fold cross-fitting Chernozhukov2018c with the final CATE estimator given as: \[\hat{\tau}(x)=\frac{1}{K} \sum^{K}_{k=1} \hat{\tau}_k(x)\] where \(\hat{\tau}_k(x)\) is the CATE estimator based on the \(k\)-th fold.\footnote{Increasing the efficiency of a sample-splitting based estimator by swapping the roles of the data samples and averaging the resulting estimates goes back to Schick1986 in the context of estimation of semi-parametric models.}
The above theoretical arguments have a direct impact on the implementation of various meta-learning algorithms. \ul{Under the double sample-splitting} the more models have to be estimated within the meta-learning algorithm, the more data splits are being implicitly induced, \ul{while} the impact \ul{thereof} in finite samples is not clear a priori as pointed out by Newey2018. \ul{As such, the} researcher faces a typical bias-variance trade-off with respect to sample-splitting. In order to illustrate the issue it is instructive to \ul{decompose} the MSE of a CATE estimator \(\hat{\tau}(x)\): \[MSE\bigg(\hat{\tau}(x)\bigg) = Var\bigg(\hat{\tau}(x)\bigg) + \bigg(Bias\big(\hat{\tau}(x)\big)\bigg)^2.\] Naively using the full data sample for estimation of both the nuisance functions as well as the CATE function leads to a higher bias due to overfitting but at the same time to lower variance as all available data is used for estimation. Using sample-splitting eliminates the overfitting bias but results in higher variance due to less data being used for estimation. In contrast, cross-fitting both removes the overfitting bias and reduces the variance by effectively using \ul{all the available information from the data} for estimation. Figure (ref) illustrates this theoretical argument by contrasting the distributions of the CATE parameter under full-sample estimation, double sample-splitting and double cross-fitting, \ul{resulting from a Monte Carlo simulation} based on a large training sample of \(32'000\) observations (further details on the meta-learner and the simulation design are provided in Sections (ref) and (ref), respectively). We observe that the \ul{theoretical} arguments can be documented in finite samples too. As such, the full sample version exhibits substantial bias due to overfitting as its distribution is shifted away from the true value of the CATE parameter, but with a rather low variance. On the contrary, the double sample-splitting version successfully eliminates the overfitting bias as the simulated distribution is centered around the true value of the CATE, however with much larger variance. Finally, the double cross-fitting version keeps the reduction in bias whilst having a much lower variance in comparison to the double sample-splitting version as the spread of the CATE distribution comes close to the full sample version, indicating the gain in efficiency of this procedure.
Apart from the illustrative example above, the empirical question remains the precise quantification of this bias-variance trade-off for various meta-learners and to what degree this might vary with different sample sizes. Different meta-learners use different nuisance functions in different ways which might have an influence on the performance under the particular estimation procedures. Even though sample-splitting and cross-fitting help to eliminate the overfitting bias, in finite samples less data available for estimation might even lead to higher bias due to errors in learning the CATE function itself, especially for small sample sizes. In this paper we address this open question via Monte Carlo simulations and compare the performance of various meta-learners under full-sample, double sample-splitting and double cross-fitting procedure for several different sample sizes to shed more light onto the finite sample properties. We follow Newey2018 and choose the double sample-splitting, respectively double cross-fitting procedure due to its theoretically faster convergence rates. Furthermore, we opt for the setting with equally sized \(K=M+1\) folds as suggested by Kennedy2020. Additionally, we always distinguish between the training and validation data. We use the training data for learning the nuisance function and the CATE function, including the double sample-splitting and double cross-fitting procedure, while we evaluate the CATEs on a set of new validation data. An illustration of the data usage under full-sample estimation, double sample-splitting and double cross-fitting procedure is provided in Figure (ref).
\ul{Further} motivation for the usage of sample-splitting and cross-fitting stems from the theoretical arguments for conducting statistical inference about the causal parameters of interest. As such, sample-splitting plays a crucial role in obtaining estimators that are not only approximately unbiased but also normally distributed which in turn allows for a valid construction of confidence intervals. In this vein, Chernozhukov2018c provide results for the estimators of average treatment effect (ATE) that rely on sample-splitting and cross-fitting procedures. Semenova2021 and Zimmert2019 extend this analysis for parametric and nonparametric estimators of group average treatment effects (GATEs), respectively. In the context of Causal Forests, Wager2018, Lechner2019, and Athey2019 also rely on sample-splitting procedures termed `honesty' to provide inference for causal effects on various levels of aggregation. Nonetheless, in the context of meta-learning estimation of causal effects, there \ul{appears to be lack of} unifying model-free theory for conducting statistical inference so far. One exception is the study by Kunzel2019 that analyses various versions of bootstrapping for estimation of standard errors for the CATEs. Recently, Jacob2021 makes use of such bootstrapping procedures to construct confidence intervals in an empirical application. Besides the computational burden, however, none of the bootstrapping procedures studied by Kunzel2019 seems to reliably provide accurate coverage rates. However, the meta-learners analyzed in Kunzel2019 do not make use of sample-splitting, which could potentially improve the performance of the bootstrapping for estimation of standard errors, given the insights from the related literature. While we do study the properties of the distribution of the CATEs within the simulation experiments in Section (ref), we do not further analyse the estimation of standard errors mainly due to computational reasons and focus primarily on the point estimators. However, apart from the computational aspects, we note that combining sample-splitting and cross-fitting with bootstrapping for statistical inference about causal effects within the meta-learning framework might be a promising avenue for future research.
\ul{In the following, we review the meta-learning algorithms for estimation of heterogeneous treatment effects and discuss their advantages and disadvantages in particular empirical settings.}
The first meta-learning algorithm we investigate is the S-learner as denoted by Kunzel2019. The early ideas of the S-learner stem from the marketing literature on uplift modelling Lo2002 and the epidemiology literature, where the S-learner is also being referred to as g-computation \parencites{Robins1986}{Snowden2011}. More recently, the S-learner has been proposed for the estimation of heterogeneous treatment effects based on the Regression Trees Athey2015 and Random Forests Lu2018 as well as their Bayesian versions, i.e. BART \parencites{Hill2011}{Green2012}. According to the naming convention of Kunzel2019, S- stands for Single as this meta-learner involves only one single model, namely the response function, \(\mu(x,w)\), that needs to be estimated. The final causal effect is, in this case, obtained as a difference \ul{between} predictions of the response function with setting the treatment indicator to $W_i=1$ and $W_i=0$, respectively.\footnote{Similar notion of predicted difference in the scope of model interpretation has been proposed by King2000.} The algorithm can be described as follows:\footnote{As a matter of notation, we refer to the training data used for model estimation with superscript \(^T\) as \(\{(X_i,Y_i,W_i)\}^T\) and the validation data used for effect prediction with superscript \(^V\) as \(\{(X_i)\}^V\).}\\ \\
\\ As can be seen from Algorithm (ref), the S-learner does not assign any special role to the treatment indicator \(W_i\) within the estimation procedure and uses it only post hoc in the computation of the causal effect. Thus, if the treatment indicator is not strongly predictive for the outcome the S-learner will tend to estimate a zero treatment effect.\footnote{Kunzel2019 argue that the S-learner is actually biased towards zero.} Nevertheless, the S-learner will perform particularly well if the true CATE function is indeed zero, i.e. if \(\tau(x)=0\), which has also been documented in the simulation experiments of Kunzel2019. For the forest based S-learner, Kunzel2019diss proposes a modification of the algorithm such that it shrinks towards the ATE instead of zero by performing a Ridge regression in the final leaves of the trees within the forest.\footnote{For a detailed explanation of this procedure see Kunzel2019diss.} In our simulations, we study a simpler modification of the forest based S-learner by always including the treatment indicator in the random subset of covariates when determining the splits. By doing so, we always give the S-learner the chance to split on the treatment indicator which might potentially alleviate the zero-bias issue. We will henceforth denote such learner as the SW-learner, where the W reflects the enforcement of the treatment indicator into the splitting set of covariates. We discuss the behaviour of the SW-learner more closely throughout the simulation results in Section (ref). Furthermore, notice that the Algorithm (ref) consists of only one nuisance function that needs to be estimated and thus does not require any sample-splitting or cross-fitting within the training sample induced by multiple nuisance functions, hence it always has access to the full sample of the training data.\footnote{Nevertheless, an optional additional sample-splitting or cross-fitting could potentially improve the performance of the S-learner by reducing the possible overfitting of the base learner as such. This is, however, beyond the scope of our analysis and is left for future research.} However, from a theoretical perspective, Alaa2018 show that modelling the CATE as a single response function as in the S-learner is not optimal in terms of achieving the fastest convergence rate as this depends on the complexity, defined as sparsity-to-smoothness ratio, of the response functions in the treated, \(\mu(x,1)\), and the control sample, \(\mu(x,0)\), and pooling these into a single response function enforces the complexity for treated and controls to be the same.
The T-learner is another common and widely used meta-learner that we investigate in our study. Similarly to the S-learner, its early applications emerged in the marketing literature on the uplift modelling \parencites{Hansotia2002}{Radcliffe2007}, while more recently it has been suggested for estimating heterogeneous treatment effects in the fields of medicine Foster2011 and econometrics Athey2015. In the \ul{econometric} literature it is sometimes also called as the basic Lechner2019, plug-in Kennedy2020 or naive Nie2021a CATE estimator. According to Kunzel2019, T- stands for Two as this meta-learner involves two models that need to be estimated, defined by the treatment indicator \(W_i\). These are namely the response function in the treated sample, \(\mu(x,1)\), and the response function in the control sample, \(\mu(x,0)\). This is in contrast to the above S-learner which pools the two response functions into a single one. However, similarly to the S-learner the causal effect is computed as a difference in predictions of the two response functions, \ul{which is motivated by the identification result as in Equation} (ref). The algorithm can be summarized as follows:\footnote{Notationwise, we refer to a subset of the data defined by a specific value of the variable as for example \(W_i=1\) by a subscript as \(\{(X_i,Y_i)\}_{W_i=1}\).}\\ \\
\\ Hence the T-learner uses the treatment indicator to split the estimation of the response function into two parts. This procedure is expected to work particularly well if the CATE function is complicated and there are no common trends in the response functions. This phenomenon finds supportive evidence in several simulation studies \parencites(see for example)(){Kunzel2019}{Jacob2020}{Curth2021}[or][]{Nie2021a}. Nonetheless, it is expected to work rather poorly if the CATE function is simple, as the response functions are not trained jointly and thus their difference might be unstable \parencites{Lechner2019}{Kennedy2020}{Nie2021a}. In terms of the estimation of the nuisance functions, the T-learner behaves similarly to the S-learner, as only the response functions need to be estimated to compute the CATE. As such no additional sample-splitting induced by multiple nuisance functions is required as the response functions are themselves estimated on separate samples defined by treated and control.\footnote{Again, this does not preclude that an optional sample-splitting or cross-fitting might be beneficial for the same reason as in the case of the S-learner Jacob2020. Using an honest forest as a base learner would also add an implicit sample-splitting procedure, however, this is not analysed herein.} Theoretically, the convergence rate of the base learners for estimating the response functions directly translates into the rate for estimating the CATE \parencites{Kunzel2019}{Kennedy2020} and thus depends on the complexity, i.e. the dimension and the smoothness, of the response functions Alaa2018.\footnote{For a detailed analysis of the optimal minimax rates for the T-learner, we refer to Alaa2018 and Kunzel2019.}
The above mentioned problems of the T-learner \ul{are aggravated} if the treatment assignment is \ul{highly} unbalanced, meaning that the vast majority of observations in the sample belongs to only one treatment status. Kunzel2019 therefore propose the X-learner which addresses this issue. The X-learner builds on the T-learner and, as such, first estimates the two response functions \(\mu(x,1)\) and \(\mu(x,0)\). It then uses these estimates to impute the unobserved individual treatment effects \ul{for the treated}, \(\tilde{\xi}_i^1\), and the control, \(\tilde{\xi}_i^0\). \ul{The imputed effects} are in turn used as pseudo-outcomes to estimate the treatment effects in the treated sample, \(\tau(x,1)\), and the control sample, \(\tau(x,0)\), respectively. The final CATE estimate \(\tau(x)\) is then a weighted average of these treatment effect estimates weighted by the propensity score, \(e(x)\).\footnote{In the original definition of the X-learner, the estimation of the propensity score is not exactly specified as it could be any weighting function in general. However, in practice the estimation of the propensity score is recommended Kunzel2019.} Thus the X-learner additionally uses the information from the treated to learn about the controls and vice-versa in a Cross regression style, hence the X term in its naming label. The learning algorithm can be detailed as follows:\\ \\
\\ According to Algorithm (ref), the X-learner, in contrast to the T-learner, firstly uses the response functions for imputing the unobserved individual treatment effects instead of directly estimating the CATE. Secondly, these imputed individual treatment effects are used for estimating the CATE and reweighted by the propensity scores. The reweighting helps to put more weight on the treatment effects which have been estimated more precisely, i.e. the ones coming from the larger treated or control sample, respectively. For this reason, the X-learner is expected to work particularly well in unbalanced settings, which is often the case in practice as the share of treated might be restricted financially or otherwise \parencites(see)(){Arceneaux2006}{Gerber2008}{Broockman2016}[or][for such unbalanced empirical settings]{Goller2020}. Furthermore, by directly estimating the treatment effects in the second step it enables the estimator to learn structural properties of the CATE function from the data and is thus expected to work well if the CATE function is approximately linear or sparse Kunzel2019. Theoretically, Kunzel2019 indeed prove that the X-learner achieves a faster convergence rate than the T-learner in cases where the treatment assignment is highly unbalanced or when the CATE function is linear. Even if no assumptions about the CATE function are imposed, the X-learner can be proven to achieve the same rate as the T-learner. This hinges on regularity conditions such as the Lipschitz continuity of the response functions.\footnote{More specifically, Kunzel2019 prove that the X-learner can indeed achieve the parametric $\sqrt{N}$ rate if the CATE function is estimated by OLS, while the nuisance functions can be estimated at any nonparametric rate. Furthermore, they prove that the X-learner achieves the minimax optimal rate if both the nuisance functions and the CATE function are estimated via $k$-NN regression.} In simulations of Kunzel2019 the X-learner performs reasonably well even in other non-favourable settings. Notice further that Algorithm (ref) requires more estimation steps than the previous two meta-learners. Additionally to the estimation of the response functions, the X-learner requires the estimation of the treatment effect functions as well as the propensity score function. This raises the question of possible overfitting and hence the need for sample-splitting and cross-fitting, respectively. However, there is theoretically no explicit requirement for sample-splitting in the case of the X-learner when estimating the nuisance functions, apart from training and validation data split Kunzel2019. Yet, it might well be that the sample-splitting and further cross-fitting have a non-negligible influence on the performance of the learner in finite samples. We address this issue by implementing the double sample-splitting and double cross-fitting version of the X-learner in the simulation study. For the case of the full-sample estimation we use the out-of-bag predictions of the underlying forest as estimates of the nuisance functions \ul{as recommended by} Lu2018 and Athey2019a. The out-of-bag predictions are based on the observations that have been left `out of the bag' when drawing bootstrap samples to estimate the trees of the forest Hastie2009. Such observations, however, randomly appear both as training as well as validation observations \ul{across the trees of the forest} and thus such out-of-bag predictions are neither the classical in-sample fitted values nor proper out-of-sample predictions \ul{and should not be confused with the `honest' predictions}. We use the out-of-bag predictions for all meta-learners within our analysis.
Although the X-learner makes use of the estimation of multiple nuisance functions, it does not provide the double robustness property which exploits the fact that the estimator remains consistent if either the response function or the propensity score function is misspecified \parencites{Kennedy2017}{Lee2017}. Recently, Kennedy2020 proposed the DR-learner where DR symbolizes the Double Robustness property of the learner. The DR-learner constructs a doubly robust score in the first estimation stage and estimates the CATE in the second stage. There have been many other versions of the DR-learner proposed in the literature, but these were restricted to a particular estimator used in the second stage and are thus not part of the meta-learning framework. For example, Semenova2021 propose a linear estimation of the CATE function, whereas a local-constant estimation is proposed by Zimmert2019 and Fan2020, which works well for the estimation of GATEs, i.e. for low-dimensional conditioning set. The main advantage of the DR-learner in comparison to the other versions lies in the general model-free second stage \ul{estimation of the CATE function}. However, common to all versions in the literature is the estimation of the doubly robust score\footnote{Also called efficient score or efficient influence function in the literature \parencites{Robins1995}{Hahn1998}.} by machine learning methods in the first stage also known as Double Machine Learning Chernozhukov2018c. For a comprehensive overview of the CATE estimators building on the doubly robust score see Knaus2020b. The specific algorithm for the DR-learner is then defined as follows:\\ \\
\\ As can be seen in Algorithm (ref) above, the DR-learner estimates the very same nuisance functions, \(\mu(x,0), \mu(x,1)\) and \(e(x)\), as the X-learner but uses them in a completely different manner. It combines the nuisance functions \ul{as well as} the outcome and treatment data in a doubly robust way to construct the pseudo-outcome \(\psi_i\), i.e. the doubly robust score. The score is then regressed on the covariates to estimate the final CATE function. Therefore, the DR-learner can also adapt to structural properties of the CATE such as smoothness or sparsity. For this reason the DR-learner is expected to work well in similar situations as the X-learner with a more balanced treatment assignment, as too extreme propensity scores might possibly yield the estimator unstable \parencites{Huber2013}{Powers2018}, especially in high dimensions DAmour2021. Moreover, it should have an additional advantage over the X-learner thanks to its double robustness property. The theoretical analysis of Kennedy2020 uses the double sample-splitting procedure in order to derive a sharp error bound that rests only on stability conditions for the estimation of the CATE function. The theoretical results then exploit the rate double robustness which allows for faster error rates for the second-stage CATE estimation under weaker rate conditions for the first-stage estimation of the nuisance functions. The simulations of Kennedy2020 also suggest a faster convergence rate of the DR-learner in comparison to the X- and T-learner. In order to achieve the optimal rates\footnote{Kennedy2020 shows that the DR-learner achieves the minimax optimal rate under smoothness or sparsity conditions for the nuisance and the CATE functions.} the DR-learner explicitly requires the double sample-splitting as defined by Newey2018, while the double cross-fitting procedure remains optional. Theoretically it is not clear how important the role of the optional cross-fitting is for the DR-learner in finite samples and how much of the efficiency loss due to sample-splitting can be thereby regained. In order to shed light on this issue we investigate the implementations of the DR-learner with double sample-splitting, double cross-fitting, as well as a version with full-sample estimation.
Yet another approach of first estimating nuisance functions and then using them to learn the treatment effects stems from the literature on \ul{partially linear model} originally developed by Robinson1988. Nie2021a build on these ideas to flexibly estimate heterogeneous treatment effects and develop the R-learner, where the R stands for the \ul{recognition of the contribution of} Robinson1988 as well as for the Residualization \ul{approach}. In the first step, the R-learner estimates the response function, \(\mu(x)\), similarly to the S-learner but without conditioning on the treatment indicator, as well as the propensity score function \(e(x)\). It then residualizes the outcome and the treatment by the predictions of the response and the propensity score function, respectively, to construct a modified outcome. In the second step, the R-learner regresses the modified outcome on the covariates, weighted by the squared residualized treatment\footnote{An estimation procedure without the weighting step is in literature referred to as the U-learner \parencites{Stadie2018}{Kunzel2019}{Nie2021a}. However, such estimator turned out to be quite unstable in the simulation experiments in Nie2021a as well as in those of Kunzel2019 and will thus not be considered further in our analysis.}, i.e. \(\big(W_i-\hat{e}(X_i)\big)^2\), to estimate the CATE function Schuler2018. Analogous transformation of the outcome is also used by the Causal Forest of Athey2019 termed local centering, or in the G-estimation for sequential trials by Robins2004. The full estimation procedure of the R-learner can be summarized as follows:\\ \\
\\ As follows from Algorithm (ref), the R-learner separates the estimation into two steps. First, it eliminates the spurious correlations between the response function \(\mu(x)\) and the propensity score function \(e(x)\) and second, it optimizes the CATE function \(\tau_R(x)\). From this standpoint the R-learner follows a related estimation scheme as the DR-learner and is expected to work well in similar settings where the nuisance functions as well as the CATE function might have \ul{a high degree of} complexity. A possible advantage of the R-learner over the DR-learner might stem from the additional weighting which reduces the impact of extreme propensity scores as pointed out by Jacob2021. In their simulation experiments, Nie2021a show good performance of the R-learner in settings with complicated nuisance functions and rather simple CATE function. The theoretical analysis of Nie2021a shows that the convergence rate of the R-learner depends on the complexity of the CATE function and is not affected by the complexity of the nuisance functions as long as they are estimated at sufficiently fast rates.\footnote{For the case of estimating the CATE function via penalized kernel regression, Nie2021a prove that the R-learner achieves the minimax optimal rate and additionally show that the X-learner does not achieve the optimal rate in general unless the treatment assignment is not highly unbalanced.} Furthermore, for these theoretical results Nie2021a explicitly require sample-splitting and cross-fitting, respectively. In particular, they advocate for a 5- or 10-fold cross-fitting procedure as defined by Chernozhukov2018c. In order to examine the importance of the cross-fitting in finite samples we compare the performance of the R-learner as in the above cases with full-sample estimation, double sample-splitting and double cross-fitting, respectively.
We study the finite sample performance of meta-learners for estimation of heterogeneous treatment effects based on Random Forests \parencites()(){Breiman2001}[see also][for a comprehensive introduction]{Biau2016}. The focus of the Monte Carlo study lies in an assessment of the influence of sample-splitting and cross-fitting in the causal effect estimation. For this purpose we compare the above discussed meta-learners estimated with full-sample, double sample-splitting, and double cross-fitting.
We rely on the Random Forest as the base learner for all meta-learners for several reasons. First, different meta-learners require estimation of different nuisance functions which involve different types of outcome variables. As such, the response functions mostly involve a continuous outcome variable whereas the propensity score function includes a binary outcome. Hence, when using Random Forests no additional adjustments need to be done in terms of estimation as it automatically estimates probabilities in case of binary outcome and expected values in case of continuous outcomes, respectively. This is in contrast to linear learners such as the Lasso Tibshirani1996, Ridge Hoerl1970 or Elastic Net Zou2005 where the estimator needs to be modified using appropriate link function for proper probability estimation Hastie2009. Second, Random Forest is a local nonparametric method which does not need any data pre-processing to flexibly learn the underlying functional form from the data Hastie2009. Thus, Random Forest is able to approximate any function with different degrees of complexity which is often the case in treatment effect estimation where the nuisance functions tend to be rather difficult complex functions while the CATE function itself is often \ul{argued to be} simple and sparse \parencites{Kunzel2019}{Kennedy2020}{Sekhon2021}. This is again an advantage in comparison to the linear learners mentioned above which become more flexible once an augmented covariate set consisting of polynomials and interactions is created and thus can be regarded as global nonparametric methods Hastie2009. Third, in contrast to other flexible state-of-the-art machine learners such as Neural Networks the theoretical properties of Random Forests are better understood which makes it less of a black-box method Meinshausen2006, Biau2012, Wager2014, Wager2014a, Scornet2015, Mentch2016, Wager2018, Athey2019. In particular, Scornet2015 prove the consistency of the original Random Forest algorithm as developed by Breiman2001 that we employ in our simulations, without relying on the `honesty' condition.\footnote{The consistency (in MSE) result holds under the condition that the number of leaves is smaller than the number of observations. In our simulations, we ensure this by growing trees with minimum leaf size bigger than 1 (see Table (ref)). We refrain from the honesty feature to avoid additional sample-splitting that would further reduce the effective sample size.} \ul{This is important as the consistency of the base learner is a fundamental condition shared across all considered meta-learners.} Scornet2015 \ul{also show that Random Forests can effectively adapt to sparsity of the underlying model, which we explicitly make use of within our simulation design.} Additionally, another reason why we do not use the Lasso and linear learners as such is due to a substantial increase in variance as they are prone to outliers as documented in the simulation studies of Jacob2020 as well as Knaus2021. Furthermore, Random Forests are a popular choice as a base-learner in empirical studies using meta-learners too \parencites{Duan2019}{Knaus2020b}. Lastly, from the practical standpoint there is a vast variety of fast and reliable software implementations of Random Forests which makes it easy to use for practitioners.\footnote{In our simulations we use the R-package ranger which provides a fast C{\nolinebreak[4]\raisebox{.4ex}{\tiny\bf ++ }} implementation of Random Forests, particularly suited for high-dimensional data Wright2017. Further options in the R language rstats include the grf package written by grf, the forestry package by forestry or the randomForest package by randomForest.}
In order to objectively evaluate the performance and the robustness of different meta-learners in estimating heterogeneous treatment effects with regard to the double sample-splitting and double cross-fitting, we design several simulation scenarios. On the one hand, for each meta-learner we construct such a data generating process (DGP) that \ul{suits the particular advantages} of the \ul{respective} meta-learner, i.e. we design a simulation scenario where each meta-learner is expected to work best. Hence, we are able to check if the particular meta-learner outperforms the others and how big the performance discrepancies are for the other meta-learners in comparison to the expected best performing meta-learner. On the other hand, we design \ul{a simulation scenario with a complex DGP} where none of the meta-learners has a priori an explicit \ul{advantage}, \ul{which serves as our main simulation design of interest}. Thus we can compare the performance of the meta-learners in an objective manner and quantify the deviations to their respective best performance cases. Furthermore, common to all DGPs is the observational study design, i.e. there is always selection into treatment and thus all considered meta-learners have to deal with confounding and not only with modelling the treatment effect itself. Moreover, in contrast to many simulation studies where the nuisances are simple low-dimensional functions \parencites{Wager2018}{Kunzel2019}{Kennedy2020}, we model all nuisance functions as highly non-linear but sparse functions with large-dimensional covariate space to test the potential of the machine learning methods, though still largely obeying the theory induced limitations. For other \ul{complex} simulation designs see also Jacob2020 or Zivich2021 as well as Lechner2019 and Knaus2021 for the Empirical Monte Carlo Studies. \ul{Importantly,} in order to study the \ul{approximate convergence rates} of the meta-learners we repeat each simulation scenario several times with increasing training sample sizes using \(N^{T}=\{ 500, 2'000, 8'000, 32'000 \}\). \ul{We emphasize that the considered sample sizes substantially exceed the ones from previous simulation studies devoted to the analysis of sample-splitting methods, which were limited to} $2'000$ Jacob2020 and $3'000$ Zivich2021 \ul{observations, respectively. Furthermore, the large samples enable us to study the performance of the meta-learners in settings in which the application of machine learning methods is arguably more relevant.} We choose to always quadruple the sample size, which allows us to easily benchmark the results with the parametric \(\sqrt{N}\) rate, \ul{in which case the estimation error is expected to halve with each increase of the sample size.} We then evaluate the performance measures on a validation set with sample size of \(N^V=10'000\) to reduce the prediction noise as is usual in many Monte Carlo studies \parencites{Janitza2016}{Hornung2017}{Okasa2019}{Jacob2020}{Knaus2021}. Lastly, in terms of the tuning parameters for the Random Forest base-learner we stick to the default, in the literature commonly used settings, corresponding to \(1'000\) trees, the number of randomly chosen split variables set to the square root of number of features, and the minimum leaf size equal to \(5\).\footnote{We refrain from cross-validation or other tuning parameter optimization procedures due to computational constraints. We recommend such optimization in the applied work as it might considerably improve the performance of the estimator Curth2021, however, for the purposes of the simulation study it would not change the relative ranking of the meta-learners as each of them uses the very same base learner.} Finally, for each DGP we simulate the training data \(R=\{2'000, 1'000, 500, 250\}\) times in total, where we use \(2'000\) replications for the smallest sample size and decrease the number of replications down to \(250\) for the largest sample size, due to computational reasons.\footnote{Notice, however, that we only halve the number of replications while quadrupling the sample size and as such we may \ul{limit} a possible deterioration of the performance in terms of the simulation error. A similar strategy for balancing the precision and the computational burden has been used in the simulations by Lechner2019, Lu2018 or Knaus2021. Detailed results on the simulation error are provided in Appendix (ref).}
For the evaluation of the performance of the considered meta-learners with regard to the sample-splitting and cross-fitting in detail, we employ several evaluation measures. First, to assess the overall estimator performance we compute the root mean squared error for each observation \(i\) from the validation sample over the \(R\) simulation replications:\footnote{We take the square root of the MSE to have the same \ul{scale} as for the other performance measures, i.e. the absolute bias and the standard deviation.} \[RMSE \big(\hat{\tau}(X_i)\big)= \sqrt{\frac{1}{R} \sum^R_{r=1} \big( \tau(X_i) - \hat{\tau}^r(X_i) \big)^2}.\] Next, we decompose the root mean squared error and evaluate the bias and variance component separately to contrast the theoretically expected asymptotic behaviour of sample-splitting and cross-fitting with the finite sample properties. Hence, we additionally compute the mean absolute bias: \[\mid BIAS \big(\hat{\tau}(X_i)\big)\mid= \frac{1}{R} \sum^R_{r=1} \mid \tau(X_i) - \hat{\tau}^r(X_i) \mid\] as well as the standard deviation of the treatment effects: \[SD \big(\hat{\tau}(X_i)\big)= \sqrt{\frac{1}{R} \sum^R_{r=1} \bigg( \hat{\tau}^r(X_i) - \frac{1}{R} \sum^R_{r=1} \hat{\tau}^r(X_i) \bigg)^2}.\] Furthermore, inspired by the simulation study of Knaus2021 we also compute the Jarque-Bera statistic \parencites{Jarque1980}{Bera1981} to test for the \ul{normality} of the treatment effect predictions:\footnote{See Thadewald2007 \ul{for a discussion of the Jarque-Bera test and its comparison to other tests for normality.}} \[JB \big(\hat{\tau}(X_i)\big)= \frac{R}{6}\bigg( S\big(\hat{\tau}(X_i)\big)^2 + \frac{1}{4} (K\big(\hat{\tau}(X_i)\big)-3)^2 \bigg)\] where \(S\big(\hat{\tau}(X_i)\big)\) and \(K\big(\hat{\tau}(X_i)\big)\) is the skewness and the kurtosis of the \(R\) treatment effect predictions \ul{for observation} \(i\), respectively. As a matter of presentation for CATEs, we report the mean values of the RMSE, absolute bias, standard deviation and the Jarque-Bera statistic over \(N^V\) validation observations.\footnote{\ul{As such, we define the average RMSE as} \(\overline{RMSE}=\frac{1}{N^V}\sum^{N^V}_{i=1} RMSE \big(\hat{\tau}(X_i)\big)\) \ul{and analogously for the remaining performance measures. Additionally, for the Jarque-Bera statistic we report also the share of CATEs from the validation sample for which the normality gets rejected at the 5% confidence level. For details, see Appendix} (ref).} \ul{Additionally, we provide evaluation of further performance measures in Appendix} (ref).
In the general simulation design we follow Kunzel2019 and specify the response functions \ul{for potential outcomes} under treatment, \(\mu_1(x)\), and control, \(\mu_0(x)\), the propensity score, \(e(x)\), and the treatment effect function, \(\tau(x)\), respectively. First, we simulate a \(p\)-dimensional matrix of covariates \(X_i \in \mathbb{R}^p\) drawing from the uniform distribution, as previously used in simulations of Wager2018, Kunzel2019 or Nie2021a among others, such that: \[X_i \sim \mathcal{U}\big([0,1]^{n\times p}\big)\] and defining the correlation structure according to Falk1999 using a random correlation matrix \(\Sigma_p\) generated by the method of Joe2006.\footnote{For a detailed correlation heat map of the covariates and descriptive statistics of the simulated datasets see Appendix (ref).} Second, we specify the response functions and simulate the potential outcomes as: \[Y_i(0)=\mu_0(X_i)+\epsilon_i(0)\] \[Y_i(1)=\mu_1(X_i)+\epsilon_i(1)\] with errors \(\epsilon_i(0), \epsilon_i(1) \stackrel{iid}{\sim} \mathcal{N}(0,1)\) that are independent of the covariates \(X_i\). Third, we define the propensity score function and simulate the treatment assignment according to: \[W_i \sim Bern\big(e(X_i)\big)\] and use the observational rule to set the observed outcomes such that: \[Y_i = W_i \cdot Y_i(1) + (1-W_i) \cdot Y_i(0)\] to complete the observable triple \(\{(X_i,Y_i,W_i)\}\). The subsequent simulation designs then differ only with respect to how the corresponding functions, namely \(\mu_0(x), \mu_1(x), e(x)\) and \(\tau(x)\) are specified. For all of our simulations we define the response function under non-treatment according to the well-known Friedman function (Friedman1991) to create a difficult yet standardized setting, which has been used also in the simulations of Biau2012 and Nie2021a, as follows:
hence effectively creating a highly non-linear but sparse response function which is difficult to estimate on its own.\footnote{Note that $\pi$ refers to the mathematical constant, i.e. $\pi \approx 3.14$.} The response function under treatment is then defined simply as: \[\mu_1(x)= \mu_0(x)+\tau(x)\] while we vary the specification of the treatment effect function \(\tau(x)\) throughout our simulation designs. Lastly, we model the propensity score function similarly to Wager2018 and Kunzel2019 using the \(\beta\) distribution with parameters \(2\) and \(4\) such that:
while the scale parameter \(\alpha\) controls the share of treated in the sample and at the same time helps to bound the resulting probabilities away from \(0\) and \(1\) and thus to avoid extreme propensity scores which might yield some meta-learners using such propensities for reweighting unstable \parencites{Huber2013}{Powers2018}. We additionally make the propensity score dependent on features \(X_i\) of dimension \(p^e\) in a non-linear fashion using the functional form specification of Nie2021a and set: \[f(x)=sin(\pi\cdot x_1 \cdot x_2 \cdot x_3 \cdot x_4)\] which creates a non-linear setting that is hard to model as opposed to, e.g. polynomial transformations alone. Similarly, such non-linear transformations for the propensity scores using the sine function have been used also in simulations by Lechner2019 and Knaus2021.
As a matter of notation we refer to \(p\) as the dimension of the covariate space, \(p^{\mu}\), \(p^{e}\) and \(p^{\tau}\) as the dimension of the signal covariates in the response function, the propensity score function, and the CATE function, respectively. We set the aforementioned dimensions as follows: \(p=100\), \(p^{\mu}=5\), \(p^e=4\) and \(p^{\tau}\) varies with forthcoming simulation designs. We note that such large-dimensional covariate set is quite unique as the majority of simulation studies relies on low-dimensional covariate sets \parencites(see e.g.)(){Kunzel2019}{Jacob2020}[or][]{Nie2021a}.\footnote{An exception is the simulation study of Powers2018 \ul{who explicitly study the estimation of heterogeneous treatment effects in high-dimensions.}} We further define the sets of covariates such that \(X^{p^{\tau}} \subset X^{p^{e}} \subset X^{p^{\mu}} \subset X^{p}\). By doing so we make it difficult for the meta-learners to accurately fit the functions and eliminate the spurious correlations between the response and propensity score functions. Moreover, it also becomes a non-trivial task to disentangle the confounding effects from the actual treatment effect heterogeneity which the herein discussed meta-learners are specifically designed for. Finally, a general overview of the simulation study is provided in Table (ref).
The first simulation design features our complicated sparse non-linear nuisance functions as defined above in Equations (ref) and (ref) in contrast to a very simple CATE function. In fact, the treatment effect here is defined as being constant and equal to zero: \[\tau(x)=0\] with a balanced treatment assignment with the scaling factor \(\alpha=\frac{1}{4}\) which results in approximately 50% treated and 50% of control units. Such DGP with zero CATE serves as a benchmark and should implicitly suit the S-learner as the treatment indicator is not predictive for the outcome. Nevertheless the other meta-learners with the exception of the T-learner should be also capable of capturing the true zero effect as this is often a showcase example when motivating the particular meta-learners as well as simulating their performance \parencites(see)(){Kunzel2019}{Kennedy2020}[and][for details]{Nie2021a}.
In the second simulation design we keep the balanced treatment allocation but feature a highly complex and non-linear CATE function resulting from a completely disjoint DGPs for the response function under treatment and under control. As such the response function under control is defined according to Equation (ref), while the response function under treatment is defined as a non-zero constant, i.e. \(\mu_1(x) = 1\). The CATE is then defined as: \[\tau(x)=\mu_1(x)-\mu_0(x).\] Such simulation setups have been previously used also in Kunzel2019 or in Nie2021a. In this case the response functions, \(\mu_0(x)\) and \(\mu_1(x)\), are uncorrelated and thus there is no advantage in pooling those two together. Rather, estimating these two functions separately is the best strategy as there is nothing additional to learn from the other treatment group. For this reason, the T-learner should perform best here, however the meta-learners which also estimate the response functions separately such as the X- and DR-learner are expected to perform well too. Clearly, other meta-learners such as the S- and R-learner which estimate the pooled response function have a disadvantage as they first need to learn the disjoint structure.
In our third simulation design we change the scaling factor in the propensity score function to \(\alpha=\frac{1}{12}\) such that we generate approximately 15% treated units.\footnote{In contrast to Kunzel2019 we do not specify the treatment imbalance as extreme as 1% treated mostly for computational reasons. Due to our smallest sample size of \(N=500\) used in the simulations and the double sample-splitting procedure, it might occasionally happen that the estimated propensity scores would be exactly zero which would prevent estimation of the DR-learner as well as the R-learner due to the division by zero when constructing the pseudo-outcomes. In our specification, even with the share of the treated being \(16.77\%\) in expectation, the aforementioned issue with zero propensity scores still might occur. In such cases, we redraw the sample to ensure at least \(15\%\) of treated. However, this occurs only a handful of times out of 2000 draws in total and only for the smallest sample size considered. Nie2021a also use similar restrictions on the propensity scores in their simulations due to the very same issue.} We then model the treatment effect as a constant as for example in Kennedy2020 or Nie2021a and thus create a scenario with highly complicated nuisance functions and very simple CATE function given as: \[\tau(x) = 1.\] Accordingly, the X-learner should perform best in this scenario given the high imbalance in the treatment assignment and the sparse CATE function at the same time. In contrast, other meta-learners using the propensity score for reweighting such as the DR- and R-learner might perform worse due to potentially extreme propensity scores close to the \(\{0,1\}\) bounds. Furthermore, the T-learner is clearly disadvantaged in this scenario due to the high treatment imbalance as well as due to the simple CATE function, whereas the S-learner is not expected to work particularly well either due to the relatively bigger effect size bounded away from zero.
In our fourth simulation design we model the CATE function similarly to the above design as a simple non-zero constant and combine it with an indicator function as also used by Kunzel2019 to add more structure to the CATE. As such we define the treatment effect as: \[\tau(x) = 1 + 1 \cdot\mathds{1}(x_1 > 0.5)\] and otherwise keep the DGP same as in the third design while only increasing the share of treated to roughly 25% as is the case in the simulations of Nie2021a by setting \(\alpha=\frac{1}{8}\). By doing so we should theoretically shift the advantage from the X-learner more onto the DR-learner as both meta-learners are motivated by complex nuisance functions and a simple CATE function with the difference of the X-learner being designed particularly for highly unbalanced treatment allocation. Also the R-learner is expected to perform relatively well in this scenario due to the less unbalanced treatment shares, whereas the S- and T-learner are not expected to perform well for the same reasons as in the above situations.
The fifth simulation design features the same nuisance functions and treatment share as the fourth design, however, here instead of the indicator function we model the treatment effect as a low-dimensional linear function as: \[\tau(x)=1+\frac{1}{2}x_1+\frac{1}{2}x_2\] as used in one of the simulation designs of Nie2021a where the R-learner performed best and as such it is also expected to have an advantage here. Yet again the X- and DR-learner should perform comparatively well in this setting while the S- and T-learner not so much for the very same reasons as stated above.
In the last simulation design we create the most complex scenario in which none of the meta-learners has an a priori advantage and thus presents our main simulation design of interest. In this case not only the nuisances but also the CATE itself is modelled as a smooth non-linear function of a slightly larger dimension than in the previous settings, i.e. \(p^{\tau}=3\). Following Wager2018 we specify the CATE function as follows: \[\tau(x)=1+\frac{4}{p^{\tau}}\sum^{p^{\tau}}_{j=1}\bigg(\frac{1}{1+e^{-12(x_j-0.5)}}-\frac{1}{2}\bigg).\] We further keep the treatment share equal to about 25% and the nuisance functions as previously specified as well. Hence, the meta-learners need to first account for the moderately unbalanced treatment shares, second accurately estimate the complex nuisance functions and disentangle their correlation, and third separate the treatment effect heterogeneity from the selection effects by precisely estimating the non-linear CATE function.
For the analysis of the simulation results we focus on the Main Simulation design as this is the most complex simulation design which does not a priori create conditions that would be advantageous for any of the considered meta-learners. \ul{Additionally, we argue that intuitively the performance in the most complex setting might be generalizable for simpler settings, too.} We then summarize the results for the rest of the simulation designs for which we provide the detailed results in Appendix (ref). Supplementary results providing additional measures, including the simulation error, bias, skewness, kurtosis, \ul{share of CATEs for which the normality has been rejected,} as well as the correlation and variance ratio of the estimated and the true CATEs are presented in Appendix (ref).
The performance of the meta-learners in the Main Simulation is depicted in Table (ref). We report the results for the average values of the RMSE, absolute bias, standard deviation and the Jarque-Bera test statistic over the \(N^V=10'000\) predicted CATEs from the validation sample. Figure (ref) details the performance of the meta-learners \ul{implemented} in the full-sample, double sample-splitting and double cross-fitting versions.
Starting with the most simple S-learner, we observe a competitive performance in terms of the average RMSE for the smaller sample sizes which, however, disappears for larger sample sizes. Taking a closer look at the results reveals that the competitive performance of the S-learner stems mainly from the very low standard deviation while being substantially biased. Indeed, the variance of the S-learner is the smallest among all meta-learners for all sample sizes. This is mainly due to its tendency to predict effects close to zero if the treatment indicator is not strongly predictive for the outcome as pointed out by Kunzel2019. This explains also the high bias of this estimator as the simulated CATEs vary between \(-1\) and \(3\) with only a small proportions of these CATEs being equal to zero (see Figure (ref) in Appendix (ref) for details). \ul{Similar pattern with high bias of the S-learner with the Random Forest as a base learner is documented also in the simulations of} Lu2018. Nevertheless, the Jarque-Bera test does not indicate evidence against the normality of the predicted CATE distribution, on average.
Considering the modified version of the S-learner with enforcement of the treatment indicator into the forest splitting set, i.e. the SW-learner, we notice almost identical performance to the one of the T-learner. This result can be explained by an observation that once the SW-learner finds the split based on the treatment indicator early within the trees it mimics the disjoint structure of the T-learner. The rest of the recursive partitioning is then very similar to the one of the T-learner which has been also documented for the case of the S-learner in the simulation experiments conducted by Kunzel2019. As a result, it seems that enforcing the treatment indicator into the splitting set helps to alleviate the high bias of the S-learner to some degree, however, it increases the variance of the estimator at the same time. Nevertheless, the bias-variance trade-off in this case results in lower average RMSE in comparison to the S-learner and the SW-learner might thus be preferred over the simple S-learner, when using the Random Forest as a base learner. Overall, the SW- and T-learner are very competitive in the smaller sample sizes both in terms of the average RMSE as well as the average absolute bias. However, with access to more training data these two learners do not improve fast enough and are outperformed by the more sophisticated learners in the largest sample size \ul{consisting of $32'000$ observations}. Concerning the distribution of the predicted CATEs there seems to be on average no statistical evidence against the normality, neither for the SW-learner nor for the T-learner.
In contrast to the above \ul{mentioned} meta-learners the X-learner makes use of the additional estimation of nuisance functions. In its full-sample version the X-learner performs best in terms of the average RMSE for all sample sizes, except the largest one. The good RMSE performance stems partly from the relatively low bias and partly from the relatively low variance of the estimator as the X-learner exhibits the smallest average absolute bias for the smaller sample sizes ($500$ and $2'000$), while having one of the lowest average standard deviations throughout all sample sizes. Interestingly, we only partly document the theoretical properties regarding the sample-splitting and cross-fitting procedures. As such, the full-sample version is the best performing one in terms of the average RMSE as well as in terms of the average absolute bias across all sample sizes, which is a pattern observed in the simulation experiments of Jacob2020 as well. Accordingly, the sample-splitting version exhibits not only higher values of the average standard deviation but also higher values of the average absolute bias. Nevertheless, we observe that the cross-fitting version successfully regains the efficiency lost due to sample-splitting as it exhibits steadily lower variance than the sample-splitting version and even lower variance than the full-sample version, while having a bias of roughly the same magnitude as the sample-splitting version. As for the distribution of the predicted CATEs, on average, we do not observe evidence for deviations from normality for any of the versions of the X-learner. Additionally, we do not observe any major differences in the speed of convergence between the different versions as can be seen in Figure (ref). Moreover, the absolute differences in the performance measures among the different versions are small in comparison to other meta-learners using nuisance functions. Albeit rather surprising at the first sight, the explanation for this phenomenon comes presumably from the different usage of the propensity score by the X-learner in comparison to the R- and DR-learner. As such, the R- and DR-learner use the propensity score \ul{together with the response functions} to construct a new pseudo-outcome which is subsequently used to estimate the CATEs. In contrast, the X-learner uses \ul{merely the response functions to create the pseudo-outcome, while} the propensity score \ul{is used only} to reweight the final CATE estimates and thus it does not enter into any additional estimation step. Therefore, \ul{the X-learner might be less prone to overfitting bias} which would \ul{partly} justify the full-sample estimation as described in Kunzel2019.\footnote{Nonetheless, this insight might still substantiate the need for sample-splitting, although only with two folds instead of three as used here.}
Assessing the performance of the DR-learner reveals some interesting insights. The first observation is that the cross-fitting version performs best of all meta-learners in terms of the average RMSE for the largest sample size of $32'000$ observations. This comes mainly from the low bias of this estimator as the average absolute bias is the lowest among all learners for the two largest sample sizes, while the average standard deviation is relatively high. However, looking at \ul{the average value of} the Jarque-Bera statistic suggests \ul{evidence against the} normality of the predicted CATEs for all but the largest sample size. Inspecting the results more closely reveals that the issue stems from heavy tails of the CATE distributions. The extreme values of the predicted CATEs are mainly caused by the propensity scores which are close to the \(\{0,1\}\) bounds. Similar issues of the DR-learner due to extreme propensity scores have also been documented in the simulation experiments of Knaus2021 as well as in the empirical application of Knaus2020b. The second observation is that for the DR-learner we clearly see how the theoretical arguments of sample-splitting and cross-fitting translate into the finite sample properties of the estimator. However, these can be documented only for large sample sizes. As such, the bias of the sample-splitting version is smaller than the one of the full-sample version in the largest sample size, while the bias of the cross-fitted version is even slightly lower than the sample-splitting version and is lower than the bias in the full-sample version for both the largest \ul{($32'000$)} and the second largest \ul{($8'000$)} sample considered. For the smaller sample sizes \ul{($500$ and $2'000$)} we see that the reduction in the overfitting bias is not large enough in comparison to the bias stemming from the estimation of the CATE function. As such, for small sample sizes the additional splitting of the sample does not leave enough observations to learn the non-linear structure of the CATE. Considering the variance of the estimator, we also observe the theoretically expected pattern. The full-sample version of the DR-learner exhibits the smallest average standard deviation throughout all sample sizes, while the standard deviation for the sample-splitting version is roughly twice as high. Nevertheless, the cross-fitting version successfully reduces the standard deviation for all sample sizes and comes close to the full-sample version, effectively regaining the lost efficiency of the estimator due to sample-splitting. Overall, in terms of the average RMSE this bias-variance trade-off results in favourable performance of the sample-splitting version in the largest and of the cross-fitting version in the two largest samples in comparison to the full-sample version. Considering the distribution of the predicted CATEs we see that the heavy tails problem due to extreme propensity scores is the worst for the sample-splitting version, where even in the second largest sample size of $8'000$ observations, the normality is rejected for more than 50% of the predicted CATEs from the validation sample (compare the supplementary results in Appendix (ref)). This stems from the smaller samples used for estimation of the propensity scores which are more likely to yield extreme values under unbalanced treatment assignment. We also observe that this issue is less pronounced for the full-sample version. The third and the last observation is yet the probably most noticeable pattern across all performance measures, namely the fast convergence of the sample-splitting and cross-fitting version of the DR-learner which is substantially faster in comparison to all other meta-learners as can be seen in Figure (ref). As such the DR-learner performs almost worst of all, both in terms of the average RMSE and average absolute bias for the smallest sample size of $500$ observations, but almost best of all for the largest sample size of $32'000$ observations. This provides evidence that the DR-learner is able to learn a \ul{highly} complex CATE function \ul{once enough data becomes available} and \ul{additionally} highlights the need for sample-splitting and cross-fitting in order to achieve the \ul{theoretically described optimal performance} Kennedy2020.
The performance of the R-learner is competitive with the other meta-learners especially in smaller samples, particularly for the full-sample version. In the smallest sample size of $500$ observations the R-learner outperforms the DR-learner in terms of the average RMSE, irrespective of the estimation procedure. However, with growing sample sizes the performance evens out and eventually for the largest sample size of $32'000$ observations the R-learner \ul{lags behind the majority of the estimators}. This is in contrast to previous results from simulations of Jacob2020 and Knaus2021 where the R-learner exhibits rather good performance, albeit studied only in smaller samples. The decomposition of the RMSE shows that while the full-sample version of the R-learner exhibits rather low bias, it suffers from a higher variance as can be seen in Figure (ref). Nonetheless, the distributions of the predicted CATEs do not show \ul{on average} deviations from the normal distribution. This is contrary to the DR-learner and illustrates the advantage of the additional weighting step. As such, even though the R-learner uses the propensity scores for reweighting to construct the modified outcome, it successfully manages to downweight the modified outcomes based on extreme propensity scores to alleviate the heavy tails issues observed in the case of the DR-learner. In particular, even for the sample-splitting version of the R-learner the share of predicted CATEs for which the normality is rejected is an order of magnitude lower in comparison to the DR-learner (see Appendix (ref) for details). In terms of the estimation procedure, we observe a similar pattern as for the X-learner in a sense that the full-sample version performs better with respect to the average RMSE and absolute bias, while the cross-fitting version helps to reduce the variance of the estimator not only in comparison to the sample-splitting version but even in comparison to the full-sample version. The sample-splitting version exhibits higher values of the average absolute bias and standard deviation for all sample sizes considered, while the convergence rates are approximately same as for the full-sample and the cross-fitting version. Hence, there is \ul{a lack of} indication that the overfitting type of bias reduction could become relevant in bigger samples. Similarly to the DR-learner, also for the R-learner the differences between the different estimation procedures seem to be higher than those for the X-learner which is again presumably due to the different usage of the propensity scores.
Inspecting the results for the rest of the simulation designs reveals further insights and helps to generalize the findings from the main and most complex simulation design discussed \ul{above}.
Within the benchmark Simulation 1 with zero constant CATE the S-learner, as expected, performs best with respect to all performance measures across all sample sizes (see Table (ref) in Appendix (ref)). However, the results reveal poor statistical properties of the S-learner as it appears to be substantially biased and inconsistent as the absolute bias as well as standard deviation increase with growing sample size.\footnote{A closer look on the estimation results reveals the reason for this phenomenon. With small sample sizes, the underlying trees of the S-learner's forest are quite shallow and barely split on the treatment indicator resulting in quite homogeneous CATE predictions which are very close to the actual zero effect. However, as the sample size increases, the chance of splits based on the treatment indicator increases which results in more heterogeneous effect predictions spread around zero. Accordingly, the bias as well as the standard deviation increase. Similar consistency issues of the forest-based S-learner seem to appear also in the simulations of Kunzel2019 where the MSE rises with growing sample size for some designs and only stabilizes with very big sample sizes.} In general, the performance of the S-learner is, in all simulation designs, plagued by the substantially higher bias than all the other meta-learners, partially accompanied by the consistency issues. The SW-learner is affected by the same issues as the S-learner in Simulation 1 but manages to substantially reduce the bias for the rest of the simulation designs and is generally close to the performance of the T-learner as seen in the Main Simulation.
In Simulation 2 with balanced treatment and complex nonlinear CATE we also observe, as expected, a very good RMSE performance of the T-learner throughout all sample sizes (see Table (ref) in Appendix (ref)). However, it exhibits quite high variance which is mostly due to the fact \ul{that it estimates} two completely disjoint response functions for estimating the CATE. Furthermore, in this design the R-learner in its full-sample version performs particularly well, which comes rather as a surprise as it pools the two disjoint response functions within the estimation procedure. Nevertheless, the R-learner achieves even lower bias than the T-learner for large samples, but with rather high variance which is a pattern observed across all simulation designs.
Simulation 3 features a highly unbalanced treatment assignment and a constant CATE for which the X-learner performs best as expected, throughout all sample sizes and irrespective of the estimation procedure (see Table (ref) in Appendix (ref)). Indeed, the differences between the particular versions, i.e. full-sample, sample-splitting and cross-fitting, are quite small which is in contrast to the R- and DR-learner confirming the insights from the Main Simulation. Within this highly unbalanced design the estimation of the propensity score function plays a key role as in this case the estimated propensity scores can get \ul{quite} often very close to zero. This, however, does not affect the performance of the X-learner as it uses the propensity scores in a fundamentally different way and even the most extreme \(\{0, 1\}\) scores would be indeed admissible as pointed out by Kunzel2019.\footnote{This would, however, violate the identification assumption of the common support.} On the contrary, the results show that such extreme propensity scores \ul{make now both the R-learner and the DR-learner unstable, with the instability being particularly pronounced in the latter meta-learner.} In the case of the DR-learner the heavy tail problem of the CATE distribution \ul{is aggravated} by more unbalanced treatment assignment as can be seen based on the Jarque-Bera statistic and also on the higher variance of the estimator. While the R-learner manages to avoid this issue by downweighting the observations with extreme propensity scores in less unbalanced settings, \ul{it is not fully able to do so} when the imbalance is very high and there is potentially a large proportion of propensity scores close to 1. This translates into the higher values of the Jarque-Bera statistic as well as to higher variance and higher bias, too. These issues lead ultimately to bad performance in terms of the average RMSE for both the R- and DR-learner.
In Simulation 4 the imbalance in the treatment assignment is less pronounced which should partly reduce the propensity score issues for the R- and DR-learner. Within this simulation design we observe similar patterns as for the Main Simulation. For the small and medium sized samples the X-learner in the full-sample version performs best in terms of the average RMSE, while it gets outperformed by the DR-learner in its cross-fitting version in the largest sample-size. While the R-learner's performance is quite competitive in smaller samples, it lags behind in larger samples as observed in other simulation designs. As a general pattern, the X-learner remains quite stable with respect to the estimation procedure whereas the DR-learner in its sample-splitting and cross-fitting version exhibits substantially faster convergence than the competing estimators. Nonetheless, based on the Jarque-Bera statistic, the heavy tail issue is less pronounced but still present as can be seen in Table (ref) in Appendix (ref).
Lastly, in Simulation 5 the CATE function gets more involved, while the treatment assignment remains unchanged. The results once again resemble the general pattern (for details see Table (ref) in Appendix (ref)). As such the R-learner is competitive mainly in the smaller sample sizes, in this case best performing in the cross-fitting version. The DR-learner in the sample-splitting and cross-fitting version exhibits faster convergence rates, however, in this case the considered sample sizes are not large enough to outperform the X-learner. The X-learner exhibits again little differences regarding the estimation procedure and outperforms the other meta-learners in all performance measures across all sample-sizes.
In order to compare the performance of the meta-learners outside a completely synthetic design of the above simulations we apply the estimators in an arguably more realistic setting using an augmented real dataset. For this purpose we use the dataset from the data challenge of the 2018 Atlantic Causal Inference Conference (2018 ACIC henceforth). This dataset is particularly suitable for a comparison of the meta-learners for two reasons. First, the dataset is based on a randomized control trial in education, namely the National Study of Learning Mindsets (NSLM) by Yeager2019, and thus provides us with a real data example. Second, the dataset has been augmented to an observational setting with measured confounding and known treatment effects Carvalho2019 which enables us to evaluate the performance of the meta-learners for the estimation of CATEs.
The dataset includes a total of \(10'391\) observations with 10 covariates, a simulated continuous outcome and a binary treatment, while the share of treated is approximately \(25\%\).\footnote{The dataset can be retrieved online from \href{https://github.com/grf-labs/grf/blob/master/experiments/acic18/synthetic_data.csv}{GitHub}. We neglect here the information about the additional school ID for simplicity and comparability reasons.} The variables are described in Table (ref) in Appendix (ref). Additionally, to create a sparse large-dimensional setting, similar to the synthetic simulations, we augment the dataset further with \(p=90\) uniformly distributed covariates, i.e. \(X_{11,...,100} \sim \mathcal{U}\big([0, 1]^{n\times p}\big)\) with the same correlation structure as used within the synthetic simulations.\footnote{For more detailed descriptive statistics of the augmented real dataset including correlation heat map of the covariates see Appendix (ref).} At a high level, we are interested in estimating the treatment effects of an intervention to foster a belief to develop intelligence in students on a measure of student achievement, conditional on observed covariates. The CATEs were generated according to the following specification: \[\tau(x)=0.228 + 0.05 \cdot\mathds{1}(x_1 < 0.07) - 0.05 \cdot\mathds{1}(x_2 < - 0.69) - 0.08 \cdot\mathds{1}(c_1 \in {1, 13, 14})\] and while the conditional independence assumption holds by construction, the confounding has a complicated functional form. For a detailed description of the data generating process used for the augmentation see Carvalho2019.
Similarly as in the \ul{synthetic} simulations we estimate the CATEs with all meta-learners and evaluate their performance \ul{with regard to the point estimates}. For this purpose we perform an empirically grounded semi-synthetic simulation study inspired, among others, by Lechner2019, Kunzel2019 and Naghi2021 where we first, set apart a validation set of size \(N=1'000\) observations, and second, sample \(R=\{2'000, 1'000, 500\}\) training sets each of sizes \(N=\{500, 2'000, 8'000\}\) observations from the remaining data. We omit the biggest sample of \(N=32'000\) observations due to the size restrictions of the dataset. We report mean performance measures in a similar fashion as in the previous simulation experiments.
The CATE results of the Semi-synthetic Simulation for all meta-learners are summarized in Table (ref), while Figure (ref) provides details on the meta-learners in the full-sample, double sample-splitting and double cross-fitting versions.
The results reveal a similar picture to the synthetic simulations in general, with the largest similarities to Simulation 3 and 5 in particular. Accordingly, the X-learner achieves the best performance in terms of the average RMSE as well as average absolute bias in all considered sample-sizes, regardless of the estimation procedure. This emphasizes the good performance of the X-learner in settings with unbalanced treatment assignment and sparse CATE function with structural properties. In the largest sample size of $8'000$ observations, also the DR- and R-learner come close to the performance of the X-learner in terms of the average RMSE, while the simpler SW- and T-learner are competitive mainly in the smaller sample sizes. We also observe a slightly faster convergence of the sample-splitting and cross-fitting version of the DR-learner as in the synthetic simulations, however, the limited sample size in this case does not allow for a sufficiently large improvement to outperform the X-learner. Given the smaller sample sizes in the Semi-synthetic Simulation, we are not able to detect the bias-variance trade-off and the sample-splitting versions always exhibit higher values of the average RMSE, average absolute bias and average standard deviation. This is particularly noticeable for the smallest sample size of $500$ observations as there is \ul{essentially} not enough data left after splitting to learn the correct CATE function. For all meta-learners the cross-fitting versions then always perform better in terms of the variance reduction and even lead to a lower bias in comparison to the sample-splitting versions. These results \ul{accentuate} the fact that the benefits of sample-splitting in removing the overfitting bias become apparent only for sufficiently large samples. Additionally, we see larger discrepancies between the estimation versions of the DR- and R-learner in comparison to very stable performance of the X-learner, \ul{similarly as in the synthetic simulations.} Lastly, the results on the distribution of the predicted CATEs resemble those of the synthetic simulations with the heavy tail problem of the DR-learner in its sample-splitting version as well as deviations from normality of the S- and SW-learner.
Given the results of our synthetic and semi-synthetic simulations there are several important findings for the estimation of heterogeneous causal effects by the meta-learners and the associated usage of sample-splitting and cross-fitting which are relevant for applied empirical work.
In general, the results suggest that meta-learners that directly model both the outcome equations and the selection process perform better, especially in larger samples, which is in line with the insights from the previous literature Knaus2021. Meta-learners modelling only the outcome equations are competitive only in smaller samples and tend to perform poorly in larger samples as they fail to properly account for the selection into treatment.
In particular, we do not recommend the usage of the S-learner for estimation of heterogeneous causal effects due to empirically documented undesirable statistical properties such as high bias and consistency issues. \ul{These results are in line with the simulation findings of} Lu2018, \ul{where the S-learner exhibits the highest bias across all simulation settings.} The herein studied modification of the S-learner, the SW-learner, alleviates the high bias of the S-learner, however, it does not solve the consistency issues. Hence, enforcing the treatment variable into the splitting set of the forest does not constitute an attractive option for estimation of causal effects. In contrast, the T-learner does not suffer from high bias or any consistency issues and has a stable performance as it uses the full data sample without the need of sample-splitting due to potential overfitting. Hence, the T-learner might be an interesting option, if a large sample is not available for the empirical analysis. Alaa2018 \ul{similarly find the T-learner to uniformly outperform the S-learner in estimating heterogeneous treatment effects.} Related simulation studies \parencites{Jacob2020}{Knaus2021} find also a relatively competitive performance of the T-learner, especially with the Random Forest as a base learner.
Among the meta-learners based on the estimation of nuisance functions, the X-learner performs very well not only in settings with highly unbalanced but also in less unbalanced treatment shares with simple CATEs and demonstrates the theoretically argued capability to learn such CATE structures Kunzel2019. Moreover, the X-learner exhibits a quite stable performance across all simulation designs with low bias and very low variance, even in small samples. Additionally, due to its particular usage of the propensity scores, the X-learner is not too sensitive to the choice of the estimation procedure. As such, both the full-sample version and the cross-fitting version of the X-learner are viable options, regardless of the sample size. For these reasons, we recommend to use the X-learner for CATE estimation if the researcher is facing a situation with very low number of treated units as well as in less unbalanced settings with potentially limited sample size. In contrast to the X-learner, the DR-learner performs particularly well in settings with nonlinear and complex CATEs if large enough samples are available. However, it tends to be unstable in small samples with unbalanced treatment assignment due to extreme propensity scores, which relates to the results of Jacob2020 and Knaus2021. Additionally, for the DR-learner the choice of the estimation procedure is crucial as its sample-splitting and cross-fitting version exhibits the fastest convergence rates of all meta-learners which highlights the theoretical arguments provided in Kennedy2020. According to the simulation evidence, we advice to employ the cross-fitting version of the DR-learner for the CATE estimation in settings with rather balanced treatment assignment and when large sample is available. Recently, Knaus2020b proposed the normalized DR-learner, that \ul{addresses} the problem of unstable CATE predictions due to extreme propensity scores which might be a viable option for smaller sample sizes and settings with unbalanced treatment shares. Lastly, the simulation evidence suggests that the R-learner is in comparison to the DR-learner less prone to unstable performance due to extreme propensity scores. However, its performance is competitive only in smaller samples, while the empirically approximated speed of convergence is slower than the one of the DR-learner and seems to depend on the CATE complexity as theoretically argued by Nie2021a. With respect to the estimation procedure we do not find a clear-cut evidence in favour of a particular version as both the full-sample as well as the cross-fitting version exhibit comparably good performance. Based on this evidence, the R-learner might be an attractive option for estimation of CATEs if the treatment assignment is not too unbalanced and if only a small sample is available. For comparable sample sizes, Knaus2021 also find the R-learner to have good performance in a variety of settings.
Overall, we point out that based on the simulation evidence, for all meta-learners the approximate convergence rates appear to be substantially slower than the parametric rate of $\sqrt{N}$. This is expected given the insights from previous literature that the estimation of more granular heterogeneous effects is a more difficult task in comparison to the estimation of average effects \parencites(compare e.g.)(){Lechner2019}[or][]{Knaus2021}. However, we note that the approximate convergence rates differ considerably among the meta-learners and their specific implementations as documented in our simulation experiments.
Our simulation evidence suggests that using the full sample for estimation of both the nuisance functions as well as the CATE function leads to a remarkably good performance in terms of both bias and variance in finite samples. Recently, Curth2021 also point out that the full-sample estimation seems to work better in practice, especially for small samples. In theory, we would expect lower variance yet higher bias due to overfitting Chernozhukov2018c. The possible reason for this phenomenon might in our case be due to the out-of-bag predictions of the forest that we use throughout the simulation experiments. Even though these predictions are not out-of-sample per se they are not directly based on the observations used for estimation and as such might help to alleviate the overfitting problem when using full sample Athey2019a. In the causal machine learning literature, such out-of-bag predictions are for example also used in the case of the Generalized Random Forest for the residualization Athey2019, similar to the one used in the R-learner. In contrast to the full-sample estimation, using the double sample-splitting for the estimation of the nuisance functions, we effectively use only one third of the available data. Theoretically, we should observe a smaller bias but higher variance of the estimators. However, in almost all cases we observe both higher bias as well as higher variance, particularly for the small sample sizes. Nonetheless, we document the expected bias-variance trade-off for the largest sample sizes. This stems mainly from the fact that using only a third of the smaller samples does not allow a sensible machine learning estimation of the highly non-linear nuisance functions featured in our simulations. However, especially for the DR-learner we do observe faster convergence rates for the sample-splitting version which is compatible with the theoretical convergence arguments \parencites{Newey2018}{Kennedy2020}. Hence, it seems to be the case that in order to benefit from the double sample-splitting the training sample must be of sufficient size, otherwise the full sample estimation achieves a better performance. Lastly, the double cross-fitting for estimation of the nuisance components effectively uses all \ul{the available information from the} data and substantially reduces the variance of the estimators, while keeping the bias low at the same time. \ul{This comes at the price of longer computation time in comparison to the sample-splitting procedure as the estimation is repeated multiple times. Nevertheless, the computation time of the cross-fitting procedure is on average comparable with the full-sample estimation} (see Appendix (ref) for details).
Based on the above simulation evidence, it seems reasonable to always use the full-sample estimation together with out-of-bag predictions (if available) \ul{when} a relatively small sample is available to the applied empirical researcher, whereas to use the double cross-fitting procedure \ul{when} a relatively large data is accessible, regardless of the choice of a meta-learner. On the contrary, the simulations do not provide any evidence for an advantageous usage of the double sample-splitting over the double cross-fitting, \ul{apart from the computational aspects}.
We investigate the finite sample performance of machine learning based meta-learners for the estimation of heterogeneous causal effects with focus on the specific estimation procedures related to data usage. In particular, we examine the properties of double sample-splitting and double cross-fitting as defined by Newey2018 in contrast to using full sample for estimation. For this purpose, we review several meta-learning algorithms for estimation of causal effects and discuss their advantages and disadvantages in particular empirical settings. We conduct an extensive simulation study with data generating processes involving highly non-linear functional forms and large-dimensional feature space, while keeping the treatment effect specifications well-structured. Furthermore, we perform a semi-synthetic simulation based on an augmented real dataset to reflect an actual empirical setting. Moreover, we repeat the simulation experiments for increasing sample sizes to empirically study the convergence properties of the meta-learners. Based on our simulation evidence, we provide a guideline for applied empirical researchers and practitioners to better inform the decisions of applying certain method and estimation procedure for their particular research objectives.
The results of our simulation study show that the choice of the estimation procedure can indeed largely impact the performance of the meta-learners in finite samples. On the one hand, we provide simulation evidence for the theoretical arguments of the bias-variance trade-off related to sample-splitting and cross-fitting which, however, become apparent only if sufficiently large samples are used. On the other hand, we document the adverse effects of these procedures in small samples, \ul{when using machine learning}. Therefore, we argue that in empirical studies based on small samples, applied empirical researchers should use the full sample for \ul{machine learning} estimation of both the nuisance functions as well as the treatment effect function as the overfitting bias is in such cases of secondary importance. However, for empirical analyses with access to large data samples, we advocate for the usage of the double cross-fitting for the estimation of treatment effects as the overfitting bias here becomes of primary importance. The double cross-fitting procedure then effectively reduces this overfitting bias and \ul{successfully preserves the full sample size efficiency of the estimator.} Moreover, \ul{if computation time is not a constraint}, we discourage applied empirical researchers to use the double sample-splitting procedure due to substantial \ul{increase in variance}, while having no benefit over the double cross-fitting in terms of bias reduction.
In contrast to the typical drawbacks of simulation studies, the particular design of our simulation experiments with varying sample size and varying treatment shares allows us to draw relevant conclusions that are not solely dependent on the particular specification of the data generating processes, but rely on the data characteristics that an applied empirical researcher can observe without imposing arbitrary assumptions. In particular, the simulation evidence implies a clear advantage for the X-learner, when the researcher is confronted with highly unbalanced treatment shares. This finding holds irrespective of the sample size at hand and as such we recommend the usage of X-learner for estimation of heterogeneous treatment effects whenever the share of treated or controls is around \(15\%\) or less. With less unbalanced treatment shares at around \(25\%\) of treated or controls, the size of the available sample becomes decisive. For smaller samples with only few hundred observations ($500$ and $2'000$), the simulation evidence again favours the usage of the X-learner. However, for bigger samples with several thousand observations ($8'000$ and $32'000$), our findings favour the DR-learner as it can successfully learn \ul{highly} complex treatment effect function if enough data is available. Finally, with perfectly balanced treatment shares, the sample size matters less. In such cases, the DR-learner as well as the R-learner are both the preferred estimators. However, we advise against the usage of these two methods in highly unbalanced settings as their performance becomes unstable due to extreme propensity scores. Finally, concerning the simpler meta-learners, we explicitly argue against the usage of the S-learner by applied empirical researchers for estimation of heterogeneous treatment effects due to the herein documented undesirable statistical properties, while the T-learner might be a reasonable choice in small samples with balanced treatment shares.
Even though we shed light on certain finite sample issues of applying different estimation procedures when using meta-learners for estimation of heterogeneous causal effects, our findings raise new relevant questions. Most importantly, the question of conducting statistical inference about the estimated heterogeneous treatment effects is worth further investigations. Based on the insights in this paper it would be of interest to investigate the performance of the bootstrapping for estimation of standard errors as studied by Kunzel2019 for meta-learners based on the double sample-splitting and double cross-fitting procedures. \ul{Moreover, a comparison of such bootstrapping inference procedure for meta-learners and the approaches used in the Causal Forest literature such as the bootstrap of little bags in the Generalized Random Forest} Athey2019 \ul{or the weight-based inference as in the Modified Causal Forest} Lechner2019 \ul{would be desirable.} Furthermore, the performance difference in \ul{the point estimation} using the out-of-bag vs. in-sample predictions could provide additional insights on the benefits of sample-splitting and cross-fitting procedures and hence to assess the robustness of our results to different types of base learners. Finally, a further simulation comparison between the X-learner, the DR-learner and its normalized version as proposed by Knaus2020b for highly unbalanced settings would be of interest.
\printbibliography