The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
113,031 characters
Comprehensive Causal Machine Learning
\begin{titlepage}
\title{\vspace{-1.5cm} \LARGE \textbf{Comprehensive Causal Machine Learning} \thanks{\scriptsize Results from two research projects that were part of the National Research Programmes ``Big Data'' (NRP 75, \url{www.nrp75.ch}, grant number 407540\_166999) and ``Digital Transformation'' (NRP 77, \url{www.nrp77.ch}, grant number 407740\_187301) of the Swiss National Science Foundation (SNSF) were the basis of this paper. The theoretical part on the mcf (contained in \ref{appendix-a}) is a revised version of the theoretical part of our unpublished ``Modified Causal Forest'' paper. We thank GPT-4 for some limited assistance in editing, and Phillip Heiler, Federica Mascolo, Fabian Muny, and Hannah Busshoff for helpful comments and suggestions. The paper was presented at research seminars at the Universities of Bern and the Collegio Carlo Alberto in Turino, at an invited session of the Annual Meeting of the German Economic Association in Berlin, at CompStat in Giessen, at American Causal Inference Conference in Seattle, at COMPIE in Amsterdam and at the Workshop on Causal Inference and Machine Learning in Groningen. We are grateful to participants for helpful remarks and interesting discussions.}}
\author{\Large \large Michael Lechner\thanks{\scriptsize Michael Lechner is also affiliated with �rebro University, CEPR, London, CESIfo, Munich, IAB, Nuremberg, IZA, Bonn, and RWI, Essen.} \\ \and \large Jana Mareckova}
\date{\large{{\today}}}
\maketitle
\thispagestyle{empty}
\begin{center}
\footnotesize
\textbf{Abstract}
\end{center}
\vspace{-25pt}
\begin{abstract} \singlespacing \footnotesize
Uncovering causal effects in multiple treatment
setting at various levels of granularity provides substantial value to
decision makers. Comprehensive machine learning approaches to causal
effect estimation allow to use a single causal machine learning approach
for estimation and inference of causal mean effects for all levels of
granularity. Focusing on selection-on-observables, this paper compares
three such approaches, the modified causal forest (\emph{mcf}), the
generalized random forest (\emph{grf}), and double machine learning
(\emph{dml}). It also compares the theoretical properties of the
approaches and provides proven theoretical guarantees for the
\emph{mcf.} The findings indicate that \emph{dml}-based methods excel
for average treatment effects at the population level (ATE) and group
level (GATE) with few groups, when selection into treatment is not too
strong. However, for finer causal heterogeneity, explicitly
outcome-centred forest-based approaches are superior. The \emph{mcf} has
three additional benefits: (i) It is the most robust estimator in cases
when \emph{dml}-based approaches underperform because of substantial
selection into treatment; (ii) it is the best estimator for GATEs when
the number of groups gets larger; and (iii), it is the only estimator
that is internally consistent, in the sense that low-dimensional causal
ATEs and GATEs are obtained as aggregates of finer-grained causal
parameters.
\end{abstract}
\vfill
\singlespacing \footnotesize
\noindent \textbf{Keywords:} Causal machine learning, statistical learning, conditional average treatment effects, individualized treatment effects, multiple treatments, selection-on-observed-variables
\vspace{5pt}
\noindent \textbf{JEL classification:} C21, C87
\vspace{5pt}
\noindent \textbf{Correspondence to}: Michael Lechner or Jana Mareckova,
Professors of Econometrics, Swiss Institute for Empirical Eco�nomic
Research (SEW), University of St.~Gallen, Switzerland,
[email removed], [email removed], \url{www.sew.unisg.ch}.
\end{titlepage}
\newpage \pagenumbering{arabic}
\section{Introduction}\label{introduction}
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 \emph{effects} of policies, medical
treatments, marketing campaigns, business decisions, etc. \citep[see, e.g.,][]{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
\emph{effects}. Building on these conditions, causal problems can be
transformed into specific prediction problems \citep[see, e.g.,][]{imbens2009recent}
for which ML methods are suitable \citep[see, e.g.,][for an overview on ML methods]{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 \citep[for an overview, see][]{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
\emph{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 \emph{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.
\citet{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
(\emph{dml}; \citetalias{chernozhukov2018dml}), the Generalized Random Forest (\emph{grf}; \citealp{athey2019grf}), and the Modified Causal Forest
(\emph{mcf}; \citealp{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,
\emph{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 \emph{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,
\emph{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
\emph{mcf} offers three additional benefits: (i) it is the most robust
estimator even for the ATE in cases when \emph{dml}-based approaches
underperform because of substantial selection into treatment that needs
to be corrected for; (ii) it outperforms both \emph{dml} and \emph{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: \autoref{2-review} 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. \autoref{3-causal-framework} introduces the notation of the potential outcome
model, defines the parameters of interest, and states the identifying
assumptions imposed under unconfoundedness. \autoref{4-ccml} discusses the
three CCMLs. The large-scale simulation study is summarized in
\autoref{5-mcs}. \autoref{6-conclusion} concludes. \ref{appendix-a} provides proven theoretical
guarantees for the \emph{mcf}, its inference procedure, a comparison to
the \emph{grf}, as well as some new results for the \emph{grf}.
Implementational details of the simulation study are collected in
\ref{appendix-b}, while \ref{appendix-c} holds its detailed results.
\section{Review of causal machine learning estimation}\label{2-review}
\subsection{Average treatment effect}\label{21-ate}
For causal inference under the unconfoundedness assumption \citep[see][]{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 \citep[see, e.g.,][]{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 \citep{westreich2010propensity}
or Boosting Models \citep{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 \citet{imai2014covariate}.
Alternatively, \citet{graham2012inverse} improve IPW estimation by
covariate balancing within the empirical likelihood framework, without
explicit modelling of the propensity score.
\citet{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. \citet{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
\autoref{3-causal-framework}, 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) \citep{hill2011bayesian} and Ensemble
Methods \citep{austin2012using} were introduced to obtain regression-based
estimators of ATEs. However, issues like regularization bias \citep{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 \citet{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., \citet{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 \citep{chernozhukov2018dml}.
Such an orthogonality property is leveraged, e.g., in \citet{belloni2014inference}
for ATE estimation with binary
treatment and in \citet{farrell2015robust} for multiple treatments. As
Neyman-orthogonal scores mitigate first order bias arising from ML
estimation of nuisance parameters \citep[p.~12]{bach2024dmlr},
Neyman-orthogonality emerged as a key element for
the Double Machine Learning (\emph{dml}) framework introduced in
\citet{chernozhukov2018dml}. \emph{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. \emph{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. \citet{chernozhukov2018dml} apply the
\emph{dml} framework to ATE estimation using the DR score of
\citet{robins1995semiparametric}\footnote{Note that the DR score of
\citet{robins1995semiparametric} is Neyman-orthogonal.} and derive
its statistical properties.
Another approach for ATE estimation is Targeted Minimum Loss Estimation
(\emph{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 \citep{van2006targeted}.
Among various \emph{tmle} adaptations\footnote{See \citet{Laan2018}
for a comprehensive overview.}, cross-validated
\emph{tmle}, \emph{cv-tmle}, proposed by \citet{Zheng2011}, shares
considerable similarities with \emph{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 \emph{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 \citet{knaus2021double}.
Automatic Debiased Machine Learning, introduced in \citet{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
\citet{robins1995semiparametric}. The estimation of the Riesz representer and
nuisance parameters requires cross-fitting as in \emph{dml}. The method
is shown to yield a consistent and asymptotically normal estimator and
provides a consistent asymptotic variance estimator.
\subsection{Conditional average treatment effect}\label{22-cate}
\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
\citep{wager2018estimation,athey2019grf,lechner2018modified}, along with
\emph{dml}-based estimators yielding estimated components of the DR
efficient score that are further regressed on the covariates as first
proposed by \citet{Laan2006}. Alternative approaches include
Meta-learners, which use multiple machine learning algorithms as ``base
learners'' to estimate large-dimensional CATEs. The \emph{X-learner},
introduced by \citet{kunzel2019metalearners}, is effective in
cases with uneven treatment groups, while the \emph{R-learner},
developed by \citet{nie2021quasi}, leverages the \citet{robinson1988root}
transformation to estimate heterogeneous causal effects. When interested
in CATEs for specific subgroups, the above-mentioned quantities of the
\emph{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 \emph{tmle} ``one-group-at-a-time'' for CATE
estimation with discrete groups might be computationally costly. As a
remedy, they propose a \emph{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 \citet{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 \citet{chernozhukov2018dml}, averaging the
estimated components of individual DR efficient scores of \citet{robins1995semiparametric}
to estimate ATE in a \emph{dml} framework, spanned
further methodological contributions, particularly in discovering
heterogeneity along a chosen set of covariates, such as gender or age
group. \emph{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
\citet{semenova2021debiased}, \citet{zimmert2019nonparametric}, and \citet{fan2022estimation}.
For a higher dimensionality of the set of the
chosen covariates, \citet{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 \citet{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, \citet{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, \citet{athey2019grf}
introduce the Generalized Random Forests (\emph{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 \citet{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, \citet{lechner2018modified} introduces a Modified Causal Forest
(\emph{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.
\emph{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 \autoref{43},
for building the forest and the estimation of the effects. A
simulation study in \citet{lechner2018modified} shows superior performance of
\emph{mcf} to the Causal Forest of \citet{wager2018estimation}. \citet{bodory2022high}
show that the \emph{mcf} works in
applications replicating results from several papers in different
fields.
\subsection{Comprehensive treatment effect estimation}\label{23-ccml}
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{tab1} 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 \emph{tmle} and Meta-learners do
not have available statistical inference for all levels and that
\emph{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, \emph{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, \emph{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, \emph{grf} and \emph{mcf}, provide in their first stage CATE
estimates at the highest level of granularity. Using direct aggregation
step of these CATEs, \emph{mcf} leverages the weighted-outcome
representation of forest predictions, while \emph{grf} uses a variant of
the AIPW score to yield higher level aggregations via a regression step.
Table \ref{tab1} summarizes the estimation at different levels of granularity for
CCMLs and other methods.
\begin{table}[h]
\centering
\sffamily
\caption{Properties of commonly used CML methods for treatment effect estimation}
\label{tab1}
\small
\begin{tabularx}{\textwidth}{lX>{\centering\arraybackslash}X>{\centering\arraybackslash}X>{\centering\arraybackslash}Xccc}
\toprule
\multicolumn{2}{l}{\textbf{Approach}} & \multicolumn{3}{c}{\textbf{Parameters}} & \textbf{Statistical} & \textbf{CATE pre-} & \textbf{Internal} \\
& & & & & \textbf{inference} & \textbf{diction de-} & \textbf{con-} \\
& & \multicolumn{2}{c}{\textit{Dimension of CATE}} & \emph{ATE} & & \textbf{pends only} & \textbf{sistency} \\
& & large & low & & & \textbf{on covari-} & \\
& & & & & & \textbf{ates} & \\
\midrule
\multirow{3}{*}{\rotatebox{90}{\textbf{CCML}}} & \textbf{\emph{dml}} & Regress & Regress & Regress & Yes & Yes & No \\
& \textbf{\emph{grf}} & LocalGMM & Regress & Regress & Yes & Yes & No \\
& \textbf{\emph{mcf}} & MCF & Aggreg & Aggreg & Yes & Yes & Yes \\
\hdashline
& \textbf{Meta-L} & Meta-L & Aggreg* & Aggreg* & No** & Yes & Yes \\
& \textbf{\emph{tmle}} & TMLE & TMLE & TMLE & Yes (low-dimen.) & No & Yes \\
& \textbf{\emph{bcf}} & BCF & BayAggr & BayAggr & Bayesian & Yes & Yes \\
\bottomrule
\end{tabularx}
\vspace{-8pt}
\captionsetup{justification=justified,singlelinecheck=false}
\caption*{\footnotesize \textsf{Note: \emph{dml}: Double / Debiased Machine Learning; \emph{LocalGMM}: Local General Method of Moments; \emph{mcf/MCF}: Modified Causal Forest; \emph{bcf/BCF}: Bayesian Causal Forest; \emph{tmle/TMLE}: Targeted Minimum Loss Estimation; \emph{Meta-L}: Meta-learner. \\
\emph{Regress}: Regression-based estimation; \emph{Aggreg}: Aggregation of estimated IATEs; \emph{BayAggr}: Bayesian Aggregation. \\
* Averaging of the finest large-dimensional CATEs for Meta-learners was implemented in, e.g., \citet{salditt2024tutorial}. \\
** However, see recent developments in constructing valid confidence intervals for Meta-learners based on conformal prediction in \citet{alaa2024conformal}.}}
\end{table}
\subsection{Simulations and empirical work}\label{24}
\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,
\citet{knaus2021machine} look at the finite-sample
performance of selected machine learning estimators, covering \emph{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 \citet{knaus2021machine}
were performed, the \emph{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
\emph{grf} with local centering. \citet{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 \citet{alaa2017bayesian},
perform best in terms of MSE. The \emph{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,
\emph{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, \emph{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 \citet{wendling2018comparing} differs from the simulation in \autoref{5-mcs}
by using fixed covariates across replications. Despite this
difference, their results are in line with the results in \autoref{5-mcs}
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
\citep{davis2017using, bertrand2021workfare, Pytka2021Understanding,
knaus2022heterogeneous, cockx2023priority, burlat2024everybody},
education \citep{knaus2021double, farbmacher2021heterogeneous},
social experiments \citep{athey2019estimating, strittmatter2023value},
energy use \citep{knittel2021machine}, natural resource rents \citep{hodler2023institutions},
finance \citepalias{audrino2024does}, and
the dating market \citep{boller2021effect} to medicine \citepalias{langenberger2023exploring}.
Regarding the comprehensive nature of the \emph{dml} framework, \citet{knaus2022double}
demonstrates its flexibility in an active labour market programme
evaluation setting. \citet{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 \citet{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.
\section{Causal framework}\label{3-causal-framework}
\subsection{The potential outcome model}\label{31}
We use Rubin's potential outcome language \citep{rubin1974estimating}
to describe a multiple
treatment model under unconfoundedness, also referred to as
selection-on-observables or conditional independence \citep{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:
\begin{align*}
& IATE(m,l;x) = E({Y^m} - {Y^l}| X = x), \\
& GATE(m,l;z) = E({Y^m} - {Y^l}| Z = z) = \int{IATE(m,l;\tilde x,z){f_{\tilde X| Z = z}}(\tilde x)}d\tilde x, \\
& ATE(m,l) = E({Y^m} - {Y^l}) = \int {IATE(m,l;x){f_X}(x)}dx.
\end{align*}
The \textbf{I}ndividualized \textbf{A}verage \textbf{T}reatment
\textbf{E}ffects (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 \textbf{A}verage
\textbf{T}reatment \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 \citet{lechner2018modified} for more details. Note that \cite{abrevaya2015estimating}
and \citet{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).
\subsection{Identifying assumptions}\label{32}
The following set of assumptions identifies the causal effects discussed
in the previous section \citep[see][]{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 \citet{imbens2000role, imbens2004nonparametric}, or \citet{Lechner2001}.}
\begin{align*}
&Y^d \rotatebox[origin=c]{180}{$\ \boldsymbol{\Pi} \ $} D | X = x, & & & & \text{(CIA)} \\
&0 < P(D=d|X=x) = p_d(x), & & \forall x \in \chi, \forall d \in \{0,\ldots, M-1\}; & & \text{(CS)} \\
&Y = \sum\limits_{d = 0}^{M-1}\underline{1}(D=d)Y^d; & & & & \text{(Observation rule)}
\end{align*}
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 \emph{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 \citep[see][]{hurwicz1950}:
\small
\begin{align*}
IATE(m,l;x) & = E({Y^m} - {Y^l}| X = x) \\
& = E({Y^m}| X = x,D = m) - E({Y^l}| X = x,D = l)\\
& = E(Y| X = x,D = m) - E(Y| X = x,D = l) \quad {\forall x \in \chi ,\forall m \ne l \in \{0,...,M -1\}.}
\end{align*}
\normalsize
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.
\section{Comprehensive approaches for estimation and inference}\label{4-ccml}
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 (\emph{dml)} and Generalized Random Forest
(\emph{grf)} are reviewed. The third subsection focuses on the Modified
Causal Forest (\emph{mcf)}. While the basic principles of the \emph{mcf}
are introduced in \citet{lechner2018modified}, that paper does not contain explicit
theoretical guarantees. Therefore, \ref{appendix-a1} provides the proofs of
consistency and asymptotic Gaussianity properties for IATEs, GATEs, and
ATEs.
\subsection{Double Machine Learning}\label{41}
\emph{dml} estimation \citep{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 \citep[as in][]{neyman1959optimal}.
A Neyman-orthogonal score that identifies ATE, GATE, and IATE for
treatments $m$ and $l$ under assumptions outlined in
\autoref{32} is the DR score of \citet{robins1995semiparametric}:
\begin{align*}
& \psi_{m,l}^{dml}(O;{\theta_{m,l}},\eta_{m,l}^{dml}(X)) = \Gamma _{m,l}^{dml}(O;\eta_{m,l}^{dml}(X)) - {\theta_{m,l}}, \\
& \Gamma_{m,l}^{dml}(O;\eta_{m,l}^{dml}(X)) = \mu_m(X) - \mu_l(X) + \frac{\underline{1}(D=m)(Y-\mu_m(X))}{p_m(X)} - \frac{\underline{1}(D=l)(Y-\mu_l(X))}{p_l(X)},
\end{align*}
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:
\begin{align*}
& E\left({\psi_{m,l}^{dml}(O;\theta_{m,l}^0,\eta _{m,l}^{dml,0})}\right) = 0,& & \theta _{m,l}^0 = ATE(m,l), \\
& E\left({\psi_{m,l}^{dml}(O;\theta_{m,l}^0(z),\eta_{m,l}^{dml,0})\left| {Z = z}\right.} \right) = 0, & & \theta _{m,l}^0(z) = GATE(m,l;z), \\
& E\left( {\psi_{m,l}^{dml}(O;\theta_{m,l}^0(x),\eta_{m,l}^{dml,0})\left| {X = x}\right.} \right) = 0,& & \theta _{m,l}^0(x) = IATE(m,l;x),
\text{ for all } x \in \chi.
\end{align*}
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 \citet{chernozhukov2018dml}. The variance estimator of the
\emph{dml} can be computed as the sample variance of the estimated DR
scores.\footnote{See Theorem 3.2 in \citet{chernozhukov2018dml}, and
Theorem 1 and Remark 1 in \citet{bach2024dmlr}.} The
\emph{dml} estimator asymptotically reaches the semiparametric
efficiency bound of \citet{hahn1998role}.\footnote{See Theorem 5.1 in
\citet{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. \citet{semenova2021debiased}
provide statistical inference for the best linear
predictor in this case, a method that is implemented in \autoref{5-mcs}.
Standard errors are estimated via heteroscedasticity-robust standard
errors for asymptotically valid pointwise and uniform confidence bands.
In general, the results in \citet{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 \citet{zimmert2019nonparametric} and \citet{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 \autoref{44} in Table \ref{tab2}, this step is called a
smoothing step.} \citet{foster2023orthogonal} and \citet{kennedy2023towards} derive
error bounds for this two-step procedure, referred to as DR-learner in
\citet{kennedy2023towards}. Particularly, the findings in \citet{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: \citet{goller2023analysing} utilizes linear
regression, while \citet{knaus2022double} opts for a Random Forest approach.
\subsection{Generalized Random Forests}\label{42}
\defcitealias{athey2019grf}{ATW19}
\emph{grf} (\citealp{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. \emph{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 \emph{intercept} term $c(x)$ in \citetalias{athey2019grf}.} Note
that the score function in \emph{dml} is based on a fully nonparametric
DR score identifying treatment effects for all pairwise combinations of
treatments, while the score function in \emph{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 \emph{local} generalized method of moments, \emph{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 \emph{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:
\begin{align*}
\hat \theta^{grf}(x) & = \widehat {IATE^{grf}}(x) \\
& =\left(\sum\limits_{i =1}^N w_i^{grf}(x)\left(\tilde d_i - \bar{\tilde{d}}_w\right)\left(\tilde d_i - \bar{\tilde{d}}_w\right)'\right)^{-1} \left(\sum\limits_{i = 1}^N w_i^{grf}(x)\left(\tilde d_i - \bar{\tilde{d}}_w\right)\left(y_i - \bar y_w \right) \right),
\end{align*}
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 \citet{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 \emph{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{appendix-c1} 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 \citet{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 \cite{sexton2009standard}. \citetalias{athey2019grf} motivate it by the
observation that ``\textit{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
\citet{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{fig:a1} in \ref{appendix-a3}
captures the \emph{grf} procedure graphically.} Consistent estimation
of the variance term yields asymptotically valid confidence intervals
for the IATE at point $x$.
Like \emph{dml}, \emph{grf} offers comprehensive estimation of treatment
effects at all levels of granularity. However, unlike \emph{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, \emph{grf} takes a more direct approach.
Specifically, \emph{grf} estimates IATEs ``directly'', while ATE and
GATEs are estimated through additional regressions. The \emph{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''.}
\begin{align*}
\widehat {AT{E^{grf}}}(m,0) = \hat \theta _{m,0}^{grf} = & \ \frac{1}{N}\sum\limits_{i= 1}^N {\Gamma_{m,0}^{grf}} ({o_i};\hat \eta_{ATE(m,0)}^{grf}({x_i})), \\
{\Gamma_{m,0}^{grf}} ({o_i};\hat \eta_{ATE(m,0)}^{grf}({x_i})) = & \ \hat{\theta}_{m,0} + \frac{\underline{1}(d_i = m)}{\hat{p}_m(x_i)} (y_i - \hat{\mu}_m(x_i)) - \frac{\underline{1}(d_i = 0)}{\hat{p}_0(x_i)} (y_i - \hat{\mu}_0(x_i)), \\
&\text{where} \quad \hat{\mu}_0(x) = \hat{\mu}(x_i) - \sum_{d>0} \hat{p}_d(x_i) \hat{\theta}_{d,0}(x_i) \\
&\text{and} \quad \hat{\mu}_m(x_i) = \hat{\mu}_0(x_i) + \hat{\theta}_{d,0}(x_i).
\end{align*}
As default, the \emph{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{appendix-a4} show the equivalence of the \emph{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 \citet{robins1995semiparametric}
for the ATE. It can therefore be expected that the \emph{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{appendix-a4}). Consistent
estimation of the ATE based on the above \emph{grf} score therefore
requires consistent estimation of the propensity scores. \ref{appendix-a4}
also shows Neyman-orthogonality of the \emph{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 \emph{dml}. Inference for the ATE and GATEs is conducted
as in \emph{dml}.
\subsection{Modified Causal Forest}\label{43}
\defcitealias{wager2018estimation}{WA18}
\defcitealias{lechner2018modified}{L18}
\citet[abbreviated as L18 in this section]{lechner2018modified} builds
on \citet[abbreviated as WA18]{wager2018estimation} and introduces the Modified Causal
Forest estimator, \emph{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.
\emph{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. \emph{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, \emph{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 \emph{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:
\begin{align*}
&{\overline{MSE}_x} =
\sum\limits_{m = 1}^{M - 1}
{\sum\limits_{l = 0}^{m - 1}
{MSE\left( {\widehat {IATE}(m,l;x)}
\right)} } , \\
&MSE\left( {\widehat
{IATE}(m,l;x)} \right) = MS{E^m}(x) +
MS{E^l}(x) - 2MC{E^{m,l}}(x),
\end{align*}
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:
\begin{equation*}
\widehat{MSE_{S(x)}^d} =
\frac{1}{N_{S(x)}^d}
\sum_{i=1}^{N} \underline{1}(x_i \in S(x)) \underline{1}(d_i = d)
\left( \hat{y}_{S(x)}^d - y_i \right)^2,
\end{equation*}
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 \emph{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{appendix-b31}.}
\begin{multline*}
\widehat{MCE_{S(x)}^{m,l}} = \\
\frac{1}{N_{S(x)}^m + N_{S(x)}^l}
\sum_{i=1}^{N} \underline{1}(x_i \in S(x))
\left( \underline{1}(d_i = m) + \underline{1}(d_i = l) \right)
\left( \hat{y}_{S(x)}^m - \tilde{y}_{(i,m)} \right)
\left( \hat{y}_{S(x)}^l - \tilde{y}_{(i,l)} \right), \\
\tilde{y}_{(i,d)} =
\begin{cases}
y_i, & \text{if } d_i = d, \\
\tilde{y}_{(i,d)}, & \text{if } d_i \neq d.
\end{cases}
\end{multline*}
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 \emph{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, \emph{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 \emph{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$ (\emph{without} expensive additional
computations) which is taken up again in \autoref{522}.
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{fig:a2} in \ref{appendix-a3} 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 \emph{grf}, the \emph{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 \emph{mcf} weights
weigh the observed outcome directly, \emph{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$,
\emph{grf} estimates $N$ weights on approximately $B/4$ trees,
while \emph{mcf} estimates $N/2$ weights on $B$ trees as
illustrated in \ref{appendix-a3} in Figures \ref{fig:a1} and \ref{fig:a2}. 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 \emph{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 \emph{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{appendix-a1} contains the (so-far missing) theoretical guarantees for
the \emph{mcf}. Under similar assumptions as in \citetalias{wager2018estimation}, the \emph{mcf}
estimators of the IATEs, GATEs, and ATEs are consistent and pointwise
asymptotically normal. The \emph{mcf} inference procedures exploit the
fact that \emph{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{appendix-a2}, while
\ref{appendix-a3} explains some of the difference to the \emph{grf} at the
implementation level.
\subsection{Overview}\label{44}
This section summarizes the statistical guarantees of the three CCMLs
for estimation of ATE, GATE, and IATE. For the ATE, \emph{dml} and
\emph{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 \autoref{41} and \autoref{42}}
The \emph{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. \emph{dml} and \emph{grf} differ from the \emph{mcf} in the
utilization of the observed data points. While \emph{dml} and \emph{grf}
compute the average of estimated components of the DR score for
observations in the corresponding group of interest, \emph{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 \autoref{43}. This will affect the efficiency
of the GATE estimators in different scenarios as captured in the
following proposition.
\begin{prop} Let $J$ denote the number of groups for
which GATEs are to be estimated and let $\mathit{\Pi}$ be a $\mathit{J \times 1}$
vector containing the non-zero population shares of observations in the
corresponding groups. Suppose that the total number of groups increases from $J$ to $J'$ and
the population share of group $j$ under $J$ groups, $\pi_j$, decreases to $\pi_{j'}$ under $J'$ groups.
Then, the asymptotic variances of the efficient \emph{dml} and \emph{grf} GATE
estimators increase for each group $j$ within $J$
proportionally to $\pi_j/\pi_{j'}$, while
the variance of \emph{mcf} estimator of the GATE does not depend on the
number of groups and remains stable. This implies $Var\left( {\widehat {GAT{E^{mcf}}}\left( {m,l;{z_j}} \right)} \right) \le Var\left( {\widehat {GAT{E^{dml/grf}}}\left( {m,l;{z_j}} \right)} \right)$
for large $J$.
\label{prop1-main}
\end{prop}
For the IATE, there are no specific statistical guarantees for
\emph{dml}. There is an upper bound on the error of the DR-learner
providing a general guarantee. \emph{grf} and \emph{mcf} guarantee
pointwise consistency and asymptotic normality and are both expected to
yield lower MSEs in the simulations than \emph{dml}-based DR-learner.
For \emph{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{tab2}.
\begin{table}[h]
\centering
\caption{Selected statistical properties of the CCMLs}
\label{tab2}
\begin{tabularx}{\textwidth}{X>{\RaggedRight\arraybackslash}X>{\RaggedRight\arraybackslash}X>{\RaggedRight\arraybackslash}X}
\toprule
\textbf{Approach} & \emph{ATE} & \emph{GATE} & \emph{IATE} \\
\midrule
\textbf{\emph{dml}} &
Consistent, asymptotically normal, asymptotically efficient &
Consistent, asymptotically normal, uniform convergence &
Depends on implementation of the smoothing (regression) step \\
\midrule
\textbf{\emph{grf}} &
Consistent, asymptotically normal, asymptotically efficient &
Consistent, asymptotically normal, uniform convergence &
Consistent, asymptotically normal, pointwise convergence known (slow) convergence rate \\
\midrule
\textbf{\emph{mcf}} &
Consistent, asymptotically normal &
Consistent, asymptotically normal, pointwise convergence &
Consistent, asymptotically normal, pointwise convergence \\
& & & \\
& & More efficient for large number of groups (compared to \emph{dml} and \emph{grf}) & \\
\bottomrule
\end{tabularx}
\end{table}
\section{Monte Carlo study}\label{5-mcs}
\subsection{Concept}\label{51}
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{tab:b2} in
\ref{appendix-b25} gives the list of the scenarios investigated. \ref{appendix-c}
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.
\subsection{Key features of the simulation study}\label{52}
\subsubsection{Data Generating Process}\label{521}
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{appendix-b}).
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{fig1} shows the distribution of the true treatment
probabilities for the strong and medium selectivity case with two
treatments and 50\% treatment share.
\begin{figure}[h!]
\caption{Distribution of treatment probabilities in treated and non-treated population}
\label{fig1}
\includegraphics[width = \textwidth]{ps.png}
\caption*{\footnotesize \textsf{Note: Two treatments and 50\% treatment share (base specification). Figures are based on 1'000'000 observations.}}
\end{figure}
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 \autoref{53}.\footnote{The base
specifications with different degrees of selectivity and types of
IATEs are captured in rows 1-15 in Table \ref{tab:b2} in \ref{appendix-b25}.} The
smooth function, when plotted as a function of the linear index,
exhibits a pattern that resembles a step function. Figure \ref{fig2} 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{appendix-b24}.
\begin{figure}[h!]
\caption{Shape of potential outcomes for different shapes of IATEs (base specification)}
\label{fig2}
\includegraphics[width = \textwidth]{pot_out_main.png}
\caption*{\raggedright \footnotesize \textsf{Note: Figures are based on 1'000'000 observations.}}
\end{figure}
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 \citet{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.
\subsubsection{Estimators}\label{522}
Below, we present the results of the Modified Causal Forest and the
Generalized Random Forest, both with the outcome variable centred prior
to estimation, \emph{mcf-cent} and \emph{grf-cent}, as well as of
Double/debiased Machine Learning with normalized weights,
\emph{dml-norm}. For IATE estimation, we also include a more efficient
version of the centred \emph{mcf} for which no inference is available
(see below), \emph{mcf-cent-eff}. The tables in Appendix C contain
additional results for the uncentred \emph{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{appendix-b3} for details.
While \ref{appendix-b31} details the implementation of the \emph{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 \emph{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 \emph{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 \emph{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 \emph{grf} procedure by switching off the half-sampling
when inference is not relevant. However, the current implementation in
the \emph{grf} package does not support this option.
As mentioned in \autoref{42}, \emph{grf} performs local centring inside
the Random Forest to remove confounding bias. The default \emph{grf}
local centring, labelled here as \emph{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 \emph{centred}. The difference between
the centred and the uncentred \emph{mcf} and \emph{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 \autoref{42}.}
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 \emph{dml}, the normalized \emph{dml} is more
robust to extreme values of the estimated propensity score. It is
obtained by normalizing the weights that are implicit in the \emph{dml}
scores. While it is straightforward to use \emph{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 \emph{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.
\subsubsection{Performance measures}\label{523}
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{appendix-c} 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.
\subsection{Results}\label{53}
In this subsection, we analyse \emph{dml-norm}, \emph{grf-cent},
and \emph{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
\emph{-norm} and \emph{-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.
\subsubsection{Average treatment effects}\label{531}
Table \ref{tab3} shows the results for the ATE. If there is no selectivity, like
in an experiment, we obtain the expected results: the efficient
estimators, \emph{dml} and \emph{grf}, are similar and outperform the
\emph{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
(\emph{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{tab:c15} in \ref{appendix-c1}, this behaviour appears for the not
normalized \emph{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 \emph{grf}, although it is not as extreme as for
\emph{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 \emph{mcf}.
One summary coming from this table is that \emph{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{appendix-c1} 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 (!) \emph{grf}
outperforms the centred \emph{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{appendix-c2}), suggests that the patterns observed in Table \ref{tab3}
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{tab:c17}, Table \ref{tab:c21}): all
estimators become biased, which increases the RMSE and leads to coverage
far below nominal levels.
\begin{table}[h!]
\sffamily
\small
\caption{Simulation results for average treatment effect (ATE)}
\label{tab3}
\resizebox{\textwidth}{!}{
\begin{tabular}{lcccccccc}
\toprule
\multicolumn{3}{c}{} & \multicolumn{4}{c}{\textbf{Estimation of effects}} & \multicolumn{2}{c}{\textbf{Inference}} \\
\midrule
\textbf{Estimator} & Selec- & Sample & Bias & Mean & Std. dev & RMSE& Bias (SE) & CovP \\
& tivity & size & & absolute & & & (SE) & (95) in \\
& & & & error& & & & \% \\
\midrule
\multicolumn{1}{c}{(1)} & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9) \\
\midrule
\textbf{dml-norm} & None & 2'500 & -0.002 & 0.034 & 0.043 & 0.043 & 0.005 & 97 \\
\textbf{grf-cent} & & & 0.003 & 0.033 & 0.042 & 0.042 & 0.000 & 95 \\
\textbf{mcf-cent} & & & 0.003 & 0.048 & 0.061 & 0.061 & -0.001 & 94 \\ \hdashline
\textbf{dml-norm} & Med- & 2'500 & 0.020 & 0.040 & 0.045 & 0.050 & 0.004 & 95 \\
\textbf{grf-cent} & ium & & 0.037 & 0.046 & 0.042 & 0.056 & -0.000 & 85 \\
\textbf{mcf-cent} & & & 0.037 & 0.057 & 0.061 & 0.072 & 0.000 & 91 \\ \hdashline
\textbf{dml-norm} & Strong & 2'500 & 0.140 & 0.140 & 0.056 & 0.150 & 0.002 & 23 \\
\textbf{grf-cent} & & & 0.090 & 0.091 & 0.046 & 0.101 & -0.006 & 41 \\
\textbf{mcf-cent} & & & 0.046 & 0.060 & 0.058 & 0.074 & 0.008 & 92 \\ \midrule
\textbf{dml-norm} & None & 10'000 & 0.000 & 0.016 & 0.020 & 0.020 & 0.003 & 98 \\
\textbf{grf-cent} & & & 0.002 & 0.017 & 0.021 & 0.021 & 0.000 & 95 \\
\textbf{mcf-cent} & & & 0.004 & 0.022 & 0.028 & 0.028 & 0.002 & 97 \\ \hdashline
\textbf{dml-norm} & Med- & 10'000 & 0.010 & 0.019 & 0.022 & 0.024 & 0.003 & 94 \\
\textbf{grf-cent} & ium & & 0.015 & 0.022 & 0.022 & 0.027 & -0.001 & 92 \\
\textbf{mcf-cent} & & & 0.014 & 0.025 & 0.027 & 0.030 & 0.003 & 96 \\ \hdashline
\textbf{dml-norm} & Strong & 10'000 & 0.079 & 0.079 & 0.032 & 0.085 & 0.000 & 28 \\
\textbf{grf-cent} & & & 0.042 & 0.044 & 0.027 & 0.050 & -0.006 & 46 \\
\textbf{mcf-cent} & & & -0.013 & 0.025 & 0.028 & 0.031 & 0.007 & 96 \\ \bottomrule
\end{tabular}
}
\vspace{-8pt}
\captionsetup{justification=justified,singlelinecheck=false}
\caption*{\footnotesize \textsf{Note: RMSE abbreviates the Root Mean Squared Error. \emph{CovP (95)} denotes the (average) probability that the true value is part of the estimated 95\% confidence interval. 1'000 / 250 replications are used for 2'500~/ 10'000 observations.}}
\end{table}
\subsubsection{Conditional average treatment effects with a small number of groups (GATEs)}\label{532}
Table \ref{tab4} 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{prop1-main}.
Table \ref{tab4} shows the results for the cases of 5 and 40
groups, while Tables \ref{tab:c28} to \ref{tab:c30} in \ref{appendix-c2} show the intermediate
cases of 10 and 20 groups as well.
\begin{table}[h!]
\sffamily
\scriptsize
\caption{Simulation results for group average treatment effects (GATEs)}
\label{tab4}
\resizebox{\textwidth}{!}{
\begin{tabular}{lcccccccc}
\toprule
\multicolumn{3}{c}{} & \multicolumn{4}{c}{\textbf{Estimation of effects}} & \multicolumn{2}{c}{\textbf{Inference}} \\
\midrule
\textbf{Estimator} & Selec- & Sample & Bias & Mean & Std. dev & RMSE& Bias (SE) & CovP \\
& tivity & size & & absolute & & & (SE) & (95) in \\
& & & & error& & & & \% \\
\midrule
\multicolumn{1}{c}{(1)} & (2) & (3) & (4) & (5) & (6) & (7) & (8) & (9) \\
\midrule
\multicolumn{9}{c}{5 Groups} \\
\midrule
\textbf{dml-norm} & None & 2'500 & -0.002 & 0.075 & 0.094 & 0.094 & 0.002 & 95 \\
\textbf{grf-cent} & & & 0.004 & 0.075 & 0.094 & 0.094 & 0.000 & 95 \\
\textbf{mcf-cent} & & & 0.002 & 0.090 & 0.091 & 0.114 & -0.010 & 83 \\ \hdashline
\textbf{dml-norm} & Med- & 2'500 & 0.020 & 0.081 & 0.081 & 0.101 & 0.001 & 95 \\
\textbf{grf-cent} & ium & & 0.037 & 0.081 & 0.093 & 0.102 & -0.000 & 93 \\
\textbf{mcf-cent} & & & 0.037 & 0.094 & 0.090 & 0.121 & -0.006 & 84 \\ \hdashline
\textbf{dml-norm} & Strong & 2'500 & 0.140 & 0.157 & 0.123 & 0.187 & -0.005 & 74 \\
\textbf{grf-cent} & & & 0.090 & 0.119 & 0.096 & 0.148 & -0.007 & 75 \\
\textbf{mcf-cent} & & & 0.046 & 0.112 & 0.085 & 0.146 & 0.009 & 81 \\ \midrule
\textbf{dml-norm} & None & 10'000 & -0.001 & 0.037 & 0.046 & 0.046 & 0.001 & 95 \\
\textbf{grf-cent} & & & 0.002 & 0.036 & 0.046 & 0.046 & 0.001 & 95 \\
\textbf{mcf-cent} & & & 0.004 & 0.045 & 0.045 & 0.058 & -0.002 & 87 \\ \hdashline
\textbf{dml-norm} & Med- & 10'000 & 0.010 & 0.041 & 0.050 & 0.051 & 0.000 & 94 \\
\textbf{grf-cent} & ium & & 0.016 & 0.041 & 0.048 & 0.051 & -0.001 & 92 \\
\textbf{mcf-cent} & & & 0.013 & 0.049 & 0.045 & 0.063 & 0.004 & 86 \\ \hdashline
\textbf{dml-norm} & Strong & 10'000 & 0.078 & 0.090 & 0.073 & 0.108 & -0.007 & 71 \\
\textbf{grf-cent} & & & 0.043 & 0.063 & 0.051 & 0.080 & -0.006 & 75 \\
\textbf{mcf-cent} & & & -0.014 & 0.085 & 0.045 & 0.100 & 0.007 & 65 \\ \midrule
\multicolumn{9}{c}{40 Groups} \\
\midrule
\textbf{dml-norm} & None & 2'500 & -0.002 & 0.214 & 0.268 & 0.269 & -0.003 & 94 \\
\textbf{grf-cent} & & & 0.004 & 0.213 & 0.269 & 0.269 & -0.001 & 95 \\
\textbf{mcf-cent} & & & 0.002 & 0.096 & 0.096 & 0.120 & -0.013 & 82 \\ \hdashline
\textbf{dml-norm} & Med- & 2'500 & 0.020 & 0.226 & 0.283 & 0.284 & -0.005 & 94 \\
\textbf{grf-cent} & ium & & 0.037 & 0.213 & 0.265 & 0.268 & -0.001 & 95 \\
\textbf{mcf-cent} & & & 0.036 & 0.099 & 0.095 & 0.128 & -0.008 & 83 \\ \hdashline
\textbf{dml-norm} & Strong & 2'500 & 0.140 & 0.293 & 0.338 & 0.367 & -0.020 & 89 \\
\textbf{grf-cent} & & & 0.090 & 0.223 & 0.255 & 0.280 & -0.003 & 92 \\
\textbf{mcf-cent} & & & 0.043 & 0.119 & 0.090 & 0.153 & 0.007 & 80 \\ \midrule
\textbf{dml-norm} & None & 10'000 & -0.001 & 0.106 & 0.132 & 0.132 & 0.000 & 95 \\
\textbf{grf-cent} & & & 0.002 & 0.104 & 0.130 & 0.130 & 0.002 & 95 \\
\textbf{mcf-cent} & & & 0.004 & 0.050 & 0.050 & 0.063 & -0.005 & 85 \\ \hdashline
\textbf{dml-norm} & Med- & 10'000 & 0.010 & 0.113 & 0.140 & 0.141 & -0.001 & 95 \\
\textbf{grf-cent} & ium & & 0.016 & 0.106 & 0.131 & 0.132 & 0.000 & 95 \\
\textbf{mcf-cent} & & & 0.012 & 0.055 & 0.051 & 0.069 & -0.002 & 84 \\ \hdashline
\textbf{dml-norm} & Strong & 10'000 & 0.079 & 0.164 & 0.192 & 0.209 & -0.015 & 89 \\
\textbf{grf-cent} & & & 0.043 & 0.115 & 0.130 & 0.145 & -0.002 & 92 \\
\textbf{mcf-cent} & & & -0.016 & 0.091 & 0.050 & 0.104 & 0.005 & 63 \\ \bottomrule
\end{tabular}
}
\vspace{-8pt}
\captionsetup{justification=justified,singlelinecheck=false}
\caption*{\footnotesize \textsf{Note: RMSE abbreviates the Root Mean Squared Error. \emph{CovP (95)} denotes the (average) probability that the true value is part of the estimated 95\% confidence interval. 1'000 / 250 replications are used for 2'500~/ 10'000 observations.}}
\end{table}
Concerning the point estimate, \emph{mcf} outperforms \emph{grf} and
\emph{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 \emph{mcf} performs better for strong
selectivity: the previous section showed that for strong selectivity the
\emph{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 \emph{mcf} than
for \emph{dml} and \emph{grf}. The reason is the way the GATE estimators
are constructed by the different methods. As mentioned above, \emph{dml}
and \emph{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 \emph{grf} and \emph{dml} estimators is exactly what is observed
in Tables \ref{tab:c28} to \ref{tab:c30} when comparing results for 5 vs. 20 groups, or 10
vs. 40 groups, respectively.
The \emph{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 \emph{mcf} increases much slower
than for \emph{grf} and \emph{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. \citet{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 \emph{mcf} estimator
across varying group sizes.} In fact, we observe that when the number
of groups increases further, the RMSE of the \emph{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: \emph{dml}
and \emph{grf} have the correct coverage rates for no and medium
selectivity even for 40 groups, while the coverage rates for \emph{mcf}
are too low.\footnote{Note that when selection into treatment gets
stronger, splitting for the \emph{mcf} (and \emph{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{tab:c28},
\ref{tab:c29}, and \ref{tab:c30} in \ref{appendix-c2} 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, \emph{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
\subsubsection{Conditional average treatment effects at the finest aggregation level (IATEs)}\label{533}
Table \ref{tab5} 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 \emph{mcf} the more efficient version is included as
well.
\begin{table}[h!]
\sffamily
\scriptsize
\caption{Simulation results for individualized average treatment effects (IATEs)}
\label{tab5}
\resizebox{\textwidth}{!}{
\begin{tabular}{lcccccc}
\toprule
\multicolumn{3}{c}{} & \multicolumn{3}{c}{\textbf{Estimation of effects}} & \multicolumn{1}{c}{\textbf{Inference}} \\
\midrule
\textbf{Estimator} & Selectivity & Sample size & Bias & Mean abso- & RMSE & CovP (95) in \\
& & & & lute error & & \% \\
\midrule
\multicolumn{1}{c}{(1)} & (2) & (3) & (4) & (5) & (7) & (9) \\
\midrule
\textbf{dml-ols} & None & 2'500 & -0.002 & 0.229 & 0.286 & 23 \\
\textbf{dml-rf} & & & -0.002 & 0.270 & 0.342 & - \\
\textbf{grf-cent} & & & 0.004 & 0.155 & 0.187 & 77 \\
\textbf{mcf-cent} & & & 0.002 & 0.162 & 0.199 & 77 \\
\textbf{mcf-cent-eff} & & & 0.004 & 0.151 & 0.184 & - \\ \hdashline
\textbf{dml-ols} & Med- & 2'500 & 0.020 & 0.236 & 0.295 & 24 \\
\textbf{dml-rf} & ium & & 0.013 & 0.284 & 0.362 & - \\
\textbf{grf-cent} & & & 0.034 & 0.175 & 0.210 & 72 \\
\textbf{mcf-cent} & & & 0.037 & 0.174 & 0.212 & 75 \\
\textbf{mcf-cent-eff} & & & 0.040 & 0.163 & 0.197 & - \\ \hdashline
\textbf{dml-ols} & Strong & 2'500 & 0.140 & 0.285 & 0.357 & 28 \\
\textbf{dml-rf} & & & 0.117 & 0.349 & 0.461 & - \\
\textbf{grf-cent} & & & 0.092 & 0.254 & 0.297 & 55 \\
\textbf{mcf-cent} & & & 0.046 & 0.212 & 0.252 & 68 \\
\textbf{mcf-cent-eff} & & & 0.048 & 0.205 & 0.240 & - \\ \midrule
\textbf{dml-ols} & None & 10'000 & -0.001 & 0.179 & 0.219 & 7 \\
\textbf{dml-rf} & & & 0.000 & 0.201 & 0.256 & - \\
\textbf{grf-cent} & & & 0.002 & 0.082 & 0.103 & 91 \\
\textbf{mcf-cent} & & & 0.004 & 0.096 & 0.121 & 83 \\
\textbf{mcf-cent-eff} & & & 0.003 & 0.088 & 0.110 & - \\ \hdashline
\textbf{dml-ols} & Med- & 10'000 & 0.010 & 0.181 & 0.222 & 8 \\
\textbf{dml-rf} & ium & & 0.005 & 0.214 & 0.274 & - \\
\textbf{grf-cent} & & & 0.016 & 0.089 & 0.111 & 89 \\
\textbf{mcf-cent} & & & 0.013 & 0.106 & 0.131 & 81 \\
\textbf{mcf-cent-eff} & & & 0.013 & 0.100 & 0.122 & - \\ \hdashline
\textbf{dml-ols} & Strong & 10'000 & 0.079 & 0.205 & 0.256 & 12 \\
\textbf{dml-rf} & & & 0.058 & 0.275 & 0.418 & - \\
\textbf{grf-cent} & & & 0.053 & 0.123 & 0.162 & 82 \\
\textbf{mcf-cent} & & & -0.014 & 0.156 & 0.185 & 66 \\
\textbf{mcf-cent-eff} & & & -0.014 & 0.151 & 0.151 & - \\ \bottomrule
\end{tabular}
}
\vspace{-8pt}
\captionsetup{justification=justified,singlelinecheck=false}
\caption*{\footnotesize \textsf{Note: \emph{dml} is the normalized \emph{dml} (denoted as \emph{dml-norm} in Tables \ref{tab3} and \ref{tab4}). RMSE abbreviates the Root Mean Squared Error. \emph{CovP (95)} denotes the (average) probability that the true value is part of the estimated 95\% confidence interval. 1'000 / 250 replications are used for 2'500 / 10'000 observations.}}
\end{table}
For the point estimates the ordering is clear-cut. \emph{grf} and
\emph{mcf} are similar and outperform both \emph{dml} versions
(\emph{ols}, \emph{rf}). In many cases, the efficient centred \emph{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 \emph{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.
\subsubsection{Summary}\label{summary}
\emph{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
\emph{dml}-based estimator in the strong selectivity setting could be
the automatic \emph{dml} that would estimate inverse probability
weights instead of propensity scores as nuisance parameters.}
\emph{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 \emph{dml}-based and
\emph{grf} -based GATE estimators to further improve their RMSE.}
However, it performs well for IATEs. It is also less affected by strong
selectivity than \emph{dml}. These results however only hold for the
modified version of the \emph{grf} in which the outcome variable in the
training data is explicitly centred before it enters the \emph{grf}
algorithm. The original uncentred version of \emph{grf} as suggested in
\citet{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 \autoref{522}, 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
\autoref{42}.}
Compared to \emph{dml} and \emph{grf}, the \emph{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 \emph{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 \emph{dml} and \emph{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.
\section{Conclusion}\label{6-conclusion}
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 (\emph{dml}), the
Generalized Random Forest (\emph{grf}), and the Modified Causal Forest
(\emph{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 \emph{dml} and the \emph{grf}, they are already known,
but not so for the \emph{mcf}. Therefore, we explicitly provide them.
The large-scale simulation study reveals scenarios in which the methods
perform well. \emph{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. \emph{grf} shows similar behaviour as
\emph{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, \emph{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 \emph{mcf} is recommended to estimate IATEs. For
large groups, \emph{mcf} is recommended for point estimation of GATEs.
When sample is large enough and slight loss of efficiency is not
detrimental, \emph{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 \emph{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: \emph{dml} is available as Python
and R packages, the g\emph{rf} is available as an R package, and the
\emph{mcf} is available as a Python package.\footnote{As examples can
serve: \citet{bach2022doubleml}, \citet{bach2024dmlr}, and \citet{knaus2022double}
for \emph{dml}; \citet{athey2019estimating} for \emph{grf}; and
\citet{bodory2022high} for \emph{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 \emph{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 \emph{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.
\bibliography{ref.bib}
\bibliographystyle{apacite}
\newpage
\begin{spacing}{1.5}