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.
113,129 characters · 26 sections · 114 citation commands
Comprehensive Causal Machine Learning
\pagenumbering{arabic}
Machine learning (ML) has paved its way into academia and industry, impacting numerous fields from health care and finance to social sciences. At the core of the ML revolution is the remarkable predictive power of the methods. However, the growing debate around ML emphasizes that prediction does not imply causation. Going beyond mere predictive associations to identify cause-and-effect relationships is at the centre of most questions concerning the effects of policies, medical treatments, marketing campaigns, business decisions, etc. athey2017beyond. The different focus of causal modelling calls for ML approaches that can reliably estimate causal effects, establishing causal machine learning (CML).
CML integrates principles of ML and causal inference. The causality literature provides conditions for identification and estimation of effects. Building on these conditions, causal problems can be transformed into specific prediction problems imbens2009recent for which ML methods are suitable hastie2009. The flexibility of ML combined with careful use of data leads to versatile causal effect estimators that are also able to uncover effect heterogeneity athey2017state.
In recent years, researchers from different disciplines contributed to method development in CML. Although such multidisciplinary work on a common subject is an advantage from a scientific point of view, it raises the question which of these many proposed methods to use in any specific empirical study. In this paper, we focus on the relative performance of a specific class of estimators for a special empirical case, which is popular in the CML literature, namely (i) when causal effects are plausibly identified under unconfoundedness (i.e., selection-on-observables) in a multiple treatment setting, and (ii) the researcher is interested in aggregated average causal effects as well as their possibly fine-grained heterogeneity.
A key condition of any causal effect estimator to be attractive in applications is that it has acceptable proven statistical properties. Such theoretical guarantees and the implied ability to conduct inference are crucial, because, different to prediction settings, realisations from the true effects, the so-called `ground truth', are unobservable. Therefore, they cannot be used to evaluate the performance of the specific estimation. Another important factor in practice is to keep the estimation of the many different causal parameters sufficiently simple, for example, by using the same, or the same type of causal machine learner for all of them. This avoids the time-intensive task of tuning and monitoring many different estimators. Additionally, methods that can predict heterogeneous treatment effects on so far unseen data containing covariates only are also an advantage in empirical settings. We will call methods that (i) have available statistical inference, (ii) can estimate average and fine-grained effects using one type of ML framework, and (iii) can predict using covariates only, as Comprehensive Causal Machine Learners (CCMLs). Finally, it is advantageous to have internal consistency of the possibly many estimated effects, in the sense that effects at the higher aggregation levels are close to appropriately aggregated lower-level effects.\footnote{A similar concept is known as coherence in hierarchical time series forecasting, where it led to development of methods that ensure coherent forecasts across collections of time series formed through aggregation. athanasopoulos2024forecast provide a review of such hierarchical time series methods.} Internal consistency avoids scenarios in which conclusions across different aggregation levels might be contradictory.
\defcitealias{chernozhukov2018dml}{Chernozhukov, Chetverikov, Demirer, Duflo, Hansen, Newey, and Robins, 2018}
In the light of these arguably desirable properties, we analyse three estimation principles that belong to the class of Comprehensive Causal Machine Learning. The methods are Debiased/Double Machine Learning (dml; \citetalias{chernozhukov2018dml}), the Generalized Random Forest (grf; athey2019grf), and the Modified Causal Forest (mcf; lechner2018modified). Only the latter fulfils the condition of internal consistency by construction.\footnote{Within the classification of methods yielding coherent predictions for hierarchical time series, mcf would be classified as a single-level, bottom-up approach.}
The intended contribution of this paper to the literature is two-fold: The first contribution is a large scale-simulation exercise to evaluate the finite-sample performance of all three CCMLs in many different settings to better understand which estimator may have advantages or disadvantages in certain situations and draw practical conclusions. The simulation assesses the CCMLs at three different levels of treatment effect granularity: (i) individual, (ii) group, and (iii) population level. Crucially, these simulations also cover the inference procedures that usually were absent from the small number of existing and far more limited simulation studies. The second contribution are the theoretical guarantees for the mcf and a comparison of the theoretical guarantees across the three different CCMLs. Additionally, similarities and differences between the three approaches are documented.
The results yield several practical recommendations: First, dml-based methods predominantly excel in estimation of average treatment effects (ATEs) or group treatment effects (GATEs) with fewer groups. Second, for finer causal heterogeneity, explicitly outcome-centred forest-based approaches are superior. Third, the mcf offers three additional benefits: (i) it is the most robust estimator even for the ATE in cases when dml-based approaches underperform because of substantial selection into treatment that needs to be corrected for; (ii) it outperforms both dml and grf in GATE estimation when the number of groups gets larger; (iii) it is the only estimator that is internally consistent.
The structure of the paper is as follows: (ref) provides a review with a particular focus on the broad contribution of CML methods to the estimation of the ATE and more fine-grained conditional average treatment effects (CATEs) under unconfoundedness. After a general overview of CML methods, the section focuses on methods able to estimate effects at all levels of granularity and identifies the CCML methods. The section ends with an overview of CML simulation comparisons and applications. (ref) introduces the notation of the potential outcome model, defines the parameters of interest, and states the identifying assumptions imposed under unconfoundedness. (ref) discusses the three CCMLs. The large-scale simulation study is summarized in (ref). (ref) concludes. (ref) provides proven theoretical guarantees for the mcf, its inference procedure, a comparison to the grf, as well as some new results for the grf. Implementational details of the simulation study are collected in (ref), while (ref) holds its detailed results.
For causal inference under the unconfoundedness assumption rubin1974estimating, the ATE represents one of the main parameters of interest capturing the overall impact of a treatment or an intervention. For ATE estimators, such as inverse probability weighting (IPW) and propensity score matching imbens2004nonparametric, propensity score estimation plays a central role. Both estimators are consistent when a parametric model for the propensity score is correctly specified. Consequently, one of the first applications of ML methods in ATE estimation involved flexible propensity score estimation. Various ML approaches, such as Neural Networks, Support Vector Machines, Decision Trees westreich2010propensity or Boosting Models mccaffrey2004propensity,westreich2010propensity, have been proposed to improve the estimation of the propensity score.
However, even flexible ML methods may exhibit poor performance if they prioritize fitting the propensity score well instead of explicitly balancing the covariates across treatment groups, which is important to reduce the bias of the estimator for the causal parameter of interest. An estimator of the propensity score that explicitly aims to maximize covariate balancing is introduced in imai2014covariate. Alternatively, graham2012inverse improve IPW estimation by covariate balancing within the empirical likelihood framework, without explicit modelling of the propensity score.
cannas2019comparison investigate the performance of matching and weighting estimators based on propensity scores obtained from Logistic Regression, Decision Tree, Bagging, Boosting, Random Forest, Neural Networks, and Na\"{i}ve Bayes. Their simulation results indicate that Random Forests consistently outperform other methods, particularly in the context of IPW. goller2020does compare Random Forests, LASSO Logit Regression, and Probit estimation of the propensity score for a matching estimator in an active labour market policy setting, revealing that LASSO may yield more credible results than conventional propensity score estimation methods, but overall results were mixed.
Under the standard unconfoundedness assumptions, to be detailed in (ref), ATE is alternatively identified as the expected difference between conditional expectations of the outcome variable for the treated and non-treated populations. This leads to regression-based estimation of the two conditional expectations. Assuming a correct specification of the outcome models, averaging their estimated differences yields a consistent ATE estimator. In the ML literature, tree-based methods such as Bayesian Additive Regression Trees (BART) hill2011bayesian and Ensemble Methods austin2012using were introduced to obtain regression-based estimators of ATEs. However, issues like regularization bias athey2018approximate and slow convergence of ML methods render the outcome-based approach less popular.
Methods that combine propensity score-based and regression-based estimators to increase robustness to misspecification have gained popularity. Two related approaches prevail in the literature: Double-robust (DR) estimators and Neyman-orthogonal scores. DR estimators utilize two nuisance functions, i.e., the propensity score and outcome models. They are consistent as long as at least one of the nuisance functions is correctly specified. Certain DR estimators were shown to be semi-parametrically efficient if both nuisance functions are correctly specified. Ideas to combine propensity score and outcome modelling in a double-robust way originally emerged in robins1994estimation, robins1995analysis, leading to the Augmented Inverse Probability Weighting (AIPW) estimator for the ATE. Nevertheless, sensitivity to extreme propensity scores, as in the case of the IPW estimator, remains an issue.\footnote{For potential benefits of trimming extreme propensity scores, see e.g., huber2013performance.}
Approaches based on Neyman-orthogonality ensure that the moment conditions identifying the target parameters are not locally affected by small perturbances in the nuisance parameter estimates chernozhukov2018dml. Such an orthogonality property is leveraged, e.g., in belloni2014inference for ATE estimation with binary treatment and in farrell2015robust for multiple treatments. As Neyman-orthogonal scores mitigate first order bias arising from ML estimation of nuisance parameters bach2024dmlr, Neyman-orthogonality emerged as a key element for the Double Machine Learning (dml) framework introduced in chernozhukov2018dml. dml is a generic framework for obtaining consistent estimators and valid inference for low-dimensional parameters. This framework combines Neyman-orthogonal scores with cross-fitting, a data partitioning technique that helps avoid further biases arising from the ML estimation of nuisance parameters. dml is well-suited for high-dimensional settings, as it operates under convergence rate conditions for nuisance parameter estimation that are achievable by many ML methods. chernozhukov2018dml apply the dml framework to ATE estimation using the DR score of robins1995semiparametric\footnote{Note that the DR score of robins1995semiparametric is Neyman-orthogonal.} and derive its statistical properties.
Another approach for ATE estimation is Targeted Minimum Loss Estimation (tmle), a general method that combines ML methods with a targeted updating step to reduce bias and variance of the estimated parameter of interest in semi-parametric and non-parametric models van2006targeted. Among various tmle adaptations\footnote{See Laan2018 for a comprehensive overview.}, cross-validated tmle, cv-tmle, proposed by Zheng2011, shares considerable similarities with dml. Both approaches require estimation of nuisance parameters that are combined in a way that yields a semi-parametrically efficient estimator of the ATE that is robust to small errors in nuisance parameter estimation under the same convergence rate conditions for nuisance parameter estimation and use of data splitting. Meanwhile dml utilizes the nuisance parameters in a Neyman-orthogonal score to mitigate the bias, \emph{cv-tmle} combines them in an iterative estimation procedure addressing the bias and yielding quantities for each observation that can be averaged into an ATE. Despite this difference, \emph{cv-tmle} and \emph{dml} are asymptotically equivalent and can be expected to yield similar results in large samples under the same choice of methods estimating the nuisance parameters, as noted by knaus2021double.
Automatic Debiased Machine Learning, introduced in chernozhukov2022automatic, offers another alternative for ATE estimation that can deal with bias stemming from using ML methods for estimation of nuisance parameters. The key idea is based on a “de-biased moment” utilizing an existence of a Riesz representer for parameters such as the ATE. The debiasing is automatic in the sense that it does not require any knowledge about the form of the Riesz representer. The Riesz representer is directly estimated within the procedure. In case of the ATE, the de-biased moment function corresponds to the DR moment function of robins1995semiparametric. The estimation of the Riesz representer and nuisance parameters requires cross-fitting as in dml. The method is shown to yield a consistent and asymptotically normal estimator and provides a consistent asymptotic variance estimator.
\defcitealias{wei2023efficient}{Wei, Petersen, van der Laan, Zheng, Wu, and Wang (2023)}
CML methods extend beyond ATE estimation and help uncover heterogeneity of the treatment effects through analysis of Conditional Average Treatment Effects (CATEs). CATEs quantify how the treatment affects an individual unit with specific low or large-dimensional characteristics. Prominent methods estimating CATEs include Causal Forests wager2018estimation,athey2019grf,lechner2018modified, along with dml-based estimators yielding estimated components of the DR efficient score that are further regressed on the covariates as first proposed by Laan2006. Alternative approaches include Meta-learners, which use multiple machine learning algorithms as “base learners” to estimate large-dimensional CATEs. The X-learner, introduced by kunzel2019metalearners, is effective in cases with uneven treatment groups, while the R-learner, developed by nie2021quasi, leverages the robinson1988root transformation to estimate heterogeneous causal effects. When interested in CATEs for specific subgroups, the above-mentioned quantities of the tmle iterative procedure can be averaged into a CATE, using the corresponding part of the sample for the estimation.\footnote{\citetalias{wei2023efficient} note that the natural implementation of tmle “one-group-at-a-time” for CATE estimation with discrete groups might be computationally costly. As a remedy, they propose a tmle method that targets multi-subgroup treatment effects simultaneously.} In case of a CATE for the finest granularity level, the \emph{tmle} estimand depends on an observed outcome and cannot serve for a prediction on a new sample containing covariates only.\footnote{Based on levy2021fundamental who introduces a \emph{tmle} method for simultaneous estimation of ATE and variance of treatment effects, where estimated CATEs serve as nuisance parameters.}
The influential paper by chernozhukov2018dml, averaging the estimated components of individual DR efficient scores of robins1995semiparametric to estimate ATE in a dml framework, spanned further methodological contributions, particularly in discovering heterogeneity along a chosen set of covariates, such as gender or age group. dml-based estimation of CATEs exploits the fact that the conditional expectation of the DR efficient score given such covariates identifies the respective CATEs. For low-dimensional sets of chosen covariates, OLS, Series, or Kernel Regressions of the estimated DR score components on the chosen set of covariates estimate the low-dimensional CATE and standard statistical inference applies, as shown in semenova2021debiased, zimmert2019nonparametric, and fan2022estimation. For a higher dimensionality of the set of the chosen covariates, kennedy2023towards introduces the DR-learner that takes the components of a classic DR estimator of ATE and regresses them on all covariates to estimate a large-dimensional CATE and derives an upper bound on the DR-learner error relative to an oracle.
The foundations for heterogeneous treatment effect estimation by Causal Forests are laid in a seminal paper by athey2016recursive. This work introduces splitting rules for treatment effect estimation and honest data splitting, i.e. using different parts of data to build the tree and to estimate the parameters, to ensure valid inference for CATEs through tree-based methods. Subsequently, wager2018estimation develop a nonparametric Causal Forest algorithm based on honest trees and splitting rules maximizing heterogeneity in treatment effects across final leaves. The last step of the estimator of the possibly large-dimensional CATE is obtained by averaging individual tree estimates. Concerning its asymptotic guarantees, the Causal Forest is consistent and pointwise asymptotically normal under further regularity conditions. In a distinct but related approach, athey2019grf introduce the Generalized Random Forests (grf). This method is based on an honest forest with a gradient-based approximation of an estimator-specific splitting rule. The forest estimation determines importance weights of each observation for the parameters of interest. In other words, these weights determine a neighbourhood of observations that will contribute to the CATE estimation. Instead of estimating the CATE directly as an average across honest trees as in wager2018estimation, the forest weights enter a set of local moment conditions identifying the CATE. This estimator is consistent and pointwise asymptotically Gaussian and generalizes beyond unconfoundedness.
Additionally, lechner2018modified introduces a Modified Causal Forest (mcf) estimation procedure for ATEs and CATEs. In a multiple treatment framework, he proposes a weights-based inference procedure utilizing the weighted-outcome representation of forest estimator. mcf further differs from the other Causal Forests in using the mean squared error of the CATE directly to find the best split. Furthermore, it introduces a two-sample honesty, detailed in (ref), for building the forest and the estimation of the effects. A simulation study in lechner2018modified shows superior performance of mcf to the Causal Forest of wager2018estimation. bodory2022high show that the mcf works in applications replicating results from several papers in different fields.
Methods that (i) can estimate average and fine-grained effects using one type of ML framework, (ii) have available statistical inference for all the effects, and (iii) can predict using covariates only, are labelled as Comprehensive Causal Machine Learners (CCMLs). Some of the above-mentioned advances provide practitioners with methods able to estimate treatment effects across all different levels of granularity. Table (ref) provides an overview of commonly used CML approaches for estimation of such treatment effects and evaluates them based on (i) available statistical inference, (ii) their ability to predict large-dimensional CATEs on new data (with covariate information only), and (iii) internal consistency\footnote{An estimator is considered being internally consistent when effects at the higher aggregation levels are close to appropriately aggregated lower-level effects.} of estimators. The overview reveals that tmle and Meta-learners do not have available statistical inference for all levels and that tmle cannot predict CATEs at the finest granularity level on a new data set that contains covariates only. Thus, these two methods are not considered to be CCMLs. On the other hand, dml and Causal Forests are approaches that fulfil the conditions for comprehensive treatment effect evaluation with available statistical inference and predictions for new data containing covariates only. For Causal Forests, there is a Bayesian alternative which will not be pursued further because its inference procedure is not based on repeated sampling inference.
Focusing on estimation at different levels of granularity via CCMLs, dml-based estimation operates via the estimated components of the DR efficient scores that are regressed on covariates to yield lower and higher-level aggregates of average treatment effects. Causal Forests, grf and mcf, provide in their first stage CATE estimates at the highest level of granularity. Using direct aggregation step of these CATEs, mcf leverages the weighted-outcome representation of forest predictions, while grf uses a variant of the AIPW score to yield higher level aggregations via a regression step. Table (ref) summarizes the estimation at different levels of granularity for CCMLs and other methods.
\defcitealias{wendling2018comparing}{Wendling, Jung, Callahan, Schuler, Shah, and Gallego (2018)} \defcitealias{audrino2024does}{Audrino, Chassot, Huang, Knaus, Lechner, and Ortega, 2024} \defcitealias{langenberger2023exploring}{Langenberger, Steinbeck, Sch\"{o}ner, Busse, Pross, and Kuklinski, 2023}
The literature on large scale simulation comparisons of finite sample properties of CML methods is very limited. Among notable contributions, knaus2021machine look at the finite-sample performance of selected machine learning estimators, covering grf and generic approaches that can be combined with any ML method, yielding large-dimensional CATE estimates based on pseudo-outcomes or modified covariates.\footnote{At the time the simulations of knaus2021machine were performed, the mcf was not yet available.} Higher level aggregates are simple averages of the corresponding CATEs. Using the Empirical Monte Carlo Study (EMCS) approach, many components of their simulations are based on real data. The results show that best-performing estimators at the large-dimensional CATE level produced most reliable estimates at higher aggregation levels, such as for the ATE, in terms of their mean squared error (MSE). In general, methods that utilized both outcome and treatment equations are among the best-performing methods, including grf with local centering. caron2022estimating evaluate several Meta-learners in an EMCS based on health data. Their results highlight variability in performance that depends on the complexity of the data generating process. In a setting with a complex CATE, multitask learners designed to estimate both outcome equations jointly, such as Multitask Gaussian Process introduced in alaa2017bayesian, perform best in terms of MSE. The X-learner performs the best in a setting with a simple CATE function and slightly unbalanced treatment groups. Another health-data-based EMCS of \citetalias{wendling2018comparing} concentrates on large-dimensional CATE estimation in a common healthcare setting when outcomes are binary and rare. The compared methods include BART, grf with local centering, Causal Boosting, and Causal Multivariate Adaptive Regression Splines (MARS). The findings show that BART and Causal Boosting perform better across the scenarios, as evaluated by the root MSE in the conditional probability (risk) difference estimates. Except for one scenario involving highly heterogeneous treatment effects, the coverage rate for BART, grf, and causal MARS was close to its nominal level.\footnote{To the best of our knowledge, this is the only study that also investigates the performance of the corresponding inference procedures. The simulation study in wendling2018comparing differs from the simulation in (ref) by using fixed covariates across replications. Despite this difference, their results are in line with the results in (ref) confirming that high treatment heterogeneity negatively affects the coverage probability. Meanwhile DGPs with low treatment heterogeneity have coverage probability closer to nominal rates.}
Empirical applications using CML for treatment estimation span across many fields and topics. The following non-exhaustive list illustrates the variety of applications ranging from active labour market policy davis2017using, bertrand2021workfare, Pytka2021Understanding, knaus2022heterogeneous, cockx2023priority, burlat2024everybody, education knaus2021double, farbmacher2021heterogeneous, social experiments athey2019estimating, strittmatter2023value, energy use knittel2021machine, natural resource rents hodler2023institutions, finance \citepalias{audrino2024does}, and the dating market boller2021effect to medicine \citepalias{langenberger2023exploring}.
Regarding the comprehensive nature of the dml framework, knaus2022double demonstrates its flexibility in an active labour market programme evaluation setting. rehill2024applied documents and evaluates some practices of using Causal Forests in applied work. Further illustrations of potential benefits of CML compared to traditional methods in empirical work is provided in baiardi2024value. They revisit estimates of average and heterogenous treatment effects from published observational studies and randomized control trials. Their results illustrate the potential of CML to capture non-linear confounding, its suitability in scenarios with large number of covariates relative to the sample size and its ability to discover treatment effect heterogeneity.
We use Rubin's potential outcome language rubin1974estimating to describe a multiple treatment model under unconfoundedness, also referred to as selection-on-observables or conditional independence imbens2000role, Lechner2001. Let $D$ denote the treatment that may take a known number of $M$ different integer values from 0 to $M-1$. The (potential) outcome of interest that realises under treatment $d$ is denoted by $Y^d$. For each observation, we observe only the particular potential outcome corresponding to the observed treatment status, i.e., ${y_i} = \sum\limits_{d = 0}^{M - 1} {\underline 1 ({d_i} = d)y_i^d}$, where $\underline 1(\cdot)$ is the indicator function which equals one if its argument is true.\footnote{If not obvious otherwise, capital letters denote random variables, and small letters their values. Small values subscripted by `$i$' denote the value of the respective variable of observation `$i$'.} There are two groups of variables to condition on, $\tilde X$ and $Z$. $\tilde X$ contains those covariates needed to correct for selection bias (confounders), while $Z$ contains variables that define (groups of) population members for which an average causal effect estimate is desired. $\tilde X$ and $Z$ may be discrete, continuous, or both. They may overlap in any way. The union of these two groups of variables is denoted by $X$, $X = (\tilde X,Z)$, $dim(X)=p$.\footnote{In all CCMLs $p$ is assumed to be finite.}
Below, we investigate the following average causal effects:
The Individualized Average Treatment Effects (IATEs) measure the mean impact of treatment $m$ compared to treatment $l$ for units with covariates $x$. The IATEs represent the causal parameters at the finest aggregation level of the covariates. On the other extreme, the Average Treatment \textbf{E}ffects (ATEs) represent the population averages. ATE is a classical parameter investigated in many empirical studies. The \textbf{G}roup \textbf{A}verage \textbf{T}reatment \textbf{E}ffect (GATE) parameters are in-between those two extremes with respect to their aggregation levels.\footnote{We presume that the analyst selects the variables $Z$ prior to estimation. However, the estimated IATEs may be analysed by methods picking $Z$ in a data-driven way to describe their dependence on certain covariates. See Section 6 in lechner2018modified for more details. Note that abrevaya2015estimating and lee2017doubly introduce similar aggregated parameters that depend on a reduced conditioning set and discuss inference in their specific settings.} IATEs and GATEs are special cases of the already mentioned \textbf{C}onditional \textbf{A}verage \textbf{T}reatment \textbf{E}ffects (CATEs).
The following set of assumptions identifies the causal effects discussed in the previous section imbens2000role, Lechner2001 (see Imbens, 2000, Lechner 2001):\footnote{To simplify the notation, we take the strongest form of these assumptions. Some parameters are identified under weaker conditions as well (for details see imbens2000role, imbens2004nonparametric, or Lechner2001.}
The conditional independence assumption (CIA) holds if there are no covariates other than $X$ that jointly influence treatment and potential outcomes (for the values of $X$ that are in the support of interest, $\chi$). The common support (CS) assumption stipulates that for each value in $\chi$, there must be the possibility to observe all treatments. The stable-unit-treatment-value assumption (SUTVA, also referred to as observation rule or consistency condition) implies that the observed value of the treatment and the outcome does not depend on the treatment allocation of the other population members (ruling out spillover and treatment scale effects). Usually, to have an interesting interpretation of the effects, it is required that $X$ is not influenced by the treatment (exogeneity). In addition to these identifying assumptions, assume that a large random sample of size $N$ from the i.i.d. random variables $Y, X, D$, $(y_{i}, x_{i}, d_{i})$, $i=1, \ldots, N$, is available and that all necessary moments of these random variables exist.
If these assumptions hold, all IATEs are identified in the sense that they can be uniquely deduced from variables that have observable sample realisations hurwicz1950:
Since the distributions used for aggregation, ${f_{\tilde X| Z = z}}(\tilde x)$ and ${f_X}(x)$, relate to observable variables ($X$) only, they are identified as well (under standard regularity conditions). This in turn implies that the GATE and ATE parameters are identified.
In this section, three comprehensive approaches for treatment effect estimation are presented. In the first two subsections, the underlying principles as well as the concrete estimation and inference algorithms for Double Machine Learning (dml) and Generalized Random Forest (grf) are reviewed. The third subsection focuses on the Modified Causal Forest (mcf). While the basic principles of the mcf are introduced in lechner2018modified, that paper does not contain explicit theoretical guarantees. Therefore, (ref) provides the proofs of consistency and asymptotic Gaussianity properties for IATEs, GATEs, and ATEs.
dml estimation chernozhukov2018dml is based on moment conditions with scores satisfying the identification condition as well as Neyman-orthogonality, yielding estimators that are robust to small estimation errors in nuisance parameters neyman1959optimal. A Neyman-orthogonal score that identifies ATE, GATE, and IATE for treatments $m$ and $l$ under assumptions outlined in (ref) is the DR score of robins1995semiparametric:
where $\mu_{d}(x)=E(Y| X=x, D=d)$, $\eta_{m,l}^{dml}(X)=(\mu_{m}(X),\mu_{l}(X), p_{m}(X), p_{l}(X))$ captures the nuisance parameters, and $\theta$ represents the treatment effect of interest. $O$ represents all observable variables, i.e., $O = (X, Y, D)$. Noting that the difference of the third and fourth term in this score has expectation zero conditional on $X$ at the true values of nuisance parameters, and letting $\eta^{0}$ and $\theta^{0}$ denote the true values of $\eta$ and $\theta$ respectively, $\psi_{m,l}^{dml}(O;{\theta_{m,l}},\eta_{m,l}^{dml}(X))$ identifies the different treatment effects:
Effectively, $ATE(m,l)$ is estimated as an average of the estimated component of the DR score, $\Gamma_{m,l}^{dml}\left({{o_i};\hat \eta _{m,l}^{dml, -k(i)}({x_i})} \right)$, across all observations, i.e., $$\widehat {AT{E^{dml}}}(m,l) = \hat \theta _{m,l}^{dml} = \frac{1}{N}\sum\limits_{i = 1}^N {\Gamma _{m,l}^{dml}\left( {{o_i};\hat \eta _{m,l}^{dml, - k(i)}({x_i})} \right)}.$$ This process involves the estimation of nuisance parameters, symbolized as $\hat \eta_{m,l}^{dml, - k(i)}({x_i})$, through $K$-fold cross-fitting. The latter ensures that for each observation $i$, the corresponding nuisance parameters are estimated without using that specific observation in their training data. Particularly, the elements of the vector $\hat \eta_{m,l}^{dml, - k(i)}({x_i})$ are estimated on $K-1$ folds not containing the observation $i$, which resides in the left-out fold $k$. This is indicated by the superscript $-k(i)$.
Neyman-orthogonal scores mitigate regularization bias in the estimation of the ATE, provided that the product of the convergence rates of the two nuisance parameter estimators within each treatment group is at least $\sqrt N $. This condition, which many ML methods achieve, allows for flexible estimation of nuisance parameters leading to estimators robust to small estimation errors in these nuisance parameters. The additional combination with cross-fitting is important to avoid overfitting when estimating the nuisance parameters and to subsequently guarantee $\sqrt N - $consistent estimation of the main parameters of interest. For technical details, refer to chernozhukov2018dml. The variance estimator of the dml can be computed as the sample variance of the estimated DR scores.\footnote{See Theorem 3.2 in chernozhukov2018dml, and Theorem 1 and Remark 1 in bach2024dmlr.} The dml estimator asymptotically reaches the semiparametric efficiency bound of hahn1998role.\footnote{See Theorem 5.1 in chernozhukov2018dml.}
Building on the identification result, $GATE(m,l;z)$ can be estimated by regressing $\Gamma_{m,l}^{dml}\left({{o_i};\hat \eta _{m,l}^{dml, -k(i)}({x_i})} \right)$ on a low-dimensional vector $Z$ of pre-specified variables. semenova2021debiased provide statistical inference for the best linear predictor in this case, a method that is implemented in (ref). Standard errors are estimated via heteroscedasticity-robust standard errors for asymptotically valid pointwise and uniform confidence bands. In general, the results in semenova2021debiased apply to linear projections onto a vector of basis functions. Under sufficient smoothness and richness of the basis functions, inference targets the $GATE(m,l;z)$ function. When basis functions are group indicators, it targets the corresponding group parameters. For continuous regressors, kernel regression offers a viable alternative, enabling the estimation of $GATE(m,l;z)$ and facilitating statistical inference, as demonstrated in zimmert2019nonparametric and fan2022estimation.
In a similar fashion, a natural approach to obtain an estimator of $IATE(m,l;x)$ is by regressing $\Gamma_{m,l}^{dml}\left({{o_i};\hat \eta _{m,l}^{dml, -k(i)}({x_i})} \right)$ on the covariates.\footnote{In (ref) in Table (ref), this step is called a smoothing step.} foster2023orthogonal and kennedy2023towards derive error bounds for this two-step procedure, referred to as DR-learner in kennedy2023towards. Particularly, the findings in kennedy2023towards provide a theoretical foundation for the validity of inferential methods used in this context for linear smoothers satisfying the required stability condition and under further small-bias condition on $\Gamma _{m,l}^{dml}\left( {{o_i};\hat \eta _{m,l}^{dml, - k(i)}({x_i})} \right)$. In practical applications, different techniques have been employed: goller2023analysing utilizes linear regression, while knaus2022double opts for a Random Forest approach.
\defcitealias{athey2019grf}{ATW19}
grf (athey2019grf, abbreviated as ATW19 in this section) extends the Random Forest into a nonparametric method estimating parameters identified by local moment conditions: $$E\left(\psi^{grf}(O;{\theta^0}(x),{\eta^{grf,0}}(x))\left| {X = x}\right. \right) = 0 \quad \text{for all } x \in \chi,$$ where ${\psi ^{grf}}(\cdot )$ is a score function identifying the true values of $\theta(x)$. As before, $\eta^{grf}(x)$ captures nuisance parameters, $\theta^{0}(x)$ and $\eta^{grf,0}(x)$ are the true values of $\theta(x)$ and $\eta^{grf}(x)$, and $O$ represents all observable variables. grf estimation of heterogeneous treatment effects is based on $\psi^{grf}(Y,\tilde D,X;\theta(x),\eta_{IATE}^{grf}(x)) = \left({Y - \tilde D'\theta (x) - \eta _{IATE}^{grf}(x)}\right) \left(1 \ \ \tilde D'\right)',$ where $\tilde D$ is an $(M-1) \times 1$ vector containing dummy variables indicating whether treatment $d \in \left\{ {1,...,M - 1} \right\}$ was received, $\theta(x)$ is a parameter vector representing corresponding treatment effects (all treatment $d$ vs 0 pairs) for point $x$, and $\eta^{grf}_{IATE}(x)$ captures the intercept and all confounders in a single nuisance parameter.\footnote{This part is coined as an intercept term $c(x)$ in \citetalias{athey2019grf}.} Note that the score function in dml is based on a fully nonparametric DR score identifying treatment effects for all pairwise combinations of treatments, while the score function in grf stems from a partial linear model in which only treatment effects with respect to a reference treatment are identified. The reference treatment is the one left out from the vector $\tilde D$, here the control group with $d = 0$. This is not too restrictive if effects that are conditioned on the treatment status are not of interest (as for IATE estimation), because the point estimates of the other treatment combinations can be obtained as $\theta_{m,l}(x)=\theta_{m,0}(x)-\theta_{l,0}(x)$.
Building on the local generalized method of moments, grf estimates IATEs, $\theta(x)$, as $$\left( {{{\hat \theta }^{grf}}(x),\hat \eta _{IATE}^{grf}(x)} \right) \in \mathop {\arg \min }\limits_{\theta (x),\eta (x)} \;{\left\| {\sum\limits_{i = 1}^N {w_i^{grf}(x){\psi ^{grf}}({o_i};\theta (x),\eta _{IATE}^{grf}(x))} } \right\|_2},$$ where the weights, $w_i^{grf}(x)$, are obtained from an honest Random Forest whose gradient-based splitting rule is designed to maximize treatment heterogeneity in the daughter leaves. The required Honesty assumption involves a data partitioning strategy to prevent overfitting. The subsample drawn for each tree is further split into two halves: one for building the tree structure and the other one for estimating the parameters of interest -- in this case, the tree weights which are subsequently aggregated into forest weights. As all observations alter between the two halves across trees when subsampled, we call this procedure `one-sample honesty'. The forest weights, summing up to 1, can be seen as a measure of relevance of the observation $i$ for the estimation of $\theta(x)$. Given the forest weights $w_i^{grf}(x)$, the solution to the given optimization problem is:
where ${\bar{\tilde{d}}_w} = \sum\limits_{i = 1}^N {w_i^{grf}(x){{\tilde d}_i}} $ and ${\bar y_w} = \sum\limits_{i = 1}^N {w_i^{grf}(x){y_i}} $.\footnote{Compare with eq. (19) in \citetalias{athey2019grf}.}
As the splits are chosen to maximize treatment heterogeneity across daughter leaves, to mitigate confounding effects at early splits, \citetalias{athey2019grf} recommend to partial out effects of the covariates $X$ on the outcome $Y$ and treatment assignment $D$ as in robinson1988root. This means building the forest using locally centred values $y_i^{cent} = {y_i} - {\hat \mu ^{( - i)}}({x_i})$ and $\tilde d_i^{cent} = {\tilde d_i} - \hat p_d^{( - i)}({x_i})$ for all observations, where $\mu(x)=E(Y| X=x)$ and the superscript $(-i)$ denotes that the prediction did not use observation $i$. The best practice is to estimate $\mu(x)$ and $p_{d}(x)$, e.g., by $K$-fold cross-fitting as mentioned in \citetalias{athey2019grf}. However, the available grf package performs local centring by subtracting out-of-bag predictions\footnote{In \citetalias{athey2019grf} out-of-bag predictions are called leave-one-out predictions as the final prediction for observation $i$ averages predictions that are obtained from trees that did not use observation $i$ for splitting.} obtained from regression forests specifically trained to predict the outcome variable $Y$ and the probabilities of treatment assignment, given the covariates $X$. Despite its computational attractiveness, this implementation violates the requirement of complete separation of training and prediction data consistent with the theoretical results. Also, the simulation results in (ref) reveal the necessity of $K$-fold cross-fitting over out-of-bag predictions for local centring to remove bias.
The asymptotic results for ${\hat \theta ^{grf}}(x)$ are based on their linear approximation motivated by the method of influence functions. \citetalias{athey2019grf} prove that the two are coupled and asymptotic properties of one apply to the other. Since the linear approximation can be interpreted as an honest Causal Forest estimator introduced in wager2018estimation, their asymptotic result can be applied to establish asymptotic Gaussianity of ${\hat \theta ^{grf}}(x)$, given that $\theta(x)$ and $\eta^{grf}_{IATE}(x)$ are consistently estimated.\footnote{This can be achieved, e.g., via honest forests as outlined in \citetalias{athey2019grf} in Theorem 3.} The linear approximation is also leveraged to derive an estimator of the pointwise standard errors for ${\hat \theta ^{grf}}(x)$. Derivation of the variance of the linear approximation yields two terms. One term can be in general consistently estimated by means of regressions and is a constant in case of $IATE(x)$. The other can be seen as an output of a regression forest with weights $w_i^{grf}(x)$ and the score function as outcome variable. \citetalias{athey2019grf} propose to estimate this term by a variant of a so-called bootstrap of little bags. This method was introduced in sexton2009standard. \citetalias{athey2019grf} motivate it by the observation that “building confidence intervals via half-sampling -- whereby evaluating an estimator on random halves of the training data to estimate its sampling error -- is closely related to the bootstrap efron1982”. A computationally efficient implementation to estimate within one forest, both the forest weights and the element of the variance by half-sampling, requires drawing for each small bag of trees a random half of the sample that is available for further subsampling. Honesty additionally requires the trees to be built on half of the subsample, and the weights $w_i^{grf}(x)$ to be computed on the half that was not used to build the tree. Thus, each $w_i^{grf}(x)$ is estimated on approximately a quarter of all the trees in the forest when the subsampling rate is close to 1.\footnote{Figure (ref) in (ref) captures the grf procedure graphically.} Consistent estimation of the variance term yields asymptotically valid confidence intervals for the IATE at point $x$.
Like dml, grf offers comprehensive estimation of treatment effects at all levels of granularity. However, unlike dml, which estimates all effects in two steps by estimating the $\Gamma _{m,l}^{dml}\left( {{o_i};\hat \eta _{m,l}^{dml, - k(i)}({x_i})} \right)$ first, and subsequently regressing them on a constant, group indicators or regressors to obtain ATE, GATE and IATE estimators, grf takes a more direct approach. Specifically, grf estimates IATEs “directly”, while ATE and GATEs are estimated through additional regressions. The grf package implements ATE estimation (treatment $m$ vs no treatment) by a variant of the AIPW estimator, plugging in the estimates from IATE estimation for the nuisance parameters:\footnote{The equations for estimating the ATE are reconstructed from the \emph{grf} R package (version 2.3.1). The derivation involves code implemented in the functions “get_scores.R” and “average_treatment_effect.R”.}
As default, the grf package estimates the nuisance parameters $\mu(x), p_{1}(x), \ldots, p_{M-1}(x)$ by out-of-bag predictions from Random Forests, and $\theta_{1,0}(x), \ldots, \theta_{M-1,0}(x)$ by a weighted local moment condition with forest-based weights. Derivations in (ref) show the equivalence of the grf score $\psi _{m,0}^{grf}(O;{\theta _{m,0}},\eta _{ATE(m,0)}^{grf}(X)) = \Gamma _{m,0}^{grf}(O;\eta _{ATE(m,0)}^{grf}(X)) - {\theta _{m,0}}$ to the DR score of robins1995semiparametric for the ATE. It can therefore be expected that the grf estimator of the ATE is asymptotically efficient when all nuisance parameters are consistently estimated. Due to a different set of nuisance parameters, consistent estimation of the ATE is robust only to misspecifications of $\mu(x)$ and $\theta_{1,0}(x), \ldots, \theta_{M-1,0}(x)$ (see (ref)). Consistent estimation of the ATE based on the above grf score therefore requires consistent estimation of the propensity scores. (ref) also shows Neyman-orthogonality of the grf score for the ATE, guaranteeing that the moment function is locally insensitive to errors in the nuisance parameter estimates. To obtain the GATE estimators, $\Gamma _{m,0}^{grf}({o_i};\hat \eta _{ATE}^{grf}({x_i}))$ are regressed on $Z$, as in dml. Inference for the ATE and GATEs is conducted as in \emph{dml}.
\defcitealias{wager2018estimation}{WA18} \defcitealias{lechner2018modified}{L18}
lechner2018modified builds on wager2018estimation and introduces the Modified Causal Forest estimator, mcf. One of the procedures proposed by \citetalias{wager2018estimation} builds trees by taking the outcome $Y$ as dependent variable and finding splits that maximize treatment effect heterogeneity. The tree estimator of an IATE at point $x$ is obtained from the final leaf containing $x$ by differencing the average outcomes of the treatment groups evaluated on the subsample that was not used to build the tree. The final forest estimator of the IATE at point $x$ is the average of all tree estimates. This procedure works well in an experimental design with heterogeneous treatment effects but may do poorly in the presence of confounding (see \citetalias{wager2018estimation}). This observation led to further extensions addressing the confounding issue.
grf introduces local centring inside its Random Forest to remove the confounding effects by transforming the data when searching for splits that maximize effect heterogeneity. mcf takes a different approach. Motivated by the fact that in the presence of selection bias the difference of the outcome means of treated and controls within all leaves will not correspond to the means of the true effects, mcf proposes a splitting criterion targeting the minimisation of the MSE of IATE directly. In addition, a penalty is added penalizing squared differences between propensity scores. Taken together, the chosen split is predictive for $Y$, the differences of $Y$ across treatment groups, and the conditional treatment assignment probabilities.
The IATE is identified as the difference of two outcome regressions $\mu_{d}(x)$ in the different treatment groups. The splitting criterion of the mcf is therefore based on minimizing the sum of the expected MSEs at a given point $x$ over all unique treatment combinations, i.e., all IATEs relating to different treatment pairs are of equal importance:
where $MS{E^d}(x) = E\left( {{{\hat \mu }_d}(x) - \mu _d^0(x)} \right)$ and $MC{E^{m,l}}(x) = E\left({{{\hat \mu }_m}(x) -\mu _m^0(x)}\right)\left( {{\hat\mu }_l}(x)\right. \allowbreak - \allowbreak \left. \mu _l^0(x)\right)$.
For finding splits based on this criterion, the $MSE^{d}$ and $MCE^{m,l}$ elements need to be estimated in the daughter leaves. Let $N_{S(x)}^d$ denote the number of observations with treatment value $d$ in a certain stratum (leaf) $S(x)$, defined by the values of the covariates $x$. Then a `natural' choice to estimate the $MSE^{d}$'s in leaf $S(x)$ is:
where $\hat y_{S(x)}^d$ is an average of the observed outcomes in treatment $d$ in leaf $S(x)$. To compute the $MCE$, the mcf uses the closest neighbour (in terms of similarity w.r.t. $x$) available in the other treatment (which is denoted by $\tilde y_{(i,d)}$ below).\footnote{Implementational details on finding the closest neighbours for the MCE estimation are in (ref).}
Analysing the estimator of $\overline{MSE}_x$ reveals that its minimization favours splits maximizing the differences of all $\hat y_{S( \cdot )}^d$ between the daughter leaves. Additionally, the splits favour large differences between $\hat y_{S( \cdot )}^m$ and $\hat y_{S( \cdot )}^l$ in the same leaf. In case of a binary treatment, this is equivalent to maximising treatment effect heterogeneity, as in \citetalias{wager2018estimation} and \citetalias{athey2019grf}, under no selection into treatment. However, for multiple treatments, the mcf asymptotically targets treatment effect heterogeneity across daughter leaves while simultaneously increasing treatment outcome variability within leaves. Overlooking this dual focus may lead to splits that, while maximizing treatment effect heterogeneity, fail to minimise the sample MSE of parameters within leaves. Thus, mcf utilizes a splitting strategy that jointly considers all treatment combinations.
Although asymptotically confounding should be taken care of by the splitting rule derived above, the mcf adds a penalty term to the splitting criterion as an additional safeguard. As before, let denote $S(x')$ and $S(x'')$ denote the values of the covariates in the daughter leaves resulting from splitting some parent leaf. \citetalias{lechner2018modified} proposes to add the following penalty to a combination of the two `final' MSEs in the daughter leaves: $$penalty(x',x'')= \lambda \left\{ {1 - \frac{1}{M}\sum\limits_{d = 0}^{M - 1} {{{\left( {P(D = d| X \in S(x')) - P(D = d| X \in S(x''))} \right)}^2}} } \right\}, \lambda >0.$$ The probabilities, local estimates of the propensity score, are computed as relative shares of the respective treatments in the potential daughter leaves. The penalty term is zero if the split leads to a perfect treatment prediction of the treatment probabilities. It reaches its maximum value, $\lambda$, when all probabilities are equal. Thus, the algorithm prefers a split that is also predictive for $P(D = d\left| {X = x} \right.)$. Of course, the choice of the exact form of this penalty function is arbitrary. Furthermore, there is the issue of how to choose $\lambda$ (without expensive additional computations) which is taken up again in (ref).
Each tree is built on a subsample drawn from one half of the data, while the other half serves for the estimation of the IATE at point $x$. We name this procedure `two-sample honesty' as data points used to place the splits (“split-construction data”) will not be used to estimate the treatment effects (“leaf-populating data”) and vice versa.\footnote{See Figure (ref) in (ref) for a graphical representation.} The final estimator based on $B$ trees takes a form: $$\widehat {IAT{E^{mcf}}}(m,l;x) = \hat \theta _{m,l}^{mcf}(x) = \frac{1}{B}\sum\limits_{b = 1}^B {\sum\limits_{i = 1}^{N/2} {w_{i,b}^{mcf}({d_i},{x_i};x,m,l){y_i} = } \sum\limits_{i = 1}^{N/2} {w_i^{mcf}({d_i},{x_i};x,m,l){y_i}} }, $$ where $w_{i,b}^{mcf}(d_i, x_i; x, m, l) = \left( \frac{\underline{1}(d_i = m)}{N_{S_b(x)}^m} - \frac{\underline{1}(d_i = l)}{N_{S_b(x)}^l} \right) \underline{1}(x_i \in S_b(x))$ represent the weights for the $IATE(x)$ estimate in tree $b$. This estimator is differencing the average outcomes of the treatment groups evaluated on the “leaf-populating” subsample in each tree leaf containing point $x$, ${S_b}(x)$, and subsequently averages these differences across trees. As in the grf, the mcf forest weight $w_{i}^{mcf}(d_{i},x_{i};x,m,l)$ captures how important observation $i$ is for estimation of the parameter of interest at point $x$. While the mcf weights weigh the observed outcome directly, grf applies the forest weights in a locally weighted moment function. There is also a difference in the number of weights estimated in each method due to different honesty concepts. For a given $B$ and $N$, grf estimates $N$ weights on approximately $B/4$ trees, while mcf estimates $N/2$ weights on $B$ trees as illustrated in (ref) in Figures (ref) and (ref). Thus, \emph{grf} has an advantage in smoothing over more observations, while \emph{mcf} has an advantage in the precision of the weights, potentially influencing the finite sample properties of both methods.
The mcf estimators of GATE and ATE are averages of corresponding group or sample $\widehat {IAT{E^{mcf}}}(m,l;{x_j})$'s. By construction, they take the form of weighted observed outcomes. The mcf GATE estimator incorporates into the weighted average not only observed outcomes of the observations in the pre-specified group but potentially also observed outcomes of close neighbouring groups \footnote{This in general depends on whether variables $Z$ are among the splitting variables, type of $Z$, and the sample size. E.g., in settings with discrete $Z$ and when some trees do not (yet) split on all groups defined by a single value of $Z$, the observed outcomes from other groups that populate the same final leaves as $x_{j}$'s from group $z$ are incorporated into the weighted average. When all trees split on all groups, only observed outcomes from group $z$ are incorporated into the weighted average.}, defined by the non-zero forest weights in the IATEs.
(ref) contains the (so-far missing) theoretical guarantees for the mcf. Under similar assumptions as in \citetalias{wager2018estimation}, the mcf estimators of the IATEs, GATEs, and ATEs are consistent and pointwise asymptotically normal. The mcf inference procedures exploit the fact that mcf predictions, like for many Random Forest based estimators, can be computed as weighted means of $Y$. The details of this weights-based inference are presented in (ref), while (ref) explains some of the difference to the grf at the implementation level.
This section summarizes the statistical guarantees of the three CCMLs for estimation of ATE, GATE, and IATE. For the ATE, dml and grf provide the strongest guarantees, as they are based on the DR efficient score, and are expected to perform well in terms of bias, MSE and inference in the simulations.\footnote{Assuming a certain quality of the nuisance parameter estimators as mentioned in (ref) and (ref)} The mcf is also expected to perform well but with a somewhat larger MSE as it does not necessarily reach the efficiency bound.
Regarding GATE, all methods guarantee consistency and asymptotic normality. dml and grf differ from the mcf in the utilization of the observed data points. While dml and grf compute the average of estimated components of the DR score for observations in the corresponding group of interest, mcf builds a weighted average of $Y$'s, potentially incorporating information from groups close to the group of interest as determined by the forest structure, as mentioned in (ref). This will affect the efficiency of the GATE estimators in different scenarios as captured in the following proposition.
For the IATE, there are no specific statistical guarantees for dml. There is an upper bound on the error of the DR-learner providing a general guarantee. grf and mcf guarantee pointwise consistency and asymptotic normality and are both expected to yield lower MSEs in the simulations than dml-based DR-learner. For grf, the convergence rate is known and correspondingly slower due to the non-parametric nature of the estimator. Its speed depends on the subsampling rate which is implicitly restricted by the forest parameters (see Theorem 5 in \citetalias{athey2019grf}). All the above-mentioned conjectures are captured in Table (ref).
In this section, we compare the finite-sample performance of CCMLs in estimating average effects across three specified aggregation levels in a Monte Carlo simulation. This comparison includes evaluating both the precision of point estimates and the robustness of inference procedures.
In executing this Monte Carlo approach, we meticulously construct artificial datasets, varying numerous elements of the data generating process (DGP). These elements include the degree of selection into treatment (referred to as selectivity henceforth), the type and quantity of covariates, sample size, diverse functional forms, varying degrees of effect heterogeneity, the influence of covariates on outcomes and heterogeneity, the share of treated within the sample, and the number of treatments resulting in 77 different DGPs.\footnote{Table (ref) in (ref) gives the list of the scenarios investigated. (ref) collects the tables that contain the corresponding results.} Despite its significant computational demands, this approach offers a comprehensive exploration of various scenarios, surpassing the more limited scope of Monte Carlo studies typically employed in papers introducing a new methodology.
The subsequent subsection provides a concise overview of the simulation designs' components. The procedure involves initially generating random data, applying the different estimators to this training dataset, predicting effects on a separate dataset of equivalent size derived from the same DGP, saving the results, and repeating these steps $R$ times. Subsequently, we calculate a range of performance metrics that capture different dimensions of the accuracy and reliability of both estimation and inference processes.
The simulated (i.i.d.) data consist of the covariates, the treatment, and the potential outcomes. The simulation of these components will be discussed in turn (for additional details see (ref)).
The simulation involves generating a range of 10 to 50 independent covariates ($p=10$, 20, 50). These covariates are normally, uniformly, or binomially (dummy variables) distributed, or as combinations of the three types. Among these $p$ covariates, the first $k$ covariates influence both the selection into treatment and the potential outcomes. The effect of these $k$ covariates is modelled to decrease linearly, ranging from 1 to $1/k$. The base specification consists of 20 uniformly and normally distributed covariates, half of them relevant.\footnote{When covariates are simulated from both, the normal and uniform distributions, the 1$^{\text{st}}$, 3$^{\text{rd}}$, 5$^{\text{th}}$, ... covariates are drawn from the uniform distribution.}
The selection process is based on a linear index function of the $k$ relevant covariates plus noise. The quantiles of this index function are used to generate the treatments. The base specification considers cases with random selectivity (experiment), medium selectivity (true $R^2$ of about 10%), and strong selectivity (true $R^2$ of about 42%). Cases of two and four treatments with equal as well as asymmetric treatment shares are considered. The base specification consists of 2 treatments with equal treatment shares. Figure (ref) shows the distribution of the true treatment probabilities for the strong and medium selectivity case with two treatments and 50% treatment share.
The non-treatment potential outcomes are obtained by simulating the expected non-treatment potential outcome as a sine-function of the linear index of the covariates plus noise. The relevance of the sine-function relative to the noise level is varied in the simulations from cases with true $R^2$ of 0 to 45% (base specification: 10%). The potential treatment outcomes are obtained by adding simulated IATEs plus noise to the expected non-treatment potential outcomes.\footnote{The 2 (2 treatments) or 4 noise terms (4 treatments) used to simulate potential outcomes are independent of each other.} The IATEs are generated as functions of the linear index or as a smooth function of the first two, most important, covariates. In the former case, a linear function, a logistic function, and a quadratic function is specified. Finally, the smooth function approach follows closely the specification of \citetalias{wager2018estimation}. In all these cases the ATE is close to one. The case of zero IATEs, and thus zero ATE, is considered as well. The specification of the IATE as a smooth function is used in the base specification presented in (ref).\footnote{The base specifications with different degrees of selectivity and types of IATEs are captured in rows 1-15 in Table (ref) in (ref).} The smooth function, when plotted as a function of the linear index, exhibits a pattern that resembles a step function. Figure (ref) shows the expected and realised potential outcomes and their relation to the linear index. Thus, it summarizes main properties of the base specification following \citetalias{wager2018estimation}, as well as of the linear specification of the IATE. Similar plots for the other specifications can be found in (ref).
The GATEs depend only on the first covariate, which is uniformly distributed if there are uniformly as well as normally distributed covariates in the DGP. This continuous covariate is split in groups with the same expected sizes and corresponding indicator variables are created. Thus, these dummy variables have no direct effect on the DGP. They are passed to the GATE estimators.
A final aspect of the DGP is the sample size. We mainly consider sample sizes of 2'500 and 10'000 as compromise between computational costs and practical relevance.\footnote{$N$=40'000 is also considered, but, for computational reasons, only for one specification.} To keep the noise from the simulations on the performance measures stable (at least for the estimators of ATE and GATEs, which can be expected to show $\sqrt N - $convergence), the number of replications, $R$, declines at the same rate as the sample increases. Thus, the results for $N$=2'500 are based on $R$=1'000, while $R$ declines to 250 for $N$=10'000.
In our simulation study, it's crucial to highlight that we adopt a methodology designed to replicate repeated sampling inference more closely, especially in scenarios involving stochastic covariates. To achieve this, the samples used for computing the out-of-training-sample effects are freshly drawn from the same Data Generating Process (DGP) for each replication. Additionally, these samples are matched in size to the training sample. This approach contrasts with the method used in knaus2021machine, where the sample for calculating effects is drawn once prior to the initial replication and remains constant across subsequent replications, even when similar DGPs are employed.
The methodology chosen for our study, however, comes with certain trade-offs. Notably, for each IATE as a function of covariate values ($IATE(x)$), there exists only a single true and estimated value across all replications. This uniqueness arises because all or at least some covariates are continuous, leading to a scenario where each specific covariate value appears only once in the simulations. Consequently, for the IATEs this aspect of our simulation methodology precludes the possibility of estimating moments other than the expectation. This is an important consideration to bear in mind when interpreting the results and the applicability of our findings.
Below, we present the results of the Modified Causal Forest and the Generalized Random Forest, both with the outcome variable centred prior to estimation, mcf-cent and grf-cent, as well as of Double/debiased Machine Learning with normalized weights, dml-norm. For IATE estimation, we also include a more efficient version of the centred mcf for which no inference is available (see below), mcf-cent-eff. The tables in Appendix C contain additional results for the uncentred mcf and \emph{grf}, standard (non-normalized) \emph{dml}, and \emph{OLS}. We describe the implementation of these estimators briefly in this section and refer the interested reader to (ref) for details.
While (ref) details the implementation of the mcf further, at least 2 points merit some more discussion. The first point concerns the penalty parameter, $\lambda$. In the simulations, $\lambda$ is set equal to $Var(Y)$. $Var(Y)$ corresponds to the MSE when the effects are estimated by the sample mean without any splits. Thus, it provides some ad-hoc benchmark for plausible values of $\lambda$. Generally, decreasing the penalty increases biases and reduces variances, et vice versa.\footnote{The sensitivity of the centred version of the mcf to the choice of the penalty parameter seems to be generally very low.} The simulations below will show that biases are more likely to occur when selectivity is strong. Thus, if a priori knowledge about the importance of selectivity is available, then the researcher might adjust the penalty term accordingly.
Secondly, if inference is not a priority, like when using the IATEs as inputs into the training of an optimal assignment algorithm, the efficiency loss inherent in the mcf's two-sample honesty approach can be avoided by cross-fitting, i.e., by repeating the estimation with exchanged roles of the two samples and averaging the two (or more) estimates. However, in such a case it is unclear how to compute the weights-based inference for the averaged estimator mcf-cent-eff as the two components of this average are correlated.\footnote{For such a cross-fitted estimator (e.g., computed as mean of the single estimators), conservative inference could be obtained by basing inference on normality with a variance taken as average over the variances of the single estimations.} A similar efficiency loss can be avoided in the grf procedure by switching off the half-sampling when inference is not relevant. However, the current implementation in the grf package does not support this option.
As mentioned in (ref), grf performs local centring inside the Random Forest to remove confounding bias. The default grf local centring, labelled here as uncentred, refers to local centring that uses out-of-bag predictions of the outcome and treatment assignment from Random Forests trained on the full sample, to calculate the residuals. The simulation study includes an additional local centring strategy, labelled as centred. The difference between the centred and the uncentred mcf and grf is that the former uses a transformed outcome variable. This transformation subtracts a Random Forest prediction of $E(Y| X)$ from the observed $Y$ (obtained with 5-fold-cross-fitting), and thus purges them from much of the influence of $X$.\footnote{Note that the \emph{grf} package does not provide centring of this kind. We added this estimator as the original version used in the \emph{grf} package performed poorly in many of our DGPs for reasons discussed in (ref).} For \emph{grf}, GATEs and ATE are estimated via linear regression in which a variant of a local AIPW estimator is regressed on group indicators or a constant.
Compared to the standard dml, the normalized dml is more robust to extreme values of the estimated propensity score. It is obtained by normalizing the weights that are implicit in the dml scores. While it is straightforward to use dml for ATE, it is less straightforward to use them for heterogeneity estimation. Here, GATEs and IATEs are obtained by regression-type approaches in which the estimated components of the dml scores serve as dependent variable. The GATEs are obtained as OLS-coefficients of a saturated regression model with the indicators for the groups defined by the discrete variable $Z$ as independent variables. IATEs are computed by using $X$ as independent variables either in a regression Random Forest or in an OLS regression. When OLS is used, inference is based on the heteroscedasticity-robust covariance matrix of the corresponding coefficients. No inference is obtained for the Random Forest based IATEs.
All estimators are not tuned as computational costs would be prohibitive given the already extensive simulations. Instead, the default values provided in the respective software packages are used.
The main performance measures are the biases of the effects and their standard errors, the standard deviation of the effects, the mean absolute error, and the root mean squared error (RMSE) of the effects. The coverage probability (CovP) for the 95% confidence interval is reported to gauge the quality of the inference.\footnote{In addition to the CovP for the 80% interval, the tables in the (ref) also report the skewness and excess kurtosis of the estimators, which are, however, in a large majority of cases in the `normal', unproblematic ranges.}
Note that, as mentioned above, the standard deviation, and thus the bias of the estimated standard error, cannot be computed for the IATE due to the simulation design. Whenever several parameters are involved (as for the GATEs and IATEs), the performance measures are computed for each parameter and then averaged.
In this subsection, we analyse dml-norm, grf-cent, and mcf-cent for the \citetalias{wager2018estimation} base specification Data Generating Processes (DGPs) and sample sizes of $N$=2,500 and $N$=10,000 (when there is no confusion, in this subsection, we will drop the -norm and -cent ending in the text). Additionally, the degree of selectivity within the data is varied.
Since it turned out that the relative performance of the estimators does not only depend on the strength of selection into treatment, but also on the specific parameter to be estimated, the results for the ATE, the GATEs, and the IATEs are discussed in turn.
Table (ref) shows the results for the ATE. If there is no selectivity, like in an experiment, we obtain the expected results: the efficient estimators, dml and grf, are similar and outperform the mcf, at least when selectivity is not too strong. All estimators are essentially unbiased, and empirical coverage is close to the nominal level. Furthermore, when the sample size quadruples, the standard deviation and RMSE halve, which is indicative of $\sqrt N - $convergence.
The analysis reveals a notable trend with increasing selectivity: the performance of the normalized double/debiased machine learning (dml-norm) estimator deteriorates, particularly in terms of bias. This leads to a large increase in the RMSE.\footnote{As can be seen from Table (ref) in (ref), this behaviour appears for the not normalized dml as well and is driven mainly by the bias of the point estimator.} One plausible explanation is that (despite the normalisation) the double-robust score becomes more problematic when propensity score values become more extreme. A similar increase in bias is visible for the grf, although it is not as extreme as for dml. The resulting bias impacts also their coverage rates, which fall to very low levels. These issues appear to a much lesser extent for the mcf.
One summary coming from this table is that mcf appears to be more robust to stronger selectivity at the cost of some additional RMSE when selectivity does not matter much. The results in (ref) show that these performance patterns with respect to the degree of selectivity also appear for linear IATEs and non-linear IATEs. However, in the case of quadratic IATEs it turns out that the uncentred (!) grf outperforms the centred grf in terms of bias even for strong selectivity. One possible explanation for this surprising behaviour is that the Random Forest estimator used to do the centring does not capture the shape of $E(Y| X=x)$ well.
Fixing selectivity to medium levels and varying the other parameters of the DGP ((ref)), suggests that the patterns observed in Table (ref) qualitatively appear almost in all other DGPs as well. The most remarkable case is when covariates are more important for the non-treatment potential outcome (Table (ref), Table (ref)): all estimators become biased, which increases the RMSE and leads to coverage far below nominal levels.
Table (ref) shows that the relative performance of the different estimators for the GATEs depends not only on the strength of selectivity, but also on the number of groups for which a GATE is computed for as captured in Proposition (ref). Table (ref) shows the results for the cases of 5 and 40 groups, while Tables (ref) to (ref) in (ref) show the intermediate cases of 10 and 20 groups as well.
Concerning the point estimate, mcf outperforms grf and dml once there are 10 and more ($N$=2'500) or 20 and more groups ($N$=10'000), independent of the strength of selectivity. It is not surprising that mcf performs better for strong selectivity: the previous section showed that for strong selectivity the mcf dominates even for a GATE with only one group, i.e., the ATE.
We observe a differential dependence of the standard deviation on the number of groups: its increase is much slower for the mcf than for dml and grf. The reason is the way the GATE estimators are constructed by the different methods. As mentioned above, dml and grf are averaging double-robust scores within the cells of the discrete $Z$. Since the estimators are $\sqrt N - $convergent, we expect that when the sample is reduced to one quarter of the original sample, the standard deviations of such estimators double. As the cells defined by the group variable in the DGPs are of approximately equal size, the number of groups and observations per group are inverse proportionally related. The resulting doubling of the standard deviation of the grf and \emph{dml} estimators is exactly what is observed in Tables (ref) to (ref) when comparing results for 5 vs. 20 groups, or 10 vs. 40 groups, respectively.
The mcf aggregates IATEs within these cells. However, as the IATEs also potentially use “leaf-populating data” outside of these cells, the standard deviation of the mcf increases much slower than for grf and dml.\footnote{The aggregation of IATEs that potentially weigh information from observations belonging to other groups into a GATE resembles smoothing across (adjacent) categories in the context of categorical regressors. heiler2021shrinkage show that optimal smoothing parameters in a non-parametric regression involving categorical regressors do not vanish asymptotically. Furthermore, they prove that the variance of the smoothed estimator is a weighted sum of the asymptotic variances from other categories. This finding provides a plausible rationale for the observed steadiness in the variance of the mcf estimator across varying group sizes.} In fact, we observe that when the number of groups increases further, the RMSE of the mcf approaches the one of the \emph{mcf} IATEs from below, while the RMSE of \emph{dml} and \emph{grf} increases substantially and sometimes exceeds the RMSE of the IATEs which are of much higher dimension.
Concerning inference, the findings are a bit more pronounced: dml and grf have the correct coverage rates for no and medium selectivity even for 40 groups, while the coverage rates for mcf are too low.\footnote{Note that when selection into treatment gets stronger, splitting for the mcf (and grf) stops earlier, as leaves must have both treated and controls after a possible split, thus bias increases.} The figures contained at the end of Tables (ref), (ref), and (ref) in (ref) indicate that this problem comes mainly from a bias of the GATEs with the smallest true values, while for the other GATEs coverage is close to the nominal level. In fact, the good performance in terms of RMSE and the problems with coverage are two sides of the same coin. Due to their aggregation from IATEs that are averaged over trees, mcf estimators share the strengths and the weaknesses of many smoothing methods: the variance is reduced at the cost of an increasing bias. To improve on the inference in finite samples, the simulations indicate that it needs some debiasing. This is however beyond the scope of this paper.
\FloatBarrier
Table (ref) contains the results for the IATEs, averaged over all $N$ IATEs. As inference is usually not the main goal when computing IATEs, for the centred mcf the more efficient version is included as well.
For the point estimates the ordering is clear-cut. grf and mcf are similar and outperform both dml versions (ols, rf). In many cases, the efficient centred mcf (\emph{mcf-cent-eff}) performs best in terms of RMSE. The only exception to this rule seems to be, again, the case of a quadratic function for the IATE, in which all methods seem to do (almost) equally bad with large RMSEs and large mean absolute errors. With respect to the other estimators, in a couple of cases OLS outperforms the more sophisticated CML estimators (but substantially underperforms in many other scenarios).
For inference, the results show that all estimators have coverage probabilities that are substantially below their nominal levels. However, given the level of granularity of the IATEs, this finding is of course not surprising. At least for the mcf, we conjecture that this problem comes from the biases of the individual IATEs. It is less likely that it comes from a too small estimated standard error, because weights-based standard error estimation is performed in a very similar way as for the ATE and the GATEs, in which it turned out to be almost unbiased.
dml performs particularly well when the target is low dimensional, such as the ATE and GATEs with few groups, and selectivity is not too strong. If selectivity is strong or the number of groups becomes too large, then its performance quickly deteriorates.\footnote{A potential alternative dml-based estimator in the strong selectivity setting could be the automatic dml that would estimate inverse probability weights instead of propensity scores as nuisance parameters.}
grf shows a similar good performance for the low dimensional parameters, as well as a similarly bad performance for GATEs with many groups.\footnote{Additional smoothing in a setting with many groups could be considered in empirical settings for dml-based and grf -based GATE estimators to further improve their RMSE.} However, it performs well for IATEs. It is also less affected by strong selectivity than dml. These results however only hold for the modified version of the grf in which the outcome variable in the training data is explicitly centred before it enters the grf algorithm. The original uncentred version of \emph{grf} as suggested in athey2019grf and implemented in their package underperforms in many scenarios due its bias problem (as does the uncentred version of the \emph{mcf}).\footnote{As described in (ref), labels centred and uncentred refer to whether the outcome variable enters the estimation procedure residualized by $K$-fold cross-fitting or not, respectively. In both cases, the \emph{grf} procedure still performs local centring using out-of-bag predictions as described in (ref).}
Compared to dml and grf, the mcf generally shows a robust and competitive performance in many scenarios. The price to pay for this robustness and competitiveness is the somewhat higher standard deviation for the very low dimensional parameters like ATE and GATEs. Depending on the sample and the resulting precision of the estimators, this price may well be worthwhile paying. A further important advantage of the mcf is that its estimates are internally consistent over aggregation levels as ATE and GATEs are computed as averages of IATEs. This advantage becomes particularly apparent when the number of groups for the GATEs increases. While dml and grf have substantially higher RMSEs than for the more fain-grained IATEs, the uncertainty (and bias) in the \emph{mcf} estimators smoothly increases with the number of groups, until it reaches the level of the most fine-grained heterogeneity parameter, the IATE.
Estimation of causal effects at different levels of granularity is of great importance for informed decision-making and tailored interventions. Lower levels of granularity capture the effect of a policy or intervention on a large population, guiding decisions on policies that cannot be targeted at individuals but must be deployed universally. Higher levels of granularity capture effects at a group or individual level that can serve for decisions how policies can be tailored more individually or for better understanding of the effects of large-scale policies at more granular levels. The complexity of such interventions requires methods that can estimate fine-grained heterogeneities of causal effects flexibly, such as some Causal Machine Learning (CML) methods.
In this paper, we investigate such CML methods subject to the restrictions that (i) they provide estimators of the causal effects at all aggregations levels, (ii) are essentially non-parametric, (iii) can predict IATE based on covariates only and (iv) that they allow for classical repeated sampling inference. Ideally, they are also internally consistent that aggregation of lower-level effects lead to the higher-level effects. Double/debiased Machine Learning (dml), the Generalized Random Forest (grf), and the Modified Causal Forest (mcf) fulfil these criteria and thus belong to the group of estimation methods which we call `Comprehensive Causal Machine Learners' (CCMLs).
Here, we describe these estimators and their proven theoretical guarantees. For dml and the grf, they are already known, but not so for the mcf. Therefore, we explicitly provide them. The large-scale simulation study reveals scenarios in which the methods perform well. dml with normalized weights performs well in terms of RMSE and coverage probability when the target is low-dimensional, i.e., ATE and GATEs with few groups, and when the selection into treatment is not strong. grf shows similar behaviour as dml, however its performance for IATEs in terms of RMSE is much better. \emph{mcf} has similar performance as \emph{grf} in case of IATEs and outperforms \emph{dml} and \emph{grf} in scenarios with many GATEs or strong selection.
The results of the simulation study offer several practical recommendations. For low-dimensional targets when selection into treatment is moderate, dml is preferred including statistical inference. For IATEs, Causal Forest based methods perform well in terms of point estimation. When inference is not a priority, the more efficient version of mcf is recommended to estimate IATEs. For large groups, mcf is recommended for point estimation of GATEs. When sample is large enough and slight loss of efficiency is not detrimental, mcf can be used to estimate all effects due to its robustness to strong selection into treatment and large GATE groups. When internal consistency of the effects is important, only mcf can guarantee it due to its aggregation strategy.
The practical use of these three CCMLs is supported by the availability of well-maintained software packages: dml is available as Python and R packages, the grf is available as an R package, and the mcf is available as a Python package.\footnote{As examples can serve: bach2022doubleml, bach2024dmlr, and knaus2022double for dml; athey2019estimating for grf; and bodory2022high for mcf. The R packages can be downloaded from CRAN. The Python packages are available on PyPI.} On a practical note, the results in this paper indicate that users of the \emph{grf} package are encouraged to perform local centring of the outcome variable via $K$-fold cross-fitting to remove any potential bias (the same holds for \emph{mcf}, where this step is already implemented in the package).
The simulation study also points to topics for further research. For example, the sensitivity of the dml estimators to strong selectivity opens a topic of how to make them more robust to extreme propensity scores without impairing its efficiency and statistical inference. Due to the smoothing character of Causal Forests, in finite samples mcf estimation of GATEs and IATEs balances the bias-variance trade-off yielding low RMSEs but impairing the coverage probability due to a small bias and low variance. Finding a way how to further de-bias the estimates would improve coverage probability in finite samples.