\section{Introduction}\label{sec:intro}

Estimating the effect of a treatment on an outcome is a fundamental problem in many fields such as medicine~\citep{shalit-icml17,laber-bio15,lakkaraju-aistats17}, public policy~\citep{grimmer-pa17} and more~\citep{li-www10,ascarza-jmr18}.
For example, doctors might be interested in how a treatment, such as a drug, affects the recovery of patients~\citep{foster-statmed11}, economists may be interested in how a job training program affects employment prospectives~\citep{lalonde-aer86}, and advertisers may want to model the average effect an advertisement has on sales~\cite{lewis-retail08}.

However, individuals may react differently to the treatment of interest, and knowing only the average treatment effect in the population is insufficient.
For example, a drug may have adverse effects on some individuals but not others~\citep{shalit-icml17}, or a person's education and background may affect how much they benefit from job training~\citep{lalonde-aer86,pearl-smr17}.
Measuring the extent to which different individuals react differently to treatment is known as \emph{heterogeneous treatment effect (HTE) estimation}.

Traditionally, HTE estimation has been done through subgroup analysis~\citep{gail-bio95,bonetti-bio04}.
However, this can lead to cherry-picking since the practitioner is the one who identifies subgroups for estimating effects.
Recently,
there has been more focus on \emph{data-driven} estimation of heterogeneous treatment effects by letting the data identify which features are important for treatment effect estimation  using machine learning techniques~\citep{johansson-icml16,louizos-neurips17,shalit-icml17,tran-aaai19}.
A straightforward approach is to create interaction terms between all covariates and use them in a regression~\citep{athey-handbook17}.
However, this approach is limited to low-dimensional datasets and is not practical for large datasets with many features.
To deal with the high dimensionality of data, many supervised learning-based methods have been developed and adapted to the problem of HTE estimation~\citep{su-jmlr09,hill-jcgs11,imai-annals13,athey-stat15,johansson-icml16,athey-pnas16,shalit-icml17,tran-aaai19,kunzel-pnas19}.

HTE estimation methods can be broadly categorized into counterfactual prediction and effect estimation methods.
A counterfactual prediction model learns to predict each unit's potential outcomes when treated and not treated.
In turn, those predictions can be used to estimate a treatment effect~\citep{foster-statmed11,shalit-icml17,kunzel-pnas19}.
A heterogeneous effect estimation method learns a model that predicts the treatment effect directly~\citep{su-jmlr09,athey-pnas16,tran-aaai19}.
Effect estimation models can be learned through specialized objective functions~\citep{su-jmlr09,athey-pnas16} or in tandem with counterfactual prediction~\citep{louizos-neurips17,kunzel-pnas19}.
Typically, counterfactual prediction and effect estimation methods use all variables in the data for estimation relying on the strong ignorability assumption~\citep{montgomery-ajps18}.
However, as we discuss in Section~\ref{sec:structural_causal_models} and show in our experiments, including all variables can lead to errors in HTE estimation.

A third category of methods, which rely on structural causal models (SCMs)~\citep{pearl-book09}, provide criteria for the selection of variables that allow for unbiased causal effects. SCMs are graphical models that model explicitly the cause-effect relationships between variables.
The selection of variables through SCMs can help data-driven HTE estimation methods better estimate causal effects, as we show in our experiments.
However, the underlying causal mechanism of data and its SCM are rarely known in real-world settings. To circumvent this problem, one can use a causal structure learning algorithm~\citep{spirtes-book00} to learn a skeleton of the underlying causal model from the data itself.
This survey provides an overview of recently developed data-driven HTE estimation methods, and an overview of SCMs and how they can be used for HTE estimation.

The rest of the paper is organized as follows.
Section~\ref{sec:background} gives basic background knowledge.
Section~\ref{sec:survey} reviews counterfactual prediction and effect estimation methods.
Section~\ref{sec:structural_causal_models} provides an overview of data-driven HTE estimation with structural causal models.
Section~\ref{sec:experiments} empirically evaluates HTE estimators on several types of datasets.
Section~\ref{sec:conclusion} concludes the paper and discusses open problems.
\section{Heterogeneous Treatment Effects}\label{sec:background}

In this section, we provide background on heterogeneous treatment effect estimation. We define the problem of heterogeneous treatment effect estimation and common assumptions to allow for estimation.
We introduce structural causal models~\citep{pearl-book09} and how they can be used to assess heterogeneity.

\subsection{Data model}\label{sec:data_model}

Before formalizing the problem and data model, we introduce the notation used in this survey.
We denote a space of random variables with a calligraphic letter (e.g., $ \mathcal{X} $).
We denote random variables with a capital letter (e.g., $Y$) and random vectors with a boldfaced capital letter (e.g., $\mathbf{X}$).
Subscripts signify the random variable associated with a specific unit (e.g., $Y_i$ or $\mathbf{X}_i$).
Lowercase letters denote an assignment of a random variable or vector (e.g., $Y_i = y$ or $\mathbf{X}_i = \mathbf{x}$).
Each unit $i$ has a vector of ``pre-treatment'' features (or covariates), $\mathbb{X}_i \in \mathcal{X}$, a binary treatment indicator, $T_i$, where $T_i=0$ and $T_i=1$ indicates $i$ is not treated and treated, respectively, and two outcomes when not treated and treated, $Y_i(T_i=0)$ and $Y_i(T_i=1)$, respectively.
For brevity, we write $Y_i(T_i=0)$ and $Y_i(T_i=1)$ as $Y_i(0)$ and $Y_i(1)$, respectively.
In this survey, we consider the treatment to be binary, as most work in estimating heterogeneous treatment effects only considers binary treatment, some recent work has considered how to estimate heterogeneous effects when the treatment is non-binary~\citep{tran-aaai19,schwab-arxiv19,kennedy-jrss17}
The underlying data distribution is represented by a tuple $(\mathbf{X}_i, T_i, Y_i(0), Y_i(1))^N$ of $N$ independent and identically distributed individuals.
The fundamental problem of causal inference is that for each individual, we can observe either the potential outcomes under treatment ($Y_i(1)$) or under control ($Y_i(0)$), but not both.
Therefore, the observed dataset as $\mathcal{D} = (\mathbf{x}_i, t_i, y_i)^N$ with $N$ individuals where $Y_i$ is defined as:
\begin{equation}
    \label{eq:potential_outcomes}
    y_i = t_i Y_i(1) + (1 - t_i) Y_i(0).
\end{equation}
We denote $ Y_i $ as the factual outcome and the counterfactual outcome, $ Y_i' $, as the outcome that would have occurred if the treatment was different.
We define the treatment group as all individual units $i$ such that $t_i = 1$ and the control group as all units $i$ such that $t_i = 0$.

\subsection{Heterogeneous Treatment Effect Estimation}\label{sec:hte}

The average treatment effect (ATE) over the population is defined as the expected difference in potential outcomes:
\begin{equation}
    \label{eq:ate}
    \text{ATE} := E[Y(1) - Y(0)].
\end{equation}
Since each unit reacts differently to treatment, we define the individual treatment effect (ITE) as the difference in the unit's potential outcomes, denoted by $\tau_i$:
\begin{equation}
    \label{eq:ite}
    \tau_i := Y_i(1) - Y_i(0).
\end{equation}
However, due to the fundamental problem of causal inference, $\tau_i$ can not be observed for any individual unit.
Instead, HTE estimation focuses on estimating the conditional average treatment effect (CATE):
\begin{equation}\label{eq:cate}
    \tau(\mathbf{x}) := \mathbb{E} \left[ Y_i(1) - Y_i(0) \mid \mathbf{X}_i = \mathbf{x} \right].
\end{equation}
The CATE describes the expected treatment effect for units in a subpopulation with features $\mathbf{x}$.
An HTE estimator, defined as $\hat{{\tau}}(\mathbf{X})$, seeks to estimate the CATE defined in equation~\eqref{eq:cate}.
The best estimator for CATE is also the best estimator for ITE in terms of the mean squared error~\citep{kunzel-pnas19}.
The goal of HTE estimation is to develop estimators of CATE which have the smallest root mean squared error for the ITE\@.
\begin{problem}\label{prob:hte_problem}
(Heterogeneous treatment effect estimation)
Given a dataset $\mathcal{D}=(\mathbf{X}_i, Y_i, T_i)$ of N instances with some true ITE, $\vect{\tau} = (\tau_i)^N$, estimate CATE $\hat{\tau}(\mathbf{X})$ such that the root-mean-squared error (RMSE) is minimized:
\begin{equation}\label{eq:cate_objective}
    \minl_{\hat{\tau}} \text{RMSE} = \frac{1}{N} \sum_{i=1}^N  (\tau_i - \hat{\tau}(\mathbf{X}_i)) ^ 2.
\end{equation}
\end{problem}
In practice, equation~\eqref{eq:cate_objective} cannot be estimated in real-world data.
HTE estimators are often evaluated on synthetic or semi-synthetic datasets~\citep{curth-neurips21}, but this may not capture performance in real-world scenarios.
Some work has been done on evaluating HTE estimators heuristically, \citep{schuler-model18,saito-cv19,schwab-perfect18}, which we describe in Section~\ref{sec:data_driven_background}.

\subsection{HTE Assumptions}\label{sec:hte_assumptions}

To estimate HTEs from observational data, some assumptions on the set of features and data given are needed.
We present some common assumptions used in heterogeneous treatment effect estimators.

\begin{assumption}\label{assumption:strong_ignorability}
    (Strong Ignorability) Given a binary treatment $T$, potential outcomes $Y(0)$, $Y(1)$, and a vector of features $\mathbf{X}$, strong ignorability states:
    \begin{align}
        \{ Y(0), Y(1) \} \perp T \mid \mathbf{X},
    \end{align}
\end{assumption}
This states the \emph{potential} outcomes are conditionally independent of the treatment conditioned on the features.
When strong ignorability holds, $\mathbf{X}$ is said to be ``admissible'', meaning treatment effects can be estimated without bias.
This assumption is also known as unconfoundedness or the assumption of no hidden confounding.

\begin{assumption}\label{assumption:overlap}
    (Overlap) For every unit $i$, $0 < P(T_i = 1) < 1$ with probability 1.
\end{assumption}
This assumption states that all units have a non-zero probability of assignment to each treatment condition.
Assumption~\ref{assumption:overlap} is also known as positivity or common support.

Another common assumption used is the Stable Unit Treatment Value Assumption (SUTVA) or the assumption of no interference~\citep{rubin-jasa80}.
\begin{assumption}\label{assumption:sutva}
    (Stable Unit Treatment Value Assumption (SUTVA)) Given any two units $i$ and $j$, the potential outcomes of $i$, $\{Y_i(0), Y_i(1) \}$ are unaffected by the treatment assignment of unit $j$:
    \begin{equation}
        \{Y_i(0), Y_i(1)\} \perp T_j, \forall j \neq i.
    \end{equation}

\end{assumption}
With the SUTVA assumption, individual outcomes are not affected by the treatment or outcomes of peers.
Though there has been some work focused on estimating HTEs without the SUTVA assumption~\citep{yuan-www21,tran-aaai22}, a majority of work in HTE estimation has focused on non-networked settings, which we focus on in this survey.
\subsection{Counterfactual prediction}\label{sec:counterfactuals}

The first category of methods predicts counterfactuals from data, the outcome when treated and control, $ \hat{Y}(1), \hat{Y}(0) $, respectively, and uses these predictions for effect estimation, $ \hat{Y}(1)-\hat{Y}(0)$.

\subsubsection{Single-model approach}

The simplest way to make counterfactual predictions is by using the \emph{single-model approach}.
The single-model approach estimates CATE using a single regression or classification model.
Given the data \( D = (\mathbf{X}_i, Y_i, T_i) \), a supervised method is trained to predict \( Y_i \) given \( \mathbf{X}_i, T_i \), \( \hat{\mu}(T_i, \mathbf{X}_i) = \hat{Y}_i \).
The CATE can be predicted as the difference between the prediction with \( T_i = 1 \) and \( T_i = 0 \).
\begin{equation}
    \hat{\tau}(\mathbf{x}) = \hat{\mu}( T = 1, \mathbf{X} = x) - \hat{\mu}(T = 0, \mathbf{X} = x).
\end{equation}
Hill used Bayesian Additive Regression Trees (BART) to predict treatment effects~\citep{hill-jcgs11}.
Athey and Imbens used regression trees as the single-model approach~\citep{athey-stat15}.
K{\"{u}}nzel et al.\ used a single-model approach with gradient-boosted trees and random forests~\citep{kunzel-pnas19} and name these single-model approaches S-Learner.

Since the single-model approach only requires regression or classification, any off-the-shelf supervised learning algorithm can be used.
However, the single-model approach may not be able to capture potential differences in treated and control groups and bias estimation.
Additionally, in some cases, it may not be possible for the treatment to be used in prediction (e.g., not used as a split in a regression tree).
This results in the estimator predicting the same outcome for different levels of treatment, which results in a zero for CATE estimation.

Besides methods using classical machine learning methods, several deep learning approaches for estimating HTEs have been proposed recently using the single-model approach.
Atan et al.\ propose a deep learning model that consists of two stages~\citep{atan-aaai18}.
In the first stage, Deep-Treat uses a bias-removing auto-encoder to learn a balanced representation, $ \Phi(\mathbf{x}) $, of treatment and control groups.
The balancing is learned through cross-entropy loss between the marginal and conditional treatment distributions given the hidden representation.
Deep-Treat trains a neural network to predict the outcome using the learned representation and the treatment in the second state.

\subsubsection{Two-model approaches}

In contrast to the single-model approach, another way to model heterogeneity is to model two potential outcomes separately, called the \emph{two-model approach}.
In the two-model approach, a separate model is built on treated and control groups.
Then, to estimate CATE, the difference between treatment group model prediction and control group model prediction is used to estimate the effect.

Formally, let \( \hat{\tau}_1(\mathbf{x}) \) and \( \hat{\tau}_2(\mathbf{x}) \) be supervised estimators for the outcome for the treated and control groups, respectively.
The CATE is estimated as:
\begin{align}
    \hat{\tau}(\mathbf{x}) & = E[Y \mid T = 1, \mathbf{X} = \mathbf{x}] - E[Y \mid T = 0, \mathbf{X} = \mathbf{x}] \\
    & = \hat{\tau}_1(\mathbf{x}) - \hat{\tau}_2(\mathbf{x}).
\end{align}
The main difference between the two-model and single-model approaches is that the two-model approach trains two separate estimators for treatment and control outcome prediction instead of one model for both.
The two model approach has been used with Linear Regression~\citep{cai-bio11}, Trees~\citep{athey-stat15,hill-jcgs11}, Random Forests~\citep{foster-statmed11}, and Gradient Boosted Trees~\citep{kunzel-pnas19}.
Figure~\refeq{fig:slearner_tlearner} contrasts how single-model and two-model approaches perform training and validation.

\begin{figure}[h]
    \centering
    \begin{subfigure}[b]{0.45\linewidth}
        \centering
        \includegraphics[width=\textwidth]{figures/ppt/SLearner.png}
        \caption{Single-model approach.}\label{fig:slearner}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.45\linewidth}
        \centering
        \includegraphics[width=\textwidth]{figures/ppt/TLearner.png}
        \caption{Two-model approach.}\label{fig:tlearner}
    \end{subfigure}
    \caption{A comparison of training and prediction of effects using single-model and two-model approaches.
    The single-model approach uses one estimator for training, and for prediction, features are appended with a treatment indicator of 1 and a treatment indicator of 0 for two separate predictions based on the single model. Then the effect is estimated as the difference.
    In the two-model approach, training is done using two separate models for treated and control groups.
    For prediction, the features are input into the two separate models, and a difference is computed as the predicted effect.}\label{fig:slearner_tlearner}
\end{figure}

In addition, recent work has applied the two-model approach with deep learning estimators.
Johansson et al.\ proposed Balancing Neural Network (BNN), which aims to leverage representation learning and balanced representations for treatment and control groups~\citep{johansson-icml16}.
They learn a neural network architecture with two connecting parts to predict the outcome.
The first part learns the representation of the features, $ \Phi(\mathbf{x}) $.
In the second part, they append the treatment to the hidden representation and learn a representation to predict the outcome.
Separately, they constrain $ \Phi(\mathbf{x}) $ to be similar between treated, $ \Phi_{t=1}(\mathbf{x}) $, and control, $ \Phi_{t=0}(\mathbf{x})$.
The intuition is that the hidden representations should be ``balanced''.
The basic BNN architecture is shown in Figure~\ref{fig:bnn}.

Shalit et al.\ extend the idea of BNNs to the Treatment-Agnostic Representation Network (TARNet) and Counterfactual Regression (CFR)~\citep{shalit-icml17}.
They extend the discrepancy metric by introducing a more flexible family of algorithms in Integral Probability Metrics (IPM) for balancing distributions.
Instead of considering a one-headed network that predicts an outcome based on appending the treatment to the hidden representation, they consider two heads predicting outcome when treated and not treated separately.
The intuition is that the treatment indicator appended to $ \Phi(\mathbf{x}) $ may lose importance in the second part of the network, and so explicitly separating the final predictions can capture differences in treatment and control.

Liu et al.\ propose Split Covariate Representation Network (SCRNet)~\citep{liu-acml20}.
In this work, they propose to divide the covariate space into four categories: confounders, instrumental variables, adjustment variables, and irrelevant variables.
Their goal is to identify the category of each feature and use only relevant features for training a neural network.
The intuition is that irrelevant features are unnecessary for estimation and instrumental variables may bias effect estimations~\citep{pearl-arxiv12,myers-aje11,middleton-pa16}.
Instead, they identify confounders, which confound or are parents of the treatment in outcome in an SCM, and parents of outcome only which are referred to as adjustment variables.
Then, they train a neural network where confounders and adjustment variables are input separately and have their own hidden representations before being concatenated and separated into a two-head network, similar TARNet.
The main difference is the splitting of the initial input into confounding and adjustment variables.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.7\textwidth]{figures/ppt/bnn.png}
    \caption{A high-level architecture of the Balancing Neural Network (BNN) from~\citep{johansson-icml16}. The first part learns a representation of the features only, $ \Phi $.
    In the second part, the treatment is appended, and the neural network is optimized for prediction performance and minimizing the representations of treated and control representations.}\label{fig:bnn}
\end{figure}

\subsubsection{Multi-Model Approaches}

Related to the TARNet approach~\citep{shalit-icml17}, Schwab et al.\ propose Perfect Match~\citep{schwab-arxiv19}.
They extend TARNet in two ways.
First, they propose using $ k $ heads to predict multiple treatment outcomes.
Second, they propose a new way to optimize the neural network model based on nearest neighbor matching and the precision in estimation of heterogeneous effect (PEHE)~\citep{hill-jcgs11}, which measures the mean-squared error in HTE estimation, called NN-PEHE\@.
They propose computing an imputed treatment effect based on nearest neighbors.
Since the true causal effect is unknown for individual units, they substitute the true counterfactual outcome with the outcome from the nearest neighbor.
Their goal is to minimize NN-PEHE for training Perfect Match.
Shi et al.\ proposed another extension of TARNet, called Dragonnet~\citep{shi-neurips19}.
Instead of a two-head neural network, they introduce an additional head which predicts the propensity score from the hidden layer.
In addition, they introduce a modification of the objective function for training neural networks called targeted regularization.

\subsection{Effect estimation}\label{sec:effect_estimation}

In contrast to methods that perform counterfactual prediction, effect estimation methods estimate the treatment effect directly.
They include methods that utilize supervised learning~\citep{athey-stat15,kunzel-pnas19,yoon-iclr18} and methods that have tailored objectives for learning HTEs~\citep{athey-pnas16,tran-aaai19,wager-jasa18}.

\subsubsection{Transformed outcomes}

In the transformed outcomes approach, the outcome variable is transformed such that a regression or classification method can be used to predict the transformed outcome.
Any off-the-shelf estimator can be directly applied since the CATE is estimated directly from the transformed outcome.

Under the assumption of unconfoundedness, the transformed outcome is~\citep{athey-stat15}:
\begin{equation}
    Y_i^* = Y_i \cdot \frac{T_i - e(\mathbf{X}_i)}{e(\mathbf{X}_i) \cdot (1 - e(\mathbf{X}_i) )},
\end{equation}
where \( Y_i^* \) is the transformed outcome and \( e(\mathbf{x}) = \Pr(T_i = 1 \mid \mathbf{X}_i = \mathbf{x}) \) is the propensity score~\citep{rosenbaum-bio83}.
In practice, the propensity score is estimated from data using a supervised estimator such as logistic regression~\citep{setoguchi-pharm08}.
In this setting, the CATE is estimated as the estimated conditional mean of the transformed outcome.
Using the transformed outcome, the conditional mean of the predicted outcome is equal to the CATE:
\begin{equation}
    \hat{\tau}_{to}(\mathbf{x}) = \mathbb{E}[Y^{*} \mid \mathbf{X} = \mathbf{x}].
\end{equation}
Transformed outcomes has been used with regression trees~\citep{athey-stat15} and for building transformed outcome forests~\citep{powers-statmed18}.

\subsubsection{Meta-learning}

Meta-learners or meta-algorithms build upon supervised learning methods, called base learners, to estimate CATE\@.
Two special cases of meta-algorithms are the single-model and two-model approaches, as they take supervised base learners and utilize their outputs to estimate CATE\@.
Xie et al.\ propose three algorithms based on propensity scores: Stratification-Multilevel (SM), Matching-Smoothing (MS), and Smoothing-Differencing (SD)~\citep{xie-soc12}.
The common trend of these algorithms is to estimate propensity scores using a supervised learning method and create strata based on propensity scores and covariates.
In SM, stratum-specific treatment effects are estimated, and a trend across strata is obtained through variance-weighted least squares regression.
For MS, treated units are matched to control units with a matching algorithm to estimate the true effect, and a nonparametric model is applied to the matched effects.
Finally, for SD, nonparametric regression is used to fit the outcome on the propensity score.
To estimate heterogeneous effects, take the difference between the nonparametric regression line between treated and untreated at different levels of the propensity score.

Powers et al.\ propose Pollinated Transformed Outcome (PTO) forests~\citep{powers-statmed18}.
First, they use a transformed outcome approach and build a random forest.
Then, predictions in trees are replaced with the estimate of $ \hat{\mu}(T=1, \mathbf{x}) $ and $ \hat{\mu}(T=0, \mathbf{x}) $ as in the single model approach.
Finally, an additional random forest is build to predict $ \hat{\tau} $.

K{\"{u}}nzel et al.\ propose X-Learner, which aims to improve upon the two-model approach~\citep{kunzel-pnas19}.
The goal of the X-Learner is to use information from the control group to learn better estimates for the treatment group and vice versa.
X-Learner consists of three steps.
The first step is to learn two supervised estimators $ \hat{\mu}_0 $ and $ \hat{\mu}_1 $.
This is the same as the two-model approach.
Step two involves imputing the treatment effects from predictions from the opposite group by taking the difference between factual and predicted counterfactual outcomes:
\begin{gather}
    \tilde{D}_i^1 = Y_i^{1} - \hat{\mu}\big(\mathbf{X}_i^{1} \big), \\
    \tilde{D}_i^0 = \hat{\mu}\big( \mathbf{X}_i^{0} \big) - Y_i^{0},
\end{gather}
where $ \tilde{D}_i^1 $ and $ \tilde{D}_i^0 $ is the imputed treatment effect for unit $ i $ if they were treated or untreated, respectively.
Then, two estimators $ \hat{\tau}_0(\mathbf{x}) $ and $ \hat{\tau}(\mathbf{x}) $, are trained on the imputed datasets, $ \tilde{D}_i^1 $ and $ \tilde{D}_i^0 $, respectively.
Finally, the CATE is estimated by a weighted average of the two estimators:
\begin{equation}
    \hat{\tau}(\mathbf{x}) = g(\mathbf{x}) \hat{\tau}_0(\mathbf{x}) + (1 - g(\mathbf{x})) \hat{\tau}_1(\mathbf{x}),
\end{equation}
where $ g \in [0,1]$ is a weight function.
K{\"{u}}nzel et al.\ suggest using an estimate of the propensity score for $ g $, so that $ g = \hat{e}(\mathbf{x}) $~\citep{kunzel-pnas19}.

Kennedy proposes a two-stage doubly robust method for HTE estimation called Double Robust Learner (DR-Learner)~\citep{kennedy-arxiv20}.
The training data is randomly separated into three sets: $ D_{1a}^k, D_{1b}^k, D_{2}^k $ each containing $ k $ observations.
In the first step, two estimators are trained on $ D_{1a}^k $ and $ D_{1b}^k $, called \emph{nuisance} training.
An estimator of the propensity scores, $ \hat{e} $, is constructed on $ D_{1a} $ and estimates of $ (\hat{\mu}_0 $ and $ \hat{\mu}_1) $ are constructed by performing prediction $ D_{1b} $.
Then, a pseudo-outcome is constructed as:
\begin{equation}
    \hat{\varphi}(\mathbf{X}, Y, T) = \frac{T - \hat{e}(\mathbf{X})}{\hat{e}(\mathbf{X}) (1 - \hat{e}(\mathbf{X})) } \Big( Y - \hat{\mu}_T(\mathbf{X}) \Big) + \hat{\mu}_1(\mathbf{X}) + \hat{\mu}_0(\mathbf{X}).
\end{equation}
Finally, regressing the pseudo-outcome on $ \mathbf{X} $ in the test sample $ D_{2}^k $ obtains estimations of the CATE\@:
\begin{equation}
    \hat{\tau}_{dr}(\mathbf{x}) = \mathbb{E}[\hat{\varphi}(\mathbf{X}, Y, T) \mid \mathbf{X}=\mathbf{x}].
\end{equation}
It is recommended to cross-fit with the three sets and use the average of the resulting CATE estimators.

\subsubsection{Tree and forest-based methods}

Tree and forest-based methods are tailored for direct HTE estimation without predicting potential outcomes.
Su et al.\ propose a tree-based method inspired by the CART algorithm for classification~\citep{brieman-cart84, su-jmlr09}.
They find binary covariate splits based on a t-test of differences between two potential children of a node.
Let \( \{\mu^L_1, \bar{y}^L_1, s_2^L, n_1\} \) be the population mean, sample mean, sample variance, and the sample size for the treated group in the left child node \(L\), and \( \{\mu^L_0, \bar{y}^L_0, s_2^L, n_2\} \) similarly for the control group.
Su et al.\ propose a measure for assessing the interaction between the split, \( s \) and the treatment by:
\begin{equation}
    g(s) = \frac{(\bar{y}^L_1 - \bar{y}^L_0) - (\bar{y}^R_1 - \bar{y}^R_0)}{\hat{\sigma} \sqrt{1/n_1 + 1/n_2 + 1/n_3 + 1/n_4}},
\end{equation}
where \( \hat{\sigma} = \sum_{i=1}^4 w_i s_i^2 \) is a pooled estimator of the constant variance, and \( w_i = (n_i - 1) / \sum_{j=1}^4 (n_j - 1) \).
For a split \( s \), \( F_t = g^2(s) \) converges to a \( \chi^2(1) \) distribution, which is used as the splitting criterion~\citep{su-jmlr09}.

Athey and Imbens introduce two splitting criteria for building trees based on mean squared error for predictions: \emph{adaptive} and \emph{honest}~\citep{athey-pnas16}.
The main goal of these trees is to maximize the observed treatment effects in subgroups to discover heterogeneity.
In a node \( \ell \), define the conditional mean outcomes under treated and control as:
\begin{gather}
    \hat{\mu}(t=1, \ell) = E[Y_i(t=1) \mid X_i \in \ell],
    \\
    \hat{\mu}(t=0, \ell) = E[Y_i(t=0) \mid X_i \in \ell].
\end{gather}
The average causal effect in a node \( \ell \) is defined as the difference in conditional means:
\begin{gather}
    \hat{\tau}(\ell) = \hat{\mu}(1, \ell) - \hat{\mu}(0, \ell).
\end{gather}
Then, the adaptive \emph{partition measure} for node \( \ell \), which is used as the score for a particular partition or node, can be defined as:
\begin{gather}
    F_a(\ell) = N_\ell \cdot \hat{\tau}^2(\ell).
\end{gather}
To build a causal tree, this partition measure is maximized at each potential split:
\begin{gather}
    \max_\ell \; \sum_{\ell} F_a.
\end{gather}

In the \emph{honest} approach, the training data is separated into two parts: \emph{training} and \emph{estimation}.
The training set is used to find splits, but causal effects are estimated using data in the estimation set.
In addition, a penalty on high variance in outcomes in a node is introduced to the splitting criteria.
Define \( V^{(1)}(\ell), V^{(0)}(\ell) \) be the variance in outcomes in node \( \ell \) and \( p \) to be the treatment probability in the training data.
The honest splitting criteria is defined as:
\begin{gather}
    F_h(\ell) = F_a - \Bigg( \frac{1}{N^{tr}} + \frac{1}{N^{est}} \Bigg)  \Bigg( \frac{V^{(1)}(\ell)}{p} + \frac{V^{(0)}(\ell)}{1-p} \Bigg)
\end{gather}
This approach balances finding heterogeneity and lowering variance in leaf nodes.

Tran and Zheleva introduced the idea of using validation sets for building causal trees~\citep{tran-aaai19}.
They propose three variants: CT-L, CT-HL, and CT-HV\@.
In CT-L, a separate validation set is used to penalize splits with estimated effects that do not match with the training data.
If we define the estimated causal effect in node \( \ell \) on the training data to be \( \hat{\tau}(\ell^{tr}) \) and the estimated causal effect on the validation data to be \( \hat{\tau}(\ell^{val}) \), then a cost term \( C \) can be defined as:
\begin{equation}
    C(\ell) =  N_\ell^{val} \cdot | \hat{\tau}(\ell^{tr}) - \hat{\tau}(\ell^{val}) |,
\end{equation}
where \( N_\ell^{val} \) is the size of the validation set in node \( \ell \).
The cost term measures the error of the estimated effect of training and validation data.
Using the cost term, the partition measure for CT-L is defined as:
\begin{gather}
    F_l(\ell): \frac{(1 - \lambda) \cdot F_a - \lambda \cdot C(\ell) }{| N_\ell^{tr} - N_\ell^{val}| + 1},
\end{gather}
where \( \lambda \) is used to control the strength of the cost function.
The other two variants, CT-HL and CT-HV, incorporate the honest penalty into building the tree.
In CT-HL, the data is split into three sets: training, estimation, and validation, and the honest penalty and validation cost are added together.
In CT-HV, the data is only split into two sets: training and validation.
The difference is that the validation set is also used as the honest estimation set, so less data needs to be split from the training data.

In addition to the introduction of CT-L, Tran and Zheleva introduce the concept of \emph{triggers} for heterogeneous treatment effects~\citep{tran-aaai19}.
When the treatment is an ordinal or monotonously increasing continuous variable, the effect depends on the amount of treatment.
For example, a clinician may be interested in understanding the minimum number of days or amount of treatment (the trigger) that patients need to take a medication to be cured~\citep{tran-aaai19}.
The goal becomes finding the threshold that maximizes the expected effect for subpopulations.

Wager and Athey proposed a random forest algorithm called Causal Forest for estimating heterogeneous treatment effects~\citep{wager-jasa18}.
Causal Forest constructs a forest of honest Causal Trees that each provides a CATE estimation. Then, an average of the predicted CATE from all trees is used as the final prediction.
They showed that the estimation of Causal Forest is asymptotically Gaussian and unbiased for the true CATE if the base learner used is \emph{honest}~\citep{wager-jasa18}.
Athey et al.\ extended this idea to Generalized Random Forests~\citep{athey-annals19}.
Generalized random forests extend random forests into a method for estimating any quantity through local moment conditions and can be used for heterogeneous treatment effect estimation.

\subsubsection{Deep learning approaches}

Several other deep learning approaches have also been developed based on different architectures.
Louizos et al.\ consider a different perspective and motivation for studying hidden representations for HTE estimation~\citep{louizos-neurips17}.
Instead of considering the given features $ \mathbf{X} $ as confounders, they consider it as noisy views on the hidden confounder $ \mathbf{Z} $.
Their goal is to learn the representation of $ \mathbf{Z} $ for causal effect adjustment using the observed features $ \mathbf{X} $.
To do this, they use variational autoencoders to learn $ \mathbf{Z} $ for the treatment and control groups, and use the hidden representation for prediction of outcomes.

Yoon et al.\ propose Generative Adversarial Networks for Individualized Treatment Effects (GANITE)~\citep{yoon-iclr18}.
GANITE has two components, a counterfactual block which generates counterfactual predictions, and an ITE block that generates CATE predictions given the counterfactual prediction.
This work leverages Generative Adversarial Networks (GANs) in each block to create a better counterfactual dataset.
After training the counterfactual block, generated data is fed to the ITE block which uses a GAN to learn to generate HTEs.
\subsection{Data-driven HTE estimation using structure learning}\label{sec:structure_learning}

As mentioned in Section~\ref{sec:hte_estimation_with_scms}, when the structural causal model is known, an adjustment set can be used to select features for HTE estimation.
However, in many real-world scenarios, the underlying causal model is unknown, and therefore, there is \emph{structure uncertainty} of the causal mechanisms.
Therefore, ignoring the underlying causal models can lead to incorrect causal effect estimations.

To use the causal model under structure uncertainty, one approach is to use causal structure learning~\citep{heinze-stat18,vowels-arxiv21}.
Many traditional causal structure learning methods can be broadly classified into constraint-based~\citep{spirtes-book00,colombo-jstor12}, score-based~\citep{chickering-jmlr02,hauser-jmlr12,nandy-annals18}, a hybrid of the two~\citep{tsamardinos-ml06}, or based on asymmetries~\citep{shimizu-jmlr06}.
These methods find a partial graph or skeleton from the data, with some unoriented edges.
With the output from structure learning methods, an adjustment set can be applied to select variables for estimation~\citep{pearl-book09,henckel-arxiv19}, and these variables are used to train an HTE estimator.

Although causal structure learning can address some problems of data-driven HTE estimation, there are some drawbacks.
First, learning the full structure can be computationally expensive for high-dimensional datasets~\citep{chickering-jmlr04,heinze-stat18,vowels-arxiv21}.
Second, causal structure learning generally cannot learn the full causal model, so learned models will have undirected edges with an ambiguous cause-effect relationship.
In addition, adjustment criteria have been developed for total causal effect estimation without heterogeneity in mind.
Therefore, adjustment sets can potentially miss heterogeneity-inducing variables, leading to incorrect HTE estimation.
In our experiments, we use randomly generated SCMs and explore how data-driven HTE estimators perform under different types of SCMs.
\section{Data driven HTE estimation}\label{sec:survey}
\section{Counterfactual prediction and effect estimation}\label{sec:survey}

This section provides a survey of existing heterogeneous treatment effect estimation methods.
We categorize the methods into two main categories: methods that perform counterfactual prediction and methods that estimate treatment effect directly.

The first category consists of methods that perform counterfactual prediction.
Counterfactual prediction methods predict counterfactual outcomes when an individual is treated and not treated~\citep{johansson-icml16,shalit-icml17,kunzel-pnas19}.
Then an estimate of an individual treatment effect is the difference in predicted outcomes.
The second category consists of methods that perform effect estimation directly.
These methods use counterfactual predictions as part of a process to learn an effect estimator~\citep{athey-pnas16,tran-aaai19,athey-annals19} or use different objective functions to estimate effects~\citep{kunzel-pnas19,louizos-neurips17}.
Table~\ref{tab:method_overview} shows subcategories of each category, together with representative papers.

\begin{table}[h]
    \centering
    \caption{An overview of HTE estimator categories and subcategories.}\label{tab:method_overview}
    \begin{tabular}[t]{M{0.47\linewidth}  M{0.47\linewidth}}
        Counterfactual Prediction & Effect Estimation \\
        \midrule
        Single-Model Approaches \citep{hill-jcgs11,kunzel-pnas19,athey-stat15, atan-aaai18} & Transformed Outcome-based Approaches \citep{rosenbaum-bio83,athey-stat15} \\
        {Two-Model Approaches} \citep{foster-statmed11,cai-bio11,athey-stat15,kunzel-pnas19, johansson-icml16,shalit-icml17}  & {Meta-Learning Approaches} \citep{xie-soc12,powers-statmed18,kunzel-pnas19,kennedy-arxiv20} \\
        {Multi-Model Approaches} \citep{schwab-arxiv19,shi-neurips19} & {Tree-based Approaches} \citep{athey-pnas16,tran-aaai19,su-jmlr09} \\
        & {Forest-based Approaches} \citep{wager-jasa18,athey-annals19,powers-statmed18} \\
        & {Deep Learning Approaches} \citep{yoon-iclr18,louizos-neurips17} \\
    \end{tabular}
\end{table}

\input{03-01-counterfactual.tex}

\input{03-02-effect_estimation.tex}

\subsection{Evaluation of HTE estimators}\label{sec:data_driven_background}

Since the true ITE is unknown in real-world scenarios, several heuristic HTE evaluation metrics have been developed for evaluating the performance of HTE estimators.
These metrics include inverse propensity weighting (IPW) validation~\citep{schuler-model18}, $ \tau $-risk~\citep{schuler-model18,nie-rlearner17}, nearest-neighbor imputation (NN-PEHE)~\citep{schwab-perfect18}, plugin validation~\citep{saito-cv19}, and counterfactual cross-validation (CFCV)~\citep{saito-cv19}.

IPW validation uses inverse propensity score weighting to a plugin estimate for the CATE, $ \tilde{\tau} $~\citep{schuler-model18}:
\begin{equation}
    \tilde{\tau}_{IPW} (\mathbf{X}_i, Y_i, T_i) = \frac{T_i}{\hat{e}(\mathbf{X}_i)}Y_i - \frac{1 - T_i}{1 - \hat{e}(\mathbf{X}_i)} Y_i.
\end{equation}
Here, $ \hat{e}(\mathbf{X}_i) $ is an estimate of the propensity score, which is the probability of receiving the treatment of interest (i.e., $\mathbb{E}[T_i|\mathbf{X}_i]$)
This plugin is used to approximate the mean-squared error for the estimator $ \hat{\tau} $:
\begin{equation}
    \hat{R}_{IPW}(\hat{\tau}) = \frac{1}{n} \sum_{i=1}^n \Big(\tilde{\tau}_{IPW} (\mathbf{X}_i, Y_i, T_i) - \hat{\tau}(\mathbf{X}_i) \Big)^2.
\end{equation}

The second metric $ \tau  $-risk~\citep{nie-rlearner17} takes as input an estimator $ \hat{\tau}  $ and approximates the objective defined in~\eqref{eq:cate_objective}.
To do this, two functions are defined: $ \hat{m}(\mathbf{X}_i) $ which is an outcome estimator for unit $i$ (i.e., $ \mathbb{E}[Y_i | \mathbf{X}_i] $) and $ \hat{e}(\mathbf{X}_i) $.
Both $ \hat{m} $ and $ \hat{e} $ are learned through supervised techniques and optimized for predictive accuracy~\citep{nie-rlearner17}.
The $ \tau  $-risk is defined as:
\begin{equation}
    \tau\text{-risk} = \frac{1}{N} \sum_{i=1}^N \Big( \big(Y_i-\hat{m}(\mathbf{X}_i)\big) - \big(T_i - \hat{e}(\mathbf{X}_i)\big)\hat{\tau}(\mathbf{X}_i) \Big)^2.
\end{equation}

The third metric, neighbor precision in estimation of heterogeneous effects (NN-PEHE)~\citep{schwab-perfect18} uses matching to approximate a treatment effect by substituting counterfactual outcomes from matched neighbors.
A unit $i$ is matched with one\footnote{Matching can be done with multiple units, and an aggregate can be computed.} control and one treated unit, denoted by NN$_0 (i)$ and NN$_1 (i)$.
Using the matched units, a substitute of the ground-truth treatment effect is computed as $ y(\text{NN}_1(i)) - y(\text{NN}_0(i)) $.
The NN-PEHE error, $ \hat{\epsilon}_{\text{NN-PEHE}} $ is computed as:
\begin{equation}
    \hat{\epsilon}_{\text{NN-PEHE}} = \frac{1}{N} \sum_{i=1}^n \bigg( \big(y(\text{NN}_1(i)) - y(\text{NN}_0(i)) \big) - \hat{\tau}(\mathbf{X}_i) \bigg)^2.
\end{equation}

Plugin validation~\citep{saito-cv19}, the fourth metric, uses machine learning algorithms to estimate a plugin $ \tau $, similar to IPW validation.
To get a plugin estimate, two predictions are needed, $ f_1(\mathbf{x}) $ and $ f_0(\mathbf{{x}}) $, which are predictions of the outcome when treated and not treated, respectively, using only the features.
Then, the plugin validation metric is defined as:
\begin{equation}
    \hat{R}_{plugin}(\hat{\tau}) = \frac{1}{N} \sum_{i=1}^N \bigg( (f_1(\mathbf{X}_i) - f_0(\mathbf{X}_i)) - \hat{\tau}(\mathbf{X}_i) \bigg)^2.
\end{equation}

Counterfactual cross-validation (CFCV) uses the plugin estimate of the CATE in a doubly robust formulation~\citep{saito-cv19}:
\begin{equation}
    \tilde{\tau}_{DR}(\mathbf{X}_i, Y_i, T_i) = \frac{T - \hat{e}(\mathbf{X}_i)}{\hat{e}(\mathbf{X}_i) (1 - \hat{e}(\mathbf{X}_i))} (Y_i - f_T(\mathbf{X}_i)) + f_1(\mathbf{X}_i) - f_0(\mathbf{X}_i),
\end{equation}
where $ \hat{e} $ is an estimate of the propensity score and $ f_t(\mathbf{x}) $ is an arbitrary regression function for a specific treatment value $ t $ that predicts the outcome given the set of features.
Then, the plugin estimate $ \tilde{\tau}_{DR} $ can be used in the mean-squared error estimate:
\begin{equation}
    \hat{R}_{CFCV} = \frac{1}{N} \sum_{i=1}^N \bigg( \tilde{\tau}_{DR}(\mathbf{X}_i, Y_i, T_i) - \hat{\tau}(\mathbf{X}_i) \bigg)^2
\end{equation}

Each metric is an estimate of the error of an HTE estimator and performance based on these heuristic metrics has been shown to correlate with true performance.
Since these metrics are heuristic, we explore their reliability in synthetic and semi-synthetic experiments in Section~\ref{sec:experiments}.

\section{Structural causal models for HTE estimation}\label{sec:structural_causal_models}

The data-driven estimation methods for counterfactual prediction and effect estimation we discussed so far do not explicitly capture how relationships between variables can affect the estimation of HTEs.
Instead, they typically rely on the strong ignorability assumption and use all features available in the data for estimation.
This can be problematic since including some variables could lead to incorrect estimation, which we discuss in more detail in Section \ref{sec:hte_estimation_with_scms}. To deal with this issue, it is helpful to know the cause-effect relationships between variables.
\emph{Structural causal models} (SCMs) are graphical models that capture these relationships and can help practitioners decide which variables should and should not be used for building HTE estimators.
One main drawback to SCMs is that they are rarely known in practice because the underlying causal mechanism that generated given data is unknown.
When the SCM is unknown, a causal structure learning algorithm can be used to learn a skeleton of the true model from the data itself.
In this section, we provide a brief introduction to SCMs and how HTEs are represented in them.

\subsection{Structural causal models}\label{sec:scm_background}

A structural causal model (SCM), $M$, consists of two sets of variables, $\mathbf{U}$ and $\mathbf{V}$, and a set of structural equations, $\mathbf{F}$, describing how values are assigned to each \emph{endogenous} variable $V^{(i)} \in \mathbf{V}$ based on the values of $\mathbf{v}$ and $\mathbf{u}$~\citep{pearl-book09}: \( v^{(i)} = f^{(i)}(\mathbf{v}, \mathbf{u}) \).
A causal graph \( G \) is a directed acyclic graph that captures the causal relationships among the variables.
Each variable, $V^{(i)}$, is represented by a node in the causal graph, an edge from $V^{(i)}$ to $V^{(j)}$ signifies that $V^{(i)}$ is a cause of $V^{(j)}$.
The edge relationship is a parent-child relationship, where the parent is a cause of the child.
The variables in $\mathbf{U}$ are considered \emph{exogenous}, and the causal graph may omit representing $\mathbf{U}$ explicitly.
Figure~\ref{fig:example_scm} shows an example of a causal graph with three exogenous variables, $ U_X, U_T, U_Y $ and three endogenous variables, $ X, T, Y $.

\begin{figure}
    \centering
    \newcommand{0.25}{0.5}

    \begin{subfigure}[b]{\columnwidth}
        \centering
        \begin{tikzpicture}
            \node (t) at (-0.75, 0) [label=below:$T$,point];
            \node (x) at (0.0, 0.75) [label=left:\(X\),point];
            \node(y) at (0.75, 0.0) [label=below:$Y$,point];
            \node (ux) at (0.0, 1.5) [label=above:\(U_X\),point];
            \node (ut) at (-1.5, 0.75) [label=above:\(U_T\),point];
            \node (uy) at (1.5, 0.75) [label=above:\(U_Y\),point];

            \path (x) edge (t);
            \path (x) edge (y);
            \path (t) edge (y);
            \path (ux) edge (x);
            \path (uy) edge (y);
            \path (ut) edge (t);
        \end{tikzpicture}
    \end{subfigure}

    \caption{An example causal model graph.
    }\label{fig:example_scm}
\end{figure}

An example set of associated structural equations is shown in eq~\eqref{eq:see1}-\eqref{eq:see3}.
\begin{gather}
    X = U_X \label{eq:see1} \\
    T = Bern(X) + U_T \label{eq:see2}  \\
    Y = \alpha \cdot T + \beta \cdot X + U_Y \label{eq:see3}
\end{gather}
In the graph, each endogenous variable has a directed edge pointing to it from an exogenous variable (e.g., $ U_X \rightarrow X $ or $ U_X $ is a parent of $ X $).
This is represented in the structural equations as being part of the structural equation (e.g., $ U_X $ in eq~\eqref{eq:see1}).
It is assumed that every endogenous variable has a parent in $ \mathbf{U} $.
We denote the parents and children of a variable, $ X $, in the graph as $ pa(X) $ and $ ch(X) $, respectively.
Ancestors of $ X $, $ an(X) $, are variables that can reach $ X $ through a directed path, and descendants of $ X $, $ de(X) $, are variables that can be reached from $ X $ through directed paths out of $ X $.

\subsubsection{Identifying effects using SCMs}\label{sec:identifying_scms}

When the SCM underlying the dataset is known, its graph $ G $ supports the use of graphical criteria for variable selection for unbiased causal effect estimation, known as adjustment set selection~\citep{shpitser-uai06,shpitser-uai10,van-uai14,perkovic-uai15,henckel-arxiv19}.
These graphical criteria rely on the concept of $d$-separation.
\begin{definition}
    ($d$-separation) A set of $ \mathbf{Z} $ nodes is said to ``block'' a path $ p $ if either
    \begin{enumerate}
        \item[$(i)$] $ p $ contains at least one arrow-emitting node that is in $ \mathbf{Z} $ or
        \item[$(ii)$] $p$ contains at least one collision node that is outside $ \mathbf{Z} $ and has no descendant in $ \mathbf{Z} $.
    \end{enumerate}
\end{definition}
If $ \mathbf{Z} $ blocks all paths from a set of nodes $ \mathbf{X} $ to a set of nodes $ \mathbf{Y} $, it is said to ``d-separate'' $ \mathbf{X} $ and $ \mathbf{Y} $, and the conditional independence, $ \mathbf{X} \perp \mathbf{Y} \mid \mathbf{Z} $, holds~\citep{bareinboim-pnas16}.
One graphical criteria for adjustment set selection is the backdoor criterion~\citep{pearl-book09}.
The backdoor criterion aims to find a set that \emph{blocks} the paths from the treatment variable to the outcome variable in the graph, starting with an arrow into the treatment variable.
Essentially, the backdoor criterion blocks all the paths going backward out of treatment.
For example, in Figure~\ref{fig:example_scm}, there is one backdoor path, namely $ T \leftarrow X \rightarrow Y $, and so the adjustment set $ \{X\} $ sufficiently satisfies the backdoor criterion.
Several adjustment criteria have been developed, such as the adjustment criterion~\citep{shpitser-uai06,van-uai14}, generalized backdoor criterion~\citep{maathuis-annals15}, generalized adjustment criterion~\citep{perkovic-uai15}, and optimal adjustment set (O-set)~\citep{henckel-arxiv19}.
These adjustment sets defined by the various adjustment criteria can be found through \emph{identification} algorithms such as the ID algorithm~\citep{shpitser-jmlr08}, the algorithmic framework by Van der Zander et al.~\citep{van-uai14} or the pruning procedure by Henckel et al.~\citep{henckel-arxiv19}.

When an adjustment set is discovered, unbiased causal effects can be estimated by conditioning on the adjustment set.
The adjustment set can be used instead of all the features for data-driven estimation.
We discuss data-driven HTE estimation with SCMs briefly in Section~\ref{sec:hte_estimation_with_scms}.

\subsection{Heterogeneous treatment effects in SCMs}

Shpitser and Pearl provide a nonparametric estimand for HTEs for some feature assignment $C=c$, called $c$-specific effects~\citep{shpitser-uai06,pearl-smr17}.
One issue with estimating $c$-specific effects is that it requires knowledge of the SCM and which variables, C, exhibit heterogeneous effects.
When the SCM is known, the O-set~\citep{henckel-arxiv19} can be used on the learned structure since it implicitly includes all potential heterogeneity-inducing variables, while some identification algorithms may miss them. When the SCM is unknown, the identifiability of \( c \)-specific effects cannot be established.
Moreover, when the set of heterogeneity-inducing variables $\mathbf{C}$ is unknown, then all combinations of features and feature values would need to be checked for heterogeneity in effect, which is exponential in the number of valid features.
Circumventing this exponential search is one of the main reasons data-driven methods are popular in practice.

Shpitser and Pearl provide a nonparametric estimand for HTEs with respect to some feature assignment $C=c$, called $c$-specific effects~\citep{shpitser-uai06,pearl-smr17}.
There are three cases when $c$-specific effects are \emph{identifiable}, i.e., there exists an unbiased estimation for $\tau(\mathbf{c}) = E[Y(1) - Y(0) | \mathbf{C}=\mathbf{c}]$.
The first case is when $\mathbf{C}$ is observed and admissible which means that it satisfies the backdoor criterion in the causal graph~\citep{pearl-book09}.
A set of variables $\mathbf{Z}$ satisfies the backdoor criterion relative to $(T,Y)$ if no node in $\mathbf{Z}$ is a descendant of $T$ and $\mathbf{Z}$ blocks every path between $T$ and $Y$ that contains an arrow into $T$.
For example, $X$ is the only variable that is admissible in all SCMs in Figure~\ref{fig:fail_examples}, since it blocks the only backdoor path.

The second case is when $\mathbf{C}$ is part of an observed admissible set.
In this case, $\mathbf{C}$ itself is not admissible, but we can observe a set of features $\mathbf{S}$ such that $\mathbf{S} \cup \mathbf{C}$ is admissible.
In Figure~\ref{fig:marloes_scm}, $\{F\}$ could induce heterogeneity but is not admissible on its own. Adding $F$ to $\{X\}$ still satisfies the backdoor criterion making $\{X, F\}$ also admissible.
The third case for identifiable $c$-specific effects is when $\mathbf{C}$ is not part of an admissible set.
The first two cases which rely on \emph{covariate adjustment}~\citep{perkovic-jmlr18} and identifying heterogeneity-inducing variables are the basis of our work.
Estimating heterogeneous effects in the third case is beyond the scope of this paper and is reserved for future work.

One issue with estimating $c$-specific effects is that it requires knowledge of the SCM and which variables exhibit heterogeneous effects.
If the SCM is unknown, the identifiability of \( c \)-specific effects cannot be established.
If $\mathbf{C}$ is unknown, then all combinations of features and feature values would need to be checked for heterogeneity in effect, which is exponential in the number of valid features~\citep{pearl-smr17}.
Circumventing this exponential search is one of the main reasons why data-driven methods are popular in practice.

\subsubsection{Detecting latent heterogeneity}

In some cases, the observed covariates that induce heterogeneity, $C$, are not known, but there is unknown heterogeneity in effects present in the data~\citep{pearl-smr17}.
Pearl discusses how heterogeneity can be detected from statistical averages alone~\citep{pearl-smr17} when the effect of treatment on the treated ($ETT$) and the effect of treatment on the untreated ($ETU$) are identifiable.
The three cases where $ETT$ can be identified are when:
\begin{itemize}
    \item The treatment is binary, and $E[Y_1]$ and $E[Y_0]$ are identifiable by some method
    \item The treatment is arbitrary, and $E[Y_x]$ is identifiable (for all $x$) by adjustment for an admissible set of covariates.
    \item $ATE$ is identified through mediating instruments.
\end{itemize}
The data-driven HTE estimation methods discussed in this article do not focus on the detection of latent heterogeneity, but using these cases when latent heterogeneity can be detected with data-driven estimators is an interesting future direction.

\subsubsection{Representation of HTEs in structural equations}\label{sec:htes_scms_background}

In SCMs, heterogeneous effects are represented by an \emph{interaction} or \emph{effect modification} through ancestors and parents of a variable.
We will use linear structural equations for illustration.
For heterogeneous effects, this interaction is specifically with the treatment variable.
Consider the causal model in Figure~\ref{fig:example_scm}, where the treatment variable of interest is $ T $ and the outcome of interest is $ Y $.
There is a potential interaction between $ X $ and $ T $.
Suppose $ Y $ is generated by the following structural equation:
\begin{equation}
    Y = \alpha \cdot T + \beta \cdot X + \gamma \cdot T \cdot X + U_Y.
\end{equation}
Here, there is a fixed effect \( \alpha \) of the treatment \( T \) and a heterogeneous effect \( \gamma \) due to the interaction between \( T \) and \( X \).
Since the effect of treatment varies with another variable, $ X $, there is a heterogeneous effect.
Heterogeneity through interaction can occur with any parent of the outcome or any parent of \emph{mediating} variables~\citep{vanderweele-epid07}.

\begin{figure}
    \centering
    \begin{tikzpicture}
        \node (t) at (-1.5, 0) [label={below left}:$T$,point];
        \node (x) at (-1.5, 1.5) [label={above}:$X$,point];

        \node (c) at (0.0, 1.5) [label={above}:$C$,point];
        \node (d) at (0.0, 0) [label={below}:$D$,point];

        \node (f) at (1.5, 1.5) [label={above}:$F$,point];
        \node (y) at (1.5, 0.0) [label={below}:$Y$,point];

        \path (x) edge (t);

        \path (x) edge (y);
        \path(t)edge(d);
        \path(c)edge(d);
        \path (d) edge (y);

        \path (f) edge (y);
    \end{tikzpicture}
    \caption{Causal model with multiple possible interactions that cause heterogeneous treatment effects.}\label{fig:multi_interaction}
\end{figure}

Consider another example causal model, shown in Figure~\ref{fig:multi_interaction}.
Here, there is a \emph{mediating} variable, $ D $, a variable that occurs on a path from the treatment, $ T $, to the outcome, $ Y $.
Interaction can occur with $ X, C$, or $ F $.
For example, an interaction with $ X $ and $ F $ can be shown as:
\begin{equation}
    Y = \alpha \cdot D \cdot X + \beta \cdot D \cdot F + U_Y.
\end{equation}
Since the treatment is mediated through $ D $ and $ X $ and $ F $ interact with $ D $, there will be a heterogeneous effect.
Since $ D $ is a mediator, if $ C $ interacts with $ T $ in the structural equation, then there will also be a heterogeneous effect, such as:
\begin{gather}
    D = \alpha \cdot T + \beta \cdot T \cdot C + U_D, \\
    Y = \gamma \cdot D + U_Y.
\end{gather}

\subsection{Data-driven HTE estimation using Structural Causal Models}\label{sec:hte_estimation_with_scms}

\begin{figure}
    \centering
    \newcommand{0.25}{0.3}
    \resizebox*{0.25\columnwidth}{!}{
        \begin{subfigure}[b]{0.325\columnwidth}
            \centering
            \begin{tikzpicture}
                \node (t) at (-0.75, 0) [label=left:$T$,point];
                \node (x) at (0.0, 0.75) [label=above:\(X\),point];
                \node(y) at (0.75, 0.0) [label=right:$Y$,point];
                \node(c) at (0.0, -0.75) [label=right:$L$,point];

                \path (x) edge (t);
                \path (x) edge (y);
                \path (t) edge (y);
                \path (t) edge (c);
                \path (y) edge (c);
            \end{tikzpicture}
            \caption{}\label{fig:collider_scm}
        \end{subfigure}
    }
    \resizebox*{0.25\columnwidth}{!}{
        \begin{subfigure}[b]{0.325\columnwidth}
            \centering
            \begin{tikzpicture}
                \node (t) at (-0.75, 0) [label=left:$T$,point];
                \node (m) at (0.0, 0) [label=below:$M$,point];
                \node (x) at (0.0, 0.75) [label=above:\(X\),point];
                \node(y) at (0.75, 0.0) [label=right:$Y$,point];

                \path (x) edge (t);
                \path (x) edge (y);
                \path (t) edge (m);
                \path (m) edge (y);
            \end{tikzpicture}
            \caption{}\label{fig:mediator_scm}

        \end{subfigure}
    }
    \resizebox*{0.25\columnwidth}{!}{
        \begin{subfigure}[b]{0.325\columnwidth}
            \centering
            \begin{tikzpicture}
                \node (a) at (-1.5, 0) [label=left:$A$,point];
                \node (t) at (-0.75, 0) [label={below left}:$T$,point];
                \node (b) at (-0.75, -0.75) [label={below}:$B$,point];
                \node (x) at (-0.75, 0.75) [label={above}:$X$,point];

                \node (c) at (0.0, 0.75) [label={above}:$C$,point];
                \node (d) at (0.0, 0) [label={below left}:$D$,point];
                \node (e) at (0.0, -0.75) [label={below}:$E$,point];

                \node (f) at (0.75, 0.75) [label={above}:$F$,point];
                \node (y) at (0.75, 0.0) [label={right}:$Y$,point];
                \node (g) at (0.75, -0.75) [label={below}:$G$,point];

                \path (a) edge (t);
                \path (x) edge (t);
                \path (t) edge (b);

                \path (x) edge (y);
                \path(t)edge(d);
                \path(c)edge(d);
                \path(d)edge(e);
                \path (d) edge (y);

                \path (f) edge (y);
                \path (y) edge (g);
            \end{tikzpicture}
            \caption{}\label{fig:marloes_scm}
        \end{subfigure}
    }
    \caption{Causal models where data-driven HTE estimation methods may perform poorly if the structure is unknown.
    Some variables (e.g., $L$, $M$, and $E$) are not valid adjustment variables for estimating the effect of $T$ on $Y$.
    }\label{fig:fail_examples}
\end{figure}

As previously mentioned, data-driven methods typically use all features available in the data for estimation.
However, not all variables in the data should be included in HTE estimation, and including incorrect variables leads to incorrect effect estimation~\citep{pearl-book09}.
Incorrect variables include descendants of treatment or instrumental variables~\citep{pearl-book09,myers-aje11} (e.g., $L, M, A, B, D, E, G$ in Figure~\ref{fig:fail_examples}).
If the underlying causal model is known, then an adjustment set can be discovered using an adjustment criterion.
Then, unbiased effects can be estimated by using the adjustment set as the features for estimation.

\subsubsection{Causal structure learning}\label{sec:structure_learning}

In many real-world scenarios, the underlying causal model is unknown, and therefore, there is \emph{structure uncertainty} of the causal mechanisms.
Therefore, ignoring the underlying causal models can lead to incorrect causal effect estimations.

To use the causal model under structure uncertainty, one approach is to use causal structure learning~\citep{heinze-stat18,vowels-arxiv21}.
Many traditional causal structure learning methods can be broadly classified into constraint-based~\citep{spirtes-book00,colombo-jstor12}, score-based~\citep{chickering-jmlr02,hauser-jmlr12,nandy-annals18}, a hybrid of the two~\citep{tsamardinos-ml06}, or based on asymmetries~\citep{shimizu-jmlr06}.
These methods find a partial graph or skeleton from the data, with some unoriented edges.
With the output from structure learning methods, an adjustment set can be applied to select variables for estimation~\citep{pearl-book09,henckel-arxiv19}, and these variables are used to train an HTE estimator.

One well-known constraint-based structure learning algorithm is the PC algorithm~\citep{spirtes-book00}.
PC conducts conditional independence tests to learn about the structure of the underlying model. It learns a completed partially directed acyclic graph (CPDAG)~\citep{chickering-jmlr02}, which is an equivalence class of directed acyclic graphs (DAG) that have uncertain edges.
PC has three main steps: (1) determining the skeleton, (2) determining v-structures, and (3), determining further edge orientations.
PC learns the skeleton by starting with a complete undirected graph.
Using a series of conditional independence tests for adjacent nodes, PC removes an edge if a conditional independence is found.
In the second step, all edges in the graph are replaced with $\circ$--$\circ$ edges.
PC considers all unshielded triples, triples $i \circ $--$\circ j$--$\circ k$ where $i$ and $k$ are not adjacent and determines whether a v-structure should exist or not.
In the final step, additional orientation rules are applied to orient as many unoriented edges as possible.
The fast causal inference (FCI) algorithm is a modification of the PC algorithm that allows hidden variables~\citep{spirtes-book00}.

A classic score-based algorithm is greedy equivalence search (GES) which also learns a CPDAG.
Score-based algorithms typically use a likelihood score of a DAG against the data, and search for a DAG or CPDAG that yields the optimal score.
GES learns the CPDAG via a greedy forward search of edge additions, followed by a backward phase with edge deletions.
A hybrid approach based on FCI and GES has also been developed, called greedy fast causal inference (GFCI)~\citep{ogarrio-pgm16}.
Heinze-Deml et al.\ provide a survey on structure discovery methods, specifically constraint-based and score-based approaches~\citep{heinze-stat18}.

Besides constraint-based approaches and score-based approaches, continuous optimization-based methods for structure discovery have been developed due to increased compute that make learning from high-dimensional datasets feasible~\citep{ma-sr14,zheng-neurips18,kyono-neurips20,vowels-arxiv21}.
Continuous optimization-based methods recast the graph-search problem into a continuous optimization problem~\citep{vowels-arxiv21}.
Vowels et al.\ provide a survey on continuous optimization-based approaches~\citep{vowels-arxiv21}.
With the learned causal graph, an \emph{identification} algorithm can be used to find an adjustment set.
Then, this adjustment set can be used as features in an HTE estimator instead of the full feature set.

Although causal structure learning can address some problems of data-driven HTE estimation, there are some drawbacks.
First, learning the full structure can be computationally expensive for high-dimensional datasets~\citep{chickering-jmlr04,heinze-stat18,vowels-arxiv21}.
Second, causal structure learning generally cannot learn the full causal model, so learned models will have undirected edges with an ambiguous cause-effect relationship.
In addition, adjustment criteria have been developed for total causal effect estimation without heterogeneity in mind.
Therefore, adjustment sets can potentially miss heterogeneity-inducing variables, leading to incorrect HTE estimation.
In our experiments, we use randomly generated SCMs and explore how data-driven HTE estimators perform under different types of SCMs.
\subsection{Synthetic dataset evaluation}\label{sec:synthetic_eval}

We first investigate how HTE estimation methods perform under different parameter settings (and thus SCMs).
Then we investigate the effect that causal feature selection has on HTE estimation.

\subsubsection{Multidimensional scaling}\label{sec:synthetic_mds}

\input{figures/mds_correlation.tex}

For each method and SCM, we compute the RMSE for the predicted effect.
This yields a (number of methods) $ \times $ (number of SCMs) matrix with RMSE values.
The correlation between pairs of SCM results is the similarity between two datasets.
Figure~\ref{fig:mds_correlation} shows a Multidimensional scaling (MDS) plot based on the correlation between SCM datasets.
The goal of the MDS plot is to visualize the similarity in performance of methods based on different SCM settings.
If two points are close together, then all learners in those experiments performed similarly.
The color and shape of a point represent the method that performed the best (lowest RMSE) for that SCM\@.

From this Figure, Causal Tree has the best performance in most settings, followed by Causal Forest and Balancing Neural Network.
Surprisingly, Causal Tree and Causal Forest are separated, meaning that the settings in which they perform the best have different RMSE results across all learners.
Balancing Neural Network occupies two regions close to Causal Tree and Causal Forest, with a majority near Causal Forest and TARNet.
This makes sense since TARNet is an extension of Balancing Neural Network.

\input{figures/pairwise_table.tex}

\subsubsection{Pairwise comparison}

Next, we investigate if there are methods that beat others consistently.
We compute the percentage of wins in Table~\ref{tab:pairwise} for which there is a 5\% decrease in RMSE\@.
For example, SLearner beats XLearner by a 5\% margin in only 0.34\% of scenarios, while XLearner beats SLearner by a 5\% margin in 50.12\% of scenarios.
From these results, there does not appear to be an overly dominant method.
For example, Causal Forest generally beats DRLearner by the margin but does not generally beat by Causal Tree or Balancing Neural Network.

From this Table, there is no clear dominant method, but there are methods which are dominated by others.
For example, SLearner generally does not perform well, especially its other meta-learning counterparts.
Causal Tree has an overall high percentage of wins but is also beaten by other methods at a relatively high percentage.
This suggests methods need to be carefully chosen for different settings and assumptions.

\subsubsection{Performance in different SCMs}

\input{figures/winner_table.tex}

In order to investigate how HTE estimators perform in different types of SCMs, we show the average rank under different parameter settings in Table~\ref{tab:mse_rankings} and the average RMSE in Table~\ref{tab:rmse_table}.
The Oracle method (Known Structure) is at the bottom of each table.
Bolded results indicate the best rank in the column (ignoring the Oracle).
In terms of ranking, structure learning (FSL) achieves the lowest average rank over most datasets.
Additionally, since some ranks are close, there is no clear winner among closely ranked methods, which supports the pairwise comparison between HTE estimation methods in Table~\ref{tab:pairwise}.

\input{figures/mse_table.tex}

We also show the average RMSE for each method in the same parameter settings in Table~\ref{tab:rmse_table}.
Here, Causal Tree and Balancing Neural Network achieve lower average RMSE.
DRLearner has high error and variance in predictions, which suggests it requires more fine-tuning of its estimators.

Looking at specific parameters, we see that SCMs with mediators generally result in higher errors than without mediators.
This makes sense since including mediators will bias the effect estimation.
Interestingly, a chain length of 1 is more challenging than a chain length of 2.
Increasing the number of HTE inducing parents also results in higher errors since 0 HTE parents imply average effect only, and more HTE inducing parents increases the complexity of the heterogeneous effect.
Causal Tree seems to be the least affected by parameter changes compared to other HTE estimators.

\subsection{Accuracy of decisions}

In practice, practitioners may be interested in how well an HTE estimator can assign a treatment decision to individuals, rather than the magnitude in the change in outcomes.
To investigate this, we explore how well an HTE estimator can decide whether it is beneficial to treat an individual or not.
To do this, we binarize the ground truth treatment effect into positive and negative effects.
For each HTE estimator, we binarize the predictions into positive and negative predictions.
Then, an accuracy is computed over multiple datasets, signifying how many examples are labeled correctly.

Table~\ref{tab:accuracy_table} shows the accuracy of decisions for each HTE estimator, using the same columns as in Tables~\ref{tab:mse_rankings} and~\ref{tab:rmse_table}.
In general, we see that Causal Tree achieves the highest accuracy of decisions, followed by Balancing Neural Network and FSL.
Meta-learners achieve very low accuracy when there are more invalid variables, such as when there are mediators.

\input{figures/accuracy_table.tex}

\subsection{Comparing evaluation metrics}

Since ground truth causal effects are generally not available in real-world datasets, several heuristic HTE evaluation metrics have been developed.
We consider inverse propensity weighting (IPW)~\citep{schuler-model18}, $\tau$-risk~\citep{schuler-model18,nie-rlearner17}, plugin validation (plugin $\tau$)~\citep{saito-cv19}, and counterfactual cross-validation (CFCV)~\citep{saito-cv19}.
Each metric uses estimators for propensity scores or outcomes to compute an error.
Each metric uses the full set of features in our experiments to compute a propensity score or outcome.
We use Logistic Regression for each metric that requires an estimate of the propensity score.
For $\tau$-risk, we use Linear Regression for outcome prediction.
For plugin validation and CFCV, we use Counterfactual Regression (CFR)~\citep{shalit-icml17} for the regression function.
In order to evaluate the reliability of these heuristic metrics, we compute the Spearman rank correlation of each metric compared to the ground truth ranking.
For each dataset, we rank each learner based on the heuristic error on the test set.
Then we compute a Spearman correlation in one of two ways.
The first way computes a correlation on the entire ranking for each dataset, and an average can be computed over all datasets.
The second way computes a correlation of the rank for each learner over all the datasets.

Table~\ref{tab:metric_dataset} shows the correlation of each heuristic metric ranking to the ground truth ranking over the synthetic datasets and the semi-synthetic (IHDP) datasets.
This table shows that most metrics are positively correlated with the ranking from the ground truth on synthetic datasets, except for plugin validation.
However, there is a large variance, and the average correlation is generally low, which suggests that having a hidden structural causal model with invalid variables may affect the metrics.
For the IHDP dataset, we see that all metrics except plugin validation positively correlate with the true ranking, with $ \tau  $-risk having the highest correlation.
Since the IHDP dataset does not contain mediators or descendants of treatment, this may improve the correlation of $ \tau  $-risk.

\begin{table*}[t]
    \centering

    \resizebox*{\linewidth}{!}{
        \begin{tabular}{l *{8}{M{1.7cm}}}
            & \multicolumn{7}{c}{Synthetic} \\
            \cmidrule(lr){2-8}

            & Overall & Has Mediator & No Mediator & Contains Confounding & No Confounding & 0 HTE Parents & 1 or 2 HTE Parents &  IHDP \\
            \midrule
            IPW &  0.10 $ \pm $ 0.55 & 0.03 $\pm$ 0.59 & 0.01 $\pm$ 0.43 & 0.02 $\pm$ 0.44 & 0.08 $\pm$ 0.45 & 0.11 $\pm$ 0.42 & 0.03 $\pm$ 0.44 & 0.10 $ \pm $ 0.43 \\
            $ \tau $-risk & 0.20 $\pm$ 0.47 & 0.01 $\pm$ 0.45 & 0.10 $\pm$ 0.43 & 0.13 $\pm$ 0.44 & 0.07 $\pm$ 0.43 & 0.13 $\pm$ 0.47 & 0.09 $\pm$ 0.43 & 0.53 $ \pm $ 0.32\\
            Plugin $\tau$ & -0.18 $\pm$ 0.55 & -0.23 $\pm$ 0.56 & -0.27 $\pm$ 0.44 & -0.27 $\pm$ 0.43 & -0.25 $\pm$ 0.44 & -0.24 $\pm$ 0.45 & -0.26 $\pm$ 0.43 & -0.34 $ \pm $ 0.38\\
            CFCV & 0.25 $\pm$ 0.59 & 0.14 $\pm$ 0.61 & 0.06 $\pm$ 0.47 & 0.03 $\pm$ 0.47 & 0.03 $\pm$ 0.45 & 0.15 $\pm$ 0.56 & 0.01 $\pm$ 0.47 & 0.20 $ \pm $ 0.44 \\
        \end{tabular}
    }
    \caption{Comparing the rank correlation of metrics with true performance.}\label{tab:metric_dataset}
    \vspace{1em}
\end{table*}
\subsection{Semi-synthetic}

The IHDP dataset is a semi-synthetic dataset, where covariates are used to generate simulated outcome so that ground truth effects are available for evaluation.
We show two results in Table~\ref{tab:ihdp_rankings}.
The first column of Table~\ref{tab:ihdp_rankings} shows the average rankings on the IHDP dataset in terms of RMSE.
From this Table, we see that FSL, TLearner, and XLearner perform consistently well in this dataset in terms of RMSE compared to all other estimator, since their rank is much higher than the next best performing estimator, SLearner.
The second column shows the $ \tau  $-risk results.
We see that the $ \tau  $-risk rank is generally correlated with the RMSE rank of each learner.
A notable difference is that Causal Forest is ranked higher based on $ \tau  $-risk, but lower based on RMSE.

\input{figures/ihdp_table.tex}
\subsection{Real-world datasets}

In real-world datasets, ground truth effects are not available.
Instead, we opt to compare HTE estimators based on $ \tau  $-risk as a proxy for methods perform in different datasets.
The average rank for each method in terms of $\tau$-risk for the LoL dataset is shown in Table~\ref{tab:ihdp_rankings}.
In this dataset, Causal Forest outranks other learners by a significant margin, followed by XLearner.
The rank in terms of $\tau$-risk for the Cannabis dataset is shown in Table~\ref{tab:ihdp_rankings}.
In the Cannabis dataset, XLearner achieves the highest $ \tau $-risk rank consistently.
Finally, the $\tau$-risk for the SHARPS dataset is shown in the last column of Table~\ref{tab:ihdp_rankings}.
In the SHARPS dataset, SLearner achieves the best rank consistently.
An interesting observation from this dataset is that the average difference in outcomes between treated and non-treated individuals is a small number (close to zero).
As pointed out by K{\"{u}}nzel et al., SLearner has a tendency to shrink estimates toward zero~\citep{kunzel-pnas19}, which may play a role in the $\tau$-risk rank in this dataset.
\section{Experiments}\label{sec:experiments}

We evaluate HTE estimators on synthetic, semi-synthetic, and real-world experiments.
We test several prominent HTE estimators in our experiments.

\subsection{Experimental setup}\label{sec:setup}

In our experiments, we compare several HTE estimators.
We employ several methods from recent work:
Causal Trees~\citep{tran-aaai19}, Causal Forest~\citep{athey-annals19},
Meta-learners (S, T, XLearners)~\citep{kunzel-pnas19}, Doubly Robust learner (DRLearner)~\citep{kennedy-arxiv20}, Balancing Neural Network~\citep{johansson-icml16}, and TARNet~\citep{shalit-icml17}.
In addition, we also use a structure learning method to find the optimal adjustment set (O-set)~\citep{henckel-arxiv19}.
We opt to use a hybrid method, greedy fast causal inference (GFCI)~\citep{ogarrio-pgm16}, and abbreviate this method as full structure learning with O-set (FSL (O-set)).
For synthetic datasets, where the causal model is known, we use an Oracle method (Known Structure) that uses the causal model for variable selection using the \emph{optimal adjustment set (O-set)}~\citep{henckel-arxiv19} and a TLearner for estimation.
We use Logistic Regression and Linear Regression for propensity scores and estimating outcomes for meta-learners on synthetic datasets.
For semi-synthetic and real-world datasets, we use Random Forest for estimating outcomes.

\subsection{Datasets}
We utilize three types of datasets.
The first type of dataset is a synthetic dataset.
The second type of dataset is a semi-synthetic dataset, the Infant Health and Development Program (IHDP) dataset~\citep{hill-jcgs11}.
The final type of dataset are real-world datasets.
We use League of Legends, Cannabis, and SHARPS.

\subsubsection{Synthetic dataset}\label{sec:setup_synthetic}

We first evaluate HTE estimators on synthetic datasets generated by random structural causal models (SCMs).
As mentioned in Section~\ref{sec:hte_estimation_with_scms}, not all variables in the data should be included in HTE estimation.
A current gap in the literature is understanding how data-driven HTE estimators perform when variables should not be included due to the underlying structure, such as mediators or descendants of treatment.
Therefore, we seek to understand how estimators perform under differing structures in our synthetic experiments.
In this setting, we assume estimators do not have access to the underlying causal graph and instead use all features in the generated dataset, which mimics the real-world use of these estimators.

We generate SCMs by following a procedure similar to one for evaluating causal structure learning algorithms~\citep{heinze-stat18}.
SCMs and their corresponding datasets differ in the following characteristics: the number of variables \( d \); the probability of edge creation \( p_e \); the variance in the noise \( \sigma \); the strength of the noise term in non-source nodes \( \rho \).
In addition, we add parameters to control confounding \( \gamma \), mediation \( m \), and heterogeneity in treatment effects \( p_h, m_p \).
We initialize a set of parameters as follows:

\begin{enumerate}
    \item Number of variables: \( d \in [10, 20, 30] \)
    \item Edge probability parameter: \( p_e \in [0.1, 0.3, 0.5] \)
    \item Noise variance: \( \sigma \in [0.2, 0.4, 0.6] \)
    \item Magnitude in noise: \( \rho \in [0.1, 0.5, 0.9] \)
    \item Confounder: \( \gamma \in [\text{True, False}] \)
    \item Mediator chain length: \( m \in [0, 1, 2] \)
    \item Number of HTE inducing parents: \( p_h \in [0, 1, 2] \)
    \item HTE from mediating parents: \( m_p \in [\text{True, False}] \)
\end{enumerate}

Here, we describe how we generate synthetic datasets used in the experiments.
We first generate an adjacency matrix, \( \mathbf{A} \) over the \( d \) variables.
Assuming that the variables \( \{1, \dots, d \} \) are causally ordered.
For each pair of nodes \( i, j \), where \( i < j \) in the causal order, we sample from Bernoulli\(( p_e )\) to determine if there is an edge from \( i \) to \( j \).
With the adjacency matrix, we sample the corresponding coefficients, \( \mathbf{C} \), from Unif\(( -1, 1 )\).
Noise terms \( \epsilon_i \) are generated by sampling from a normal distribution with mean zero and variance \( \sigma \).
For a non-source node, $ i $ (nodes with parents), the noise term is multiplied by the magnitude in noise, $ \gamma $.

Next, we choose the treatment and outcome nodes based on two parameters, \( \gamma \) and \( m \). \( \gamma \) controls whether there is confounding in the generated dataset and \( m \) controls the length of the mediating chain (\( 0 \) means no mediator).
First, we select all pairs of nodes with a \emph{directed} \( m \) hop path.
Then, if \( \gamma \) is True, we filter all pairs that do not contain a \emph{backdoor path}.
Otherwise, we filter out pairs that contain a backdoor path.
Finally, we sample a pair of nodes from the filtered set to obtain the treatment and outcome nodes.
If no pairs satisfy the parameter setting (e.g., no confounding), the graph is resampled up to 100 times.
Otherwise, we discard the parameter setting.

A source node \( i \) (node with no parents) takes the value of the sampled noise terms, \( \epsilon_i \).
For a non-source node, \( i \), we consider several additive factors.
First is the base term from the coefficient matrix: \( C_{pa(i), i} X_{pa(i)} \), where \( C_{pa(i), i} \) are the coefficients from \( pa(i) \) to \( i \) and \( X_{pa(i)} \) are the values generated for \( pa(i) \).
Next are the interaction terms, which are determined by the mediating chain length \(m\), number of HTE inducing parents \( p_h \), and whether there is HTE from parents of mediators \( m_p \).
For mediators, we have: \( \mathbbm{1}_{m}(i) \mathbbm{1}_{m_p} X_{pa_m(i)} X_{pa_n(i)} \), where \( \mathbbm{1}_{m} \) is an indicator for if \( i \) is a mediator, \( \mathbbm{1}_{m_p} \) is an indicator for if there is HTE induced from mediator parents, and \( pa_m \) and \( pa_n \) represent mediating and non-mediating parents.
For outcomes, the interaction is defined similarly: \( \mathbbm{1}_{y}(i) X_{pa_m(i)} X_{pa_n(i)} \), where \( \mathbbm{1}_{y}(i) \) indicates whether the node is the outcome node.
The number of HTE parents determines how many non-mediator parents (\( pa_n \)) interact with the treatment.
A value of zero means there is no heterogeneity.
So the final data generating equation is:
\begin{align}
    X_i & = C_{pa(i), i} X_{pa(i)}
    \nonumber \\
    & + \mathbbm{1}_{m}(i) \mathbbm{1}_{m_p} X_{pa_m(i)} X_{pa_n(i)}
    \nonumber \\
    & + \mathbbm{1}_{y}(i) X_{pa_m(i)} X_{pa_n(i)} + \epsilon_i
\end{align}
For each parameter, we randomly sample from all other parameter values \( 500 \) times to generate SCMs.
For each of the $11,000$ SCMs, we create one dataset of size \( 10,000 \).

\subsubsection{IHDP dataset}
The covariates in the Infant Health and Development Program (IHDP) dataset come from a randomized experiment studying the effects of specialist home visits on future cognitive test scores~\citep{hill-jcgs11}.
Outcomes can be simulated using the covariates.
We simulate 1,000 instances of simulated outcomes and average over all results\footnote{Outcomes are generated using the JustCause package: https://github.com/inovex/justcause}.

\subsubsection{Real-world datasets}

We additionally explore three real-world datasets.
The first dataset is a League of Legends (LoL) dataset~\citep{he-fdg21}.
LoL is one of the most popular multiplayer online battle arena (MOBA) games.
Players compete in teams to capture the opposing team's base.
Riot Games, the developer of LoL, regularly updates the game's software through \emph{patches}.
Here, we use patches as the treatment variable and measure how it affects number of \emph{kills} on a team-level.

\emph{Cannabis} is a dataset consisting of tweets about the e-cigarette Juul and cannabis-related topics~\citep{adhikari-icwsm21}.
English tweets related to the e-cigarette Juul are collected from 2016 to 2018.
For the users who mention Juul, cannabis-related tweets are collected from 2014 to 2018.
In this dataset, stance towards Juul and cannabis in tweets are crowdsourced, and a model is trained for stance detection.
In this work, we are interested in how users' stance on Juul affects their stance on cannabis.
To do this, we find users who have Juul-related tweets occurring before a cannabis-related tweet and detect their stance on Juul and cannabis using the classifier from~\citep{adhikari-icwsm21}, where features are the word embeddings for tweets using BERT~\citep{devlin-nacl19}.
The treatment variable of interest is a user's stance on Juul, in favor (1) or against (0).
The outcome variable is their stance on cannabis, in favor (1) or against (0).
Since this dataset contains only one group, we bootstrap the estimators to compute standard deviation.

The third dataset is about a study of interventions that may improve performance among undergraduates at the University of Illinois at Urbana Champaign (UIUC): strengthening human adaptive reasoning and problem solving (SHARPS).
In this dataset, students were given a battery of cognitive, physical, and diet evaluations.
In subsequent weeks, they were given one of four interventions: no training (active control), fitness training, fitness and cognitive training, and fitness, cognitive, and mindfulness training.
Since most HTE estimators can only handle binary treatment variables, we measure the effect of fitness training compared to active training.
The outcome variable we are interested in is the change in the number of correct answers on a digit-symbol substitution task.
In order to compute standard deviation, we bootstrap estimators.

\subsection{Evaluation}

When the ground truth is known, a common metric is the \emph{precision in estimation of heterogeneous effects} (PEHE):
\begin{equation}
    PEHE = \sqrt{\frac{1}{n} \sum_i^n \Big( (Y_i(1) - Y_i(0)) - E[Y_i(1) - Y_i(0) \mid \mathbf{X}_i = \mathbf{x}_i] \Big)^2},
\end{equation}
where the expectation is the prediction of an HTE estimator.
This is the root mean squared error (RMSE) of the true effect and the predicted effect.

In practice, we may be more interested in how well an HTE estimator can predict whether an individual will benefit from treatment or not.
For synthetic datasets, we compute an \emph{accuracy of decisions}.
For an individual, $ i $, let the true individual effect be $ \tau_i $, and the predicted effect be $ \hat{\tau}_i $.
Let the indicator $ \mathbbm{1}\big[ sign(\tau_i) = sign(\hat{\tau}_i) \big] $ determine if the sign of the predicted effect matches the sign of the true effect.
The accuracy is defined as:
\begin{equation}
    Acc = \frac{1}{N} \sum_{i}^{N} \mathbbm{1}[sign(\tau_i) = sign(\hat{\tau}_i)].
\end{equation}

If the ground truth is unknown, several metrics for heuristically evaluating HTE estimators have been proposed~\citep{saito-cv19,schuler-model18}.
Therefore, we additionally report the $ \tau  $-risk for evaluation:
\begin{equation}
    \tau\text{-risk} = \frac{1}{n} = \sum_i^n \Big( (Y_i - \hat{m}(\mathbf{X}_i)) - (T_i - \hat{p}(\mathbf{X}_i)) \hat{\tau}(\mathbf{X_i}) \Big)^2,
\end{equation}

When comparing estimators across different parameter settings and SCMs, we also rank all methods based on their error for a particular SCM\@.
Then we average the rank of each method across SCMs and report on the average rank.
To study performance under specific settings, we also report on average ranks for subsets of parameter settings.

\input{05-01_synthetic.tex}

\input{05-02_semi.tex}

\input{05-03_real.tex}

\section{Discussion}

\section{Summary and discussion}\label{sec:conclusion}

In this paper, we provide an introduction to data-driven heterogeneous treatment effect estimation and an overview of data-driven estimators that have been developed.
We broadly categorized data-driven estimators into counterfactual prediction methods and effect estimation methods.
To explore how these data-driven methods perform against each other, we explored their use in synthetic, semi-synthetic, and real-world datasets.
In general, we found that the performance of HTE estimators varied depending on the datasets, and there is no straightforward optimal method to use in all scenarios.

A current limitation of data-driven HTE estimators is the assumption of strong ignorability.
Our experiments show that without knowledge of the underlying model, the error of HTE estimations will increase, as shown in synthetic experiments.
This weakness can be more pronounced in real-world datasets, and addressing how to use the underlying model in HTE estimation is an essential avenue for future research.
Besides the assumption of strong ignorability, selection bias~\citep{heckman-jes79,cortes-alt08,bareinboim-pnas16} in the data can impact the accuracy of estimation. Some work has been done to address accounting for selection bias in causal models~\citep{bareinboim-pnas16}, which further emphasizes the importance of using causal models with data-driven estimation.
In addition, understanding the causal model, evaluation of current models currently has a large gap due to a lack of benchmark datasets.
A recent paper points out this gap and suggests alternatives for metrics and benchmark datasets~\citep{curth-neurips21} and more research in this area is an important future direction.
Finally, as we have explored in our experiments, heuristic evaluation metrics are not always reliable.
Understanding when these metrics can be used and how to improve them is another fruitful avenue of research.

\section{Discussion}

\section{Summary and discussion}\label{sec:conclusion}

In this paper, we provide an introduction to data-driven heterogeneous treatment effect estimation and an overview of data-driven estimators that have been developed.
We broadly categorized data-driven estimators into counterfactual prediction methods and effect estimation methods.
To explore how these data-driven methods perform against each other, we explored their use in synthetic, semi-synthetic, and real-world datasets.
In general, we found that the performance of HTE estimators varied depending on the datasets, and there is no straightforward optimal method to use in all scenarios.

A current limitation of data-driven HTE estimators is the assumption of strong ignorability.
Our experiments show that without knowledge of the underlying model, the error of HTE estimations will increase, as shown in synthetic experiments.
This weakness can be more pronounced in real-world datasets, and addressing how to use the underlying model in HTE estimation is an essential avenue for future research.
Besides the assumption of strong ignorability, selection bias~\citep{heckman-jes79,cortes-alt08,bareinboim-pnas16} in the data can impact the accuracy of estimation. Some work has been done to address accounting for selection bias in causal models~\citep{bareinboim-pnas16}, which further emphasizes the importance of using causal models with data-driven estimation.
In addition, understanding the causal model, evaluation of current models currently has a large gap due to a lack of benchmark datasets.
A recent paper points out this gap and suggests alternatives for metrics and benchmark datasets~\citep{curth-neurips21} and more research in this area is an important future direction.
Finally, as we have explored in our experiments, heuristic evaluation metrics are not always reliable.
Understanding when these metrics can be used and how to improve them is another fruitful avenue of research.

\subsection{Real-world datasets}

In real-world datasets, ground truth effects are not available.
Instead, we opt to compare HTE estimators based on $ \tau  $-risk as a proxy for methods perform in different datasets.
The average rank for each method in terms of $\tau$-risk for the LoL dataset is shown in Table~\ref{tab:ihdp_rankings}.
In this dataset, Causal Forest outranks other learners by a significant margin, followed by XLearner.
The rank in terms of $\tau$-risk for the Cannabis dataset is shown in Table~\ref{tab:ihdp_rankings}.
In the Cannabis dataset, XLearner achieves the highest $ \tau $-risk rank consistently.
Finally, the $\tau$-risk for the SHARPS dataset is shown in the last column of Table~\ref{tab:ihdp_rankings}.
In the SHARPS dataset, SLearner achieves the best rank consistently.
An interesting observation from this dataset is that the average difference in outcomes between treated and non-treated individuals is a small number (close to zero).
As pointed out by K{\"{u}}nzel et al., SLearner has a tendency to shrink estimates toward zero~\citep{kunzel-pnas19}, which may play a role in the $\tau$-risk rank in this dataset.
\section{Structural causal models for HTE estimation}\label{sec:structural_causal_models}

The data-driven estimation methods for counterfactual prediction and effect estimation we discussed so far do not explicitly capture how relationships between variables can affect the estimation of HTEs.
Instead, they typically rely on the strong ignorability assumption and use all features available in the data for estimation.
This can be problematic since including some variables could lead to incorrect estimation, which we discuss in more detail in Section \ref{sec:hte_estimation_with_scms}. To deal with this issue, it is helpful to know the cause-effect relationships between variables.
\emph{Structural causal models} (SCMs) are graphical models that capture these relationships and can help practitioners decide which variables should and should not be used for building HTE estimators.
One main drawback to SCMs is that they are rarely known in practice because the underlying causal mechanism that generated given data is unknown.
When the SCM is unknown, a causal structure learning algorithm can be used to learn a skeleton of the true model from the data itself.
In this section, we provide a brief introduction to SCMs and how HTEs are represented in them.

\subsection{Structural causal models}\label{sec:scm_background}

A structural causal model (SCM), $M$, consists of two sets of variables, $\mathbf{U}$ and $\mathbf{V}$, and a set of structural equations, $\mathbf{F}$, describing how values are assigned to each \emph{endogenous} variable $V^{(i)} \in \mathbf{V}$ based on the values of $\mathbf{v}$ and $\mathbf{u}$~\citep{pearl-book09}: \( v^{(i)} = f^{(i)}(\mathbf{v}, \mathbf{u}) \).
A causal graph \( G \) is a directed acyclic graph that captures the causal relationships among the variables.
Each variable, $V^{(i)}$, is represented by a node in the causal graph, an edge from $V^{(i)}$ to $V^{(j)}$ signifies that $V^{(i)}$ is a cause of $V^{(j)}$.
The edge relationship is a parent-child relationship, where the parent is a cause of the child.
The variables in $\mathbf{U}$ are considered \emph{exogenous}, and the causal graph may omit representing $\mathbf{U}$ explicitly.
Figure~\ref{fig:example_scm} shows an example of a causal graph with three exogenous variables, $ U_X, U_T, U_Y $ and three endogenous variables, $ X, T, Y $.

\begin{figure}
    \centering
    \newcommand{0.25}{0.5}

    \begin{subfigure}[b]{\columnwidth}
        \centering
        \begin{tikzpicture}
            \node (t) at (-0.75, 0) [label=below:$T$,point];
            \node (x) at (0.0, 0.75) [label=left:\(X\),point];
            \node(y) at (0.75, 0.0) [label=below:$Y$,point];
            \node (ux) at (0.0, 1.5) [label=above:\(U_X\),point];
            \node (ut) at (-1.5, 0.75) [label=above:\(U_T\),point];
            \node (uy) at (1.5, 0.75) [label=above:\(U_Y\),point];

            \path (x) edge (t);
            \path (x) edge (y);
            \path (t) edge (y);
            \path (ux) edge (x);
            \path (uy) edge (y);
            \path (ut) edge (t);
        \end{tikzpicture}
    \end{subfigure}

    \caption{An example causal model graph.
    }\label{fig:example_scm}
\end{figure}

An example set of associated structural equations is shown in eq~\eqref{eq:see1}-\eqref{eq:see3}.
\begin{gather}
    X = U_X \label{eq:see1} \\
    T = Bern(X) + U_T \label{eq:see2}  \\
    Y = \alpha \cdot T + \beta \cdot X + U_Y \label{eq:see3}
\end{gather}
In the graph, each endogenous variable has a directed edge pointing to it from an exogenous variable (e.g., $ U_X \rightarrow X $ or $ U_X $ is a parent of $ X $).
This is represented in the structural equations as being part of the structural equation (e.g., $ U_X $ in eq~\eqref{eq:see1}).
It is assumed that every endogenous variable has a parent in $ \mathbf{U} $.
We denote the parents and children of a variable, $ X $, in the graph as $ pa(X) $ and $ ch(X) $, respectively.
Ancestors of $ X $, $ an(X) $, are variables that can reach $ X $ through a directed path, and descendants of $ X $, $ de(X) $, are variables that can be reached from $ X $ through directed paths out of $ X $.

\subsubsection{Identifying effects using SCMs}\label{sec:identifying_scms}

When the SCM underlying the dataset is known, its graph $ G $ supports the use of graphical criteria for variable selection for unbiased causal effect estimation, known as adjustment set selection~\citep{shpitser-uai06,shpitser-uai10,van-uai14,perkovic-uai15,henckel-arxiv19}.
These graphical criteria rely on the concept of $d$-separation.
\begin{definition}
    ($d$-separation) A set of $ \mathbf{Z} $ nodes is said to ``block'' a path $ p $ if either
    \begin{enumerate}
        \item[$(i)$] $ p $ contains at least one arrow-emitting node that is in $ \mathbf{Z} $ or
        \item[$(ii)$] $p$ contains at least one collision node that is outside $ \mathbf{Z} $ and has no descendant in $ \mathbf{Z} $.
    \end{enumerate}
\end{definition}
If $ \mathbf{Z} $ blocks all paths from a set of nodes $ \mathbf{X} $ to a set of nodes $ \mathbf{Y} $, it is said to ``d-separate'' $ \mathbf{X} $ and $ \mathbf{Y} $, and the conditional independence, $ \mathbf{X} \perp \mathbf{Y} \mid \mathbf{Z} $, holds~\citep{bareinboim-pnas16}.
One graphical criteria for adjustment set selection is the backdoor criterion~\citep{pearl-book09}.
The backdoor criterion aims to find a set that \emph{blocks} the paths from the treatment variable to the outcome variable in the graph, starting with an arrow into the treatment variable.
Essentially, the backdoor criterion blocks all the paths going backward out of treatment.
For example, in Figure~\ref{fig:example_scm}, there is one backdoor path, namely $ T \leftarrow X \rightarrow Y $, and so the adjustment set $ \{X\} $ sufficiently satisfies the backdoor criterion.
Several adjustment criteria have been developed, such as the adjustment criterion~\citep{shpitser-uai06,van-uai14}, generalized backdoor criterion~\citep{maathuis-annals15}, generalized adjustment criterion~\citep{perkovic-uai15}, and optimal adjustment set (O-set)~\citep{henckel-arxiv19}.
These adjustment sets defined by the various adjustment criteria can be found through \emph{identification} algorithms such as the ID algorithm~\citep{shpitser-jmlr08}, the algorithmic framework by Van der Zander et al.~\citep{van-uai14} or the pruning procedure by Henckel et al.~\citep{henckel-arxiv19}.

When an adjustment set is discovered, unbiased causal effects can be estimated by conditioning on the adjustment set.
The adjustment set can be used instead of all the features for data-driven estimation.
We discuss data-driven HTE estimation with SCMs briefly in Section~\ref{sec:hte_estimation_with_scms}.

\subsection{Heterogeneous treatment effects in SCMs}

Shpitser and Pearl provide a nonparametric estimand for HTEs for some feature assignment $C=c$, called $c$-specific effects~\citep{shpitser-uai06,pearl-smr17}.
One issue with estimating $c$-specific effects is that it requires knowledge of the SCM and which variables, C, exhibit heterogeneous effects.
When the SCM is known, the O-set~\citep{henckel-arxiv19} can be used on the learned structure since it implicitly includes all potential heterogeneity-inducing variables, while some identification algorithms may miss them. When the SCM is unknown, the identifiability of \( c \)-specific effects cannot be established.
Moreover, when the set of heterogeneity-inducing variables $\mathbf{C}$ is unknown, then all combinations of features and feature values would need to be checked for heterogeneity in effect, which is exponential in the number of valid features.
Circumventing this exponential search is one of the main reasons data-driven methods are popular in practice.

Shpitser and Pearl provide a nonparametric estimand for HTEs with respect to some feature assignment $C=c$, called $c$-specific effects~\citep{shpitser-uai06,pearl-smr17}.
There are three cases when $c$-specific effects are \emph{identifiable}, i.e., there exists an unbiased estimation for $\tau(\mathbf{c}) = E[Y(1) - Y(0) | \mathbf{C}=\mathbf{c}]$.
The first case is when $\mathbf{C}$ is observed and admissible which means that it satisfies the backdoor criterion in the causal graph~\citep{pearl-book09}.
A set of variables $\mathbf{Z}$ satisfies the backdoor criterion relative to $(T,Y)$ if no node in $\mathbf{Z}$ is a descendant of $T$ and $\mathbf{Z}$ blocks every path between $T$ and $Y$ that contains an arrow into $T$.
For example, $X$ is the only variable that is admissible in all SCMs in Figure~\ref{fig:fail_examples}, since it blocks the only backdoor path.

The second case is when $\mathbf{C}$ is part of an observed admissible set.
In this case, $\mathbf{C}$ itself is not admissible, but we can observe a set of features $\mathbf{S}$ such that $\mathbf{S} \cup \mathbf{C}$ is admissible.
In Figure~\ref{fig:marloes_scm}, $\{F\}$ could induce heterogeneity but is not admissible on its own. Adding $F$ to $\{X\}$ still satisfies the backdoor criterion making $\{X, F\}$ also admissible.
The third case for identifiable $c$-specific effects is when $\mathbf{C}$ is not part of an admissible set.
The first two cases which rely on \emph{covariate adjustment}~\citep{perkovic-jmlr18} and identifying heterogeneity-inducing variables are the basis of our work.
Estimating heterogeneous effects in the third case is beyond the scope of this paper and is reserved for future work.

One issue with estimating $c$-specific effects is that it requires knowledge of the SCM and which variables exhibit heterogeneous effects.
If the SCM is unknown, the identifiability of \( c \)-specific effects cannot be established.
If $\mathbf{C}$ is unknown, then all combinations of features and feature values would need to be checked for heterogeneity in effect, which is exponential in the number of valid features~\citep{pearl-smr17}.
Circumventing this exponential search is one of the main reasons why data-driven methods are popular in practice.

\subsubsection{Detecting latent heterogeneity}

In some cases, the observed covariates that induce heterogeneity, $C$, are not known, but there is unknown heterogeneity in effects present in the data~\citep{pearl-smr17}.
Pearl discusses how heterogeneity can be detected from statistical averages alone~\citep{pearl-smr17} when the effect of treatment on the treated ($ETT$) and the effect of treatment on the untreated ($ETU$) are identifiable.
The three cases where $ETT$ can be identified are when:
\begin{itemize}
    \item The treatment is binary, and $E[Y_1]$ and $E[Y_0]$ are identifiable by some method
    \item The treatment is arbitrary, and $E[Y_x]$ is identifiable (for all $x$) by adjustment for an admissible set of covariates.
    \item $ATE$ is identified through mediating instruments.
\end{itemize}
The data-driven HTE estimation methods discussed in this article do not focus on the detection of latent heterogeneity, but using these cases when latent heterogeneity can be detected with data-driven estimators is an interesting future direction.

\subsubsection{Representation of HTEs in structural equations}\label{sec:htes_scms_background}

In SCMs, heterogeneous effects are represented by an \emph{interaction} or \emph{effect modification} through ancestors and parents of a variable.
We will use linear structural equations for illustration.
For heterogeneous effects, this interaction is specifically with the treatment variable.
Consider the causal model in Figure~\ref{fig:example_scm}, where the treatment variable of interest is $ T $ and the outcome of interest is $ Y $.
There is a potential interaction between $ X $ and $ T $.
Suppose $ Y $ is generated by the following structural equation:
\begin{equation}
    Y = \alpha \cdot T + \beta \cdot X + \gamma \cdot T \cdot X + U_Y.
\end{equation}
Here, there is a fixed effect \( \alpha \) of the treatment \( T \) and a heterogeneous effect \( \gamma \) due to the interaction between \( T \) and \( X \).
Since the effect of treatment varies with another variable, $ X $, there is a heterogeneous effect.
Heterogeneity through interaction can occur with any parent of the outcome or any parent of \emph{mediating} variables~\citep{vanderweele-epid07}.

\begin{figure}
    \centering
    \begin{tikzpicture}
        \node (t) at (-1.5, 0) [label={below left}:$T$,point];
        \node (x) at (-1.5, 1.5) [label={above}:$X$,point];

        \node (c) at (0.0, 1.5) [label={above}:$C$,point];
        \node (d) at (0.0, 0) [label={below}:$D$,point];

        \node (f) at (1.5, 1.5) [label={above}:$F$,point];
        \node (y) at (1.5, 0.0) [label={below}:$Y$,point];

        \path (x) edge (t);

        \path (x) edge (y);
        \path(t)edge(d);
        \path(c)edge(d);
        \path (d) edge (y);

        \path (f) edge (y);
    \end{tikzpicture}
    \caption{Causal model with multiple possible interactions that cause heterogeneous treatment effects.}\label{fig:multi_interaction}
\end{figure}

Consider another example causal model, shown in Figure~\ref{fig:multi_interaction}.
Here, there is a \emph{mediating} variable, $ D $, a variable that occurs on a path from the treatment, $ T $, to the outcome, $ Y $.
Interaction can occur with $ X, C$, or $ F $.
For example, an interaction with $ X $ and $ F $ can be shown as:
\begin{equation}
    Y = \alpha \cdot D \cdot X + \beta \cdot D \cdot F + U_Y.
\end{equation}
Since the treatment is mediated through $ D $ and $ X $ and $ F $ interact with $ D $, there will be a heterogeneous effect.
Since $ D $ is a mediator, if $ C $ interacts with $ T $ in the structural equation, then there will also be a heterogeneous effect, such as:
\begin{gather}
    D = \alpha \cdot T + \beta \cdot T \cdot C + U_D, \\
    Y = \gamma \cdot D + U_Y.
\end{gather}

\subsection{Data-driven HTE estimation using Structural Causal Models}\label{sec:hte_estimation_with_scms}

\begin{figure}
    \centering
    \newcommand{0.25}{0.3}
    \resizebox*{0.25\columnwidth}{!}{
        \begin{subfigure}[b]{0.325\columnwidth}
            \centering
            \begin{tikzpicture}
                \node (t) at (-0.75, 0) [label=left:$T$,point];
                \node (x) at (0.0, 0.75) [label=above:\(X\),point];
                \node(y) at (0.75, 0.0) [label=right:$Y$,point];
                \node(c) at (0.0, -0.75) [label=right:$L$,point];

                \path (x) edge (t);
                \path (x) edge (y);
                \path (t) edge (y);
                \path (t) edge (c);
                \path (y) edge (c);
            \end{tikzpicture}
            \caption{}\label{fig:collider_scm}
        \end{subfigure}
    }
    \resizebox*{0.25\columnwidth}{!}{
        \begin{subfigure}[b]{0.325\columnwidth}
            \centering
            \begin{tikzpicture}
                \node (t) at (-0.75, 0) [label=left:$T$,point];
                \node (m) at (0.0, 0) [label=below:$M$,point];
                \node (x) at (0.0, 0.75) [label=above:\(X\),point];
                \node(y) at (0.75, 0.0) [label=right:$Y$,point];

                \path (x) edge (t);
                \path (x) edge (y);
                \path (t) edge (m);
                \path (m) edge (y);
            \end{tikzpicture}
            \caption{}\label{fig:mediator_scm}

        \end{subfigure}
    }
    \resizebox*{0.25\columnwidth}{!}{
        \begin{subfigure}[b]{0.325\columnwidth}
            \centering
            \begin{tikzpicture}
                \node (a) at (-1.5, 0) [label=left:$A$,point];
                \node (t) at (-0.75, 0) [label={below left}:$T$,point];
                \node (b) at (-0.75, -0.75) [label={below}:$B$,point];
                \node (x) at (-0.75, 0.75) [label={above}:$X$,point];

                \node (c) at (0.0, 0.75) [label={above}:$C$,point];
                \node (d) at (0.0, 0) [label={below left}:$D$,point];
                \node (e) at (0.0, -0.75) [label={below}:$E$,point];

                \node (f) at (0.75, 0.75) [label={above}:$F$,point];
                \node (y) at (0.75, 0.0) [label={right}:$Y$,point];
                \node (g) at (0.75, -0.75) [label={below}:$G$,point];

                \path (a) edge (t);
                \path (x) edge (t);
                \path (t) edge (b);

                \path (x) edge (y);
                \path(t)edge(d);
                \path(c)edge(d);
                \path(d)edge(e);
                \path (d) edge (y);

                \path (f) edge (y);
                \path (y) edge (g);
            \end{tikzpicture}
            \caption{}\label{fig:marloes_scm}
        \end{subfigure}
    }
    \caption{Causal models where data-driven HTE estimation methods may perform poorly if the structure is unknown.
    Some variables (e.g., $L$, $M$, and $E$) are not valid adjustment variables for estimating the effect of $T$ on $Y$.
    }\label{fig:fail_examples}
\end{figure}

As previously mentioned, data-driven methods typically use all features available in the data for estimation.
However, not all variables in the data should be included in HTE estimation, and including incorrect variables leads to incorrect effect estimation~\citep{pearl-book09}.
Incorrect variables include descendants of treatment or instrumental variables~\citep{pearl-book09,myers-aje11} (e.g., $L, M, A, B, D, E, G$ in Figure~\ref{fig:fail_examples}).
If the underlying causal model is known, then an adjustment set can be discovered using an adjustment criterion.
Then, unbiased effects can be estimated by using the adjustment set as the features for estimation.

\subsubsection{Causal structure learning}\label{sec:structure_learning}

In many real-world scenarios, the underlying causal model is unknown, and therefore, there is \emph{structure uncertainty} of the causal mechanisms.
Therefore, ignoring the underlying causal models can lead to incorrect causal effect estimations.

To use the causal model under structure uncertainty, one approach is to use causal structure learning~\citep{heinze-stat18,vowels-arxiv21}.
Many traditional causal structure learning methods can be broadly classified into constraint-based~\citep{spirtes-book00,colombo-jstor12}, score-based~\citep{chickering-jmlr02,hauser-jmlr12,nandy-annals18}, a hybrid of the two~\citep{tsamardinos-ml06}, or based on asymmetries~\citep{shimizu-jmlr06}.
These methods find a partial graph or skeleton from the data, with some unoriented edges.
With the output from structure learning methods, an adjustment set can be applied to select variables for estimation~\citep{pearl-book09,henckel-arxiv19}, and these variables are used to train an HTE estimator.

One well-known constraint-based structure learning algorithm is the PC algorithm~\citep{spirtes-book00}.
PC conducts conditional independence tests to learn about the structure of the underlying model. It learns a completed partially directed acyclic graph (CPDAG)~\citep{chickering-jmlr02}, which is an equivalence class of directed acyclic graphs (DAG) that have uncertain edges.
PC has three main steps: (1) determining the skeleton, (2) determining v-structures, and (3), determining further edge orientations.
PC learns the skeleton by starting with a complete undirected graph.
Using a series of conditional independence tests for adjacent nodes, PC removes an edge if a conditional independence is found.
In the second step, all edges in the graph are replaced with $\circ$--$\circ$ edges.
PC considers all unshielded triples, triples $i \circ $--$\circ j$--$\circ k$ where $i$ and $k$ are not adjacent and determines whether a v-structure should exist or not.
In the final step, additional orientation rules are applied to orient as many unoriented edges as possible.
The fast causal inference (FCI) algorithm is a modification of the PC algorithm that allows hidden variables~\citep{spirtes-book00}.

A classic score-based algorithm is greedy equivalence search (GES) which also learns a CPDAG.
Score-based algorithms typically use a likelihood score of a DAG against the data, and search for a DAG or CPDAG that yields the optimal score.
GES learns the CPDAG via a greedy forward search of edge additions, followed by a backward phase with edge deletions.
A hybrid approach based on FCI and GES has also been developed, called greedy fast causal inference (GFCI)~\citep{ogarrio-pgm16}.
Heinze-Deml et al.\ provide a survey on structure discovery methods, specifically constraint-based and score-based approaches~\citep{heinze-stat18}.

Besides constraint-based approaches and score-based approaches, continuous optimization-based methods for structure discovery have been developed due to increased compute that make learning from high-dimensional datasets feasible~\citep{ma-sr14,zheng-neurips18,kyono-neurips20,vowels-arxiv21}.
Continuous optimization-based methods recast the graph-search problem into a continuous optimization problem~\citep{vowels-arxiv21}.
Vowels et al.\ provide a survey on continuous optimization-based approaches~\citep{vowels-arxiv21}.
With the learned causal graph, an \emph{identification} algorithm can be used to find an adjustment set.
Then, this adjustment set can be used as features in an HTE estimator instead of the full feature set.

Although causal structure learning can address some problems of data-driven HTE estimation, there are some drawbacks.
First, learning the full structure can be computationally expensive for high-dimensional datasets~\citep{chickering-jmlr04,heinze-stat18,vowels-arxiv21}.
Second, causal structure learning generally cannot learn the full causal model, so learned models will have undirected edges with an ambiguous cause-effect relationship.
In addition, adjustment criteria have been developed for total causal effect estimation without heterogeneity in mind.
Therefore, adjustment sets can potentially miss heterogeneity-inducing variables, leading to incorrect HTE estimation.
In our experiments, we use randomly generated SCMs and explore how data-driven HTE estimators perform under different types of SCMs.
\subsection{Effect estimation}\label{sec:effect_estimation}

In contrast to methods that perform counterfactual prediction, effect estimation methods estimate the treatment effect directly.
They include methods that utilize supervised learning~\citep{athey-stat15,kunzel-pnas19,yoon-iclr18} and methods that have tailored objectives for learning HTEs~\citep{athey-pnas16,tran-aaai19,wager-jasa18}.

\subsubsection{Transformed outcomes}

In the transformed outcomes approach, the outcome variable is transformed such that a regression or classification method can be used to predict the transformed outcome.
Any off-the-shelf estimator can be directly applied since the CATE is estimated directly from the transformed outcome.

Under the assumption of unconfoundedness, the transformed outcome is~\citep{athey-stat15}:
\begin{equation}
    Y_i^* = Y_i \cdot \frac{T_i - e(\mathbf{X}_i)}{e(\mathbf{X}_i) \cdot (1 - e(\mathbf{X}_i) )},
\end{equation}
where \( Y_i^* \) is the transformed outcome and \( e(\mathbf{x}) = \Pr(T_i = 1 \mid \mathbf{X}_i = \mathbf{x}) \) is the propensity score~\citep{rosenbaum-bio83}.
In practice, the propensity score is estimated from data using a supervised estimator such as logistic regression~\citep{setoguchi-pharm08}.
In this setting, the CATE is estimated as the estimated conditional mean of the transformed outcome.
Using the transformed outcome, the conditional mean of the predicted outcome is equal to the CATE:
\begin{equation}
    \hat{\tau}_{to}(\mathbf{x}) = \mathbb{E}[Y^{*} \mid \mathbf{X} = \mathbf{x}].
\end{equation}
Transformed outcomes has been used with regression trees~\citep{athey-stat15} and for building transformed outcome forests~\citep{powers-statmed18}.

\subsubsection{Meta-learning}

Meta-learners or meta-algorithms build upon supervised learning methods, called base learners, to estimate CATE\@.
Two special cases of meta-algorithms are the single-model and two-model approaches, as they take supervised base learners and utilize their outputs to estimate CATE\@.
Xie et al.\ propose three algorithms based on propensity scores: Stratification-Multilevel (SM), Matching-Smoothing (MS), and Smoothing-Differencing (SD)~\citep{xie-soc12}.
The common trend of these algorithms is to estimate propensity scores using a supervised learning method and create strata based on propensity scores and covariates.
In SM, stratum-specific treatment effects are estimated, and a trend across strata is obtained through variance-weighted least squares regression.
For MS, treated units are matched to control units with a matching algorithm to estimate the true effect, and a nonparametric model is applied to the matched effects.
Finally, for SD, nonparametric regression is used to fit the outcome on the propensity score.
To estimate heterogeneous effects, take the difference between the nonparametric regression line between treated and untreated at different levels of the propensity score.

Powers et al.\ propose Pollinated Transformed Outcome (PTO) forests~\citep{powers-statmed18}.
First, they use a transformed outcome approach and build a random forest.
Then, predictions in trees are replaced with the estimate of $ \hat{\mu}(T=1, \mathbf{x}) $ and $ \hat{\mu}(T=0, \mathbf{x}) $ as in the single model approach.
Finally, an additional random forest is build to predict $ \hat{\tau} $.

K{\"{u}}nzel et al.\ propose X-Learner, which aims to improve upon the two-model approach~\citep{kunzel-pnas19}.
The goal of the X-Learner is to use information from the control group to learn better estimates for the treatment group and vice versa.
X-Learner consists of three steps.
The first step is to learn two supervised estimators $ \hat{\mu}_0 $ and $ \hat{\mu}_1 $.
This is the same as the two-model approach.
Step two involves imputing the treatment effects from predictions from the opposite group by taking the difference between factual and predicted counterfactual outcomes:
\begin{gather}
    \tilde{D}_i^1 = Y_i^{1} - \hat{\mu}\big(\mathbf{X}_i^{1} \big), \\
    \tilde{D}_i^0 = \hat{\mu}\big( \mathbf{X}_i^{0} \big) - Y_i^{0},
\end{gather}
where $ \tilde{D}_i^1 $ and $ \tilde{D}_i^0 $ is the imputed treatment effect for unit $ i $ if they were treated or untreated, respectively.
Then, two estimators $ \hat{\tau}_0(\mathbf{x}) $ and $ \hat{\tau}(\mathbf{x}) $, are trained on the imputed datasets, $ \tilde{D}_i^1 $ and $ \tilde{D}_i^0 $, respectively.
Finally, the CATE is estimated by a weighted average of the two estimators:
\begin{equation}
    \hat{\tau}(\mathbf{x}) = g(\mathbf{x}) \hat{\tau}_0(\mathbf{x}) + (1 - g(\mathbf{x})) \hat{\tau}_1(\mathbf{x}),
\end{equation}
where $ g \in [0,1]$ is a weight function.
K{\"{u}}nzel et al.\ suggest using an estimate of the propensity score for $ g $, so that $ g = \hat{e}(\mathbf{x}) $~\citep{kunzel-pnas19}.

Kennedy proposes a two-stage doubly robust method for HTE estimation called Double Robust Learner (DR-Learner)~\citep{kennedy-arxiv20}.
The training data is randomly separated into three sets: $ D_{1a}^k, D_{1b}^k, D_{2}^k $ each containing $ k $ observations.
In the first step, two estimators are trained on $ D_{1a}^k $ and $ D_{1b}^k $, called \emph{nuisance} training.
An estimator of the propensity scores, $ \hat{e} $, is constructed on $ D_{1a} $ and estimates of $ (\hat{\mu}_0 $ and $ \hat{\mu}_1) $ are constructed by performing prediction $ D_{1b} $.
Then, a pseudo-outcome is constructed as:
\begin{equation}
    \hat{\varphi}(\mathbf{X}, Y, T) = \frac{T - \hat{e}(\mathbf{X})}{\hat{e}(\mathbf{X}) (1 - \hat{e}(\mathbf{X})) } \Big( Y - \hat{\mu}_T(\mathbf{X}) \Big) + \hat{\mu}_1(\mathbf{X}) + \hat{\mu}_0(\mathbf{X}).
\end{equation}
Finally, regressing the pseudo-outcome on $ \mathbf{X} $ in the test sample $ D_{2}^k $ obtains estimations of the CATE\@:
\begin{equation}
    \hat{\tau}_{dr}(\mathbf{x}) = \mathbb{E}[\hat{\varphi}(\mathbf{X}, Y, T) \mid \mathbf{X}=\mathbf{x}].
\end{equation}
It is recommended to cross-fit with the three sets and use the average of the resulting CATE estimators.

\subsubsection{Tree and forest-based methods}

Tree and forest-based methods are tailored for direct HTE estimation without predicting potential outcomes.
Su et al.\ propose a tree-based method inspired by the CART algorithm for classification~\citep{brieman-cart84, su-jmlr09}.
They find binary covariate splits based on a t-test of differences between two potential children of a node.
Let \( \{\mu^L_1, \bar{y}^L_1, s_2^L, n_1\} \) be the population mean, sample mean, sample variance, and the sample size for the treated group in the left child node \(L\), and \( \{\mu^L_0, \bar{y}^L_0, s_2^L, n_2\} \) similarly for the control group.
Su et al.\ propose a measure for assessing the interaction between the split, \( s \) and the treatment by:
\begin{equation}
    g(s) = \frac{(\bar{y}^L_1 - \bar{y}^L_0) - (\bar{y}^R_1 - \bar{y}^R_0)}{\hat{\sigma} \sqrt{1/n_1 + 1/n_2 + 1/n_3 + 1/n_4}},
\end{equation}
where \( \hat{\sigma} = \sum_{i=1}^4 w_i s_i^2 \) is a pooled estimator of the constant variance, and \( w_i = (n_i - 1) / \sum_{j=1}^4 (n_j - 1) \).
For a split \( s \), \( F_t = g^2(s) \) converges to a \( \chi^2(1) \) distribution, which is used as the splitting criterion~\citep{su-jmlr09}.

Athey and Imbens introduce two splitting criteria for building trees based on mean squared error for predictions: \emph{adaptive} and \emph{honest}~\citep{athey-pnas16}.
The main goal of these trees is to maximize the observed treatment effects in subgroups to discover heterogeneity.
In a node \( \ell \), define the conditional mean outcomes under treated and control as:
\begin{gather}
    \hat{\mu}(t=1, \ell) = E[Y_i(t=1) \mid X_i \in \ell],
    \\
    \hat{\mu}(t=0, \ell) = E[Y_i(t=0) \mid X_i \in \ell].
\end{gather}
The average causal effect in a node \( \ell \) is defined as the difference in conditional means:
\begin{gather}
    \hat{\tau}(\ell) = \hat{\mu}(1, \ell) - \hat{\mu}(0, \ell).
\end{gather}
Then, the adaptive \emph{partition measure} for node \( \ell \), which is used as the score for a particular partition or node, can be defined as:
\begin{gather}
    F_a(\ell) = N_\ell \cdot \hat{\tau}^2(\ell).
\end{gather}
To build a causal tree, this partition measure is maximized at each potential split:
\begin{gather}
    \max_\ell \; \sum_{\ell} F_a.
\end{gather}

In the \emph{honest} approach, the training data is separated into two parts: \emph{training} and \emph{estimation}.
The training set is used to find splits, but causal effects are estimated using data in the estimation set.
In addition, a penalty on high variance in outcomes in a node is introduced to the splitting criteria.
Define \( V^{(1)}(\ell), V^{(0)}(\ell) \) be the variance in outcomes in node \( \ell \) and \( p \) to be the treatment probability in the training data.
The honest splitting criteria is defined as:
\begin{gather}
    F_h(\ell) = F_a - \Bigg( \frac{1}{N^{tr}} + \frac{1}{N^{est}} \Bigg)  \Bigg( \frac{V^{(1)}(\ell)}{p} + \frac{V^{(0)}(\ell)}{1-p} \Bigg)
\end{gather}
This approach balances finding heterogeneity and lowering variance in leaf nodes.

Tran and Zheleva introduced the idea of using validation sets for building causal trees~\citep{tran-aaai19}.
They propose three variants: CT-L, CT-HL, and CT-HV\@.
In CT-L, a separate validation set is used to penalize splits with estimated effects that do not match with the training data.
If we define the estimated causal effect in node \( \ell \) on the training data to be \( \hat{\tau}(\ell^{tr}) \) and the estimated causal effect on the validation data to be \( \hat{\tau}(\ell^{val}) \), then a cost term \( C \) can be defined as:
\begin{equation}
    C(\ell) =  N_\ell^{val} \cdot | \hat{\tau}(\ell^{tr}) - \hat{\tau}(\ell^{val}) |,
\end{equation}
where \( N_\ell^{val} \) is the size of the validation set in node \( \ell \).
The cost term measures the error of the estimated effect of training and validation data.
Using the cost term, the partition measure for CT-L is defined as:
\begin{gather}
    F_l(\ell): \frac{(1 - \lambda) \cdot F_a - \lambda \cdot C(\ell) }{| N_\ell^{tr} - N_\ell^{val}| + 1},
\end{gather}
where \( \lambda \) is used to control the strength of the cost function.
The other two variants, CT-HL and CT-HV, incorporate the honest penalty into building the tree.
In CT-HL, the data is split into three sets: training, estimation, and validation, and the honest penalty and validation cost are added together.
In CT-HV, the data is only split into two sets: training and validation.
The difference is that the validation set is also used as the honest estimation set, so less data needs to be split from the training data.

In addition to the introduction of CT-L, Tran and Zheleva introduce the concept of \emph{triggers} for heterogeneous treatment effects~\citep{tran-aaai19}.
When the treatment is an ordinal or monotonously increasing continuous variable, the effect depends on the amount of treatment.
For example, a clinician may be interested in understanding the minimum number of days or amount of treatment (the trigger) that patients need to take a medication to be cured~\citep{tran-aaai19}.
The goal becomes finding the threshold that maximizes the expected effect for subpopulations.

Wager and Athey proposed a random forest algorithm called Causal Forest for estimating heterogeneous treatment effects~\citep{wager-jasa18}.
Causal Forest constructs a forest of honest Causal Trees that each provides a CATE estimation. Then, an average of the predicted CATE from all trees is used as the final prediction.
They showed that the estimation of Causal Forest is asymptotically Gaussian and unbiased for the true CATE if the base learner used is \emph{honest}~\citep{wager-jasa18}.
Athey et al.\ extended this idea to Generalized Random Forests~\citep{athey-annals19}.
Generalized random forests extend random forests into a method for estimating any quantity through local moment conditions and can be used for heterogeneous treatment effect estimation.

\subsubsection{Deep learning approaches}

Several other deep learning approaches have also been developed based on different architectures.
Louizos et al.\ consider a different perspective and motivation for studying hidden representations for HTE estimation~\citep{louizos-neurips17}.
Instead of considering the given features $ \mathbf{X} $ as confounders, they consider it as noisy views on the hidden confounder $ \mathbf{Z} $.
Their goal is to learn the representation of $ \mathbf{Z} $ for causal effect adjustment using the observed features $ \mathbf{X} $.
To do this, they use variational autoencoders to learn $ \mathbf{Z} $ for the treatment and control groups, and use the hidden representation for prediction of outcomes.

Yoon et al.\ propose Generative Adversarial Networks for Individualized Treatment Effects (GANITE)~\citep{yoon-iclr18}.
GANITE has two components, a counterfactual block which generates counterfactual predictions, and an ITE block that generates CATE predictions given the counterfactual prediction.
This work leverages Generative Adversarial Networks (GANs) in each block to create a better counterfactual dataset.
After training the counterfactual block, generated data is fed to the ITE block which uses a GAN to learn to generate HTEs.
\section{Introduction}\label{sec:intro}

Estimating the effect of a treatment on an outcome is a fundamental problem in many fields such as medicine~\citep{shalit-icml17,laber-bio15,lakkaraju-aistats17}, public policy~\citep{grimmer-pa17} and more~\citep{li-www10,ascarza-jmr18}.
For example, doctors might be interested in how a treatment, such as a drug, affects the recovery of patients~\citep{foster-statmed11}, economists may be interested in how a job training program affects employment prospectives~\citep{lalonde-aer86}, and advertisers may want to model the average effect an advertisement has on sales~\cite{lewis-retail08}.

However, individuals may react differently to the treatment of interest, and knowing only the average treatment effect in the population is insufficient.
For example, a drug may have adverse effects on some individuals but not others~\citep{shalit-icml17}, or a person's education and background may affect how much they benefit from job training~\citep{lalonde-aer86,pearl-smr17}.
Measuring the extent to which different individuals react differently to treatment is known as \emph{heterogeneous treatment effect (HTE) estimation}.

Traditionally, HTE estimation has been done through subgroup analysis~\citep{gail-bio95,bonetti-bio04}.
However, this can lead to cherry-picking since the practitioner is the one who identifies subgroups for estimating effects.
Recently,
there has been more focus on \emph{data-driven} estimation of heterogeneous treatment effects by letting the data identify which features are important for treatment effect estimation  using machine learning techniques~\citep{johansson-icml16,louizos-neurips17,shalit-icml17,tran-aaai19}.
A straightforward approach is to create interaction terms between all covariates and use them in a regression~\citep{athey-handbook17}.
However, this approach is limited to low-dimensional datasets and is not practical for large datasets with many features.
To deal with the high dimensionality of data, many supervised learning-based methods have been developed and adapted to the problem of HTE estimation~\citep{su-jmlr09,hill-jcgs11,imai-annals13,athey-stat15,johansson-icml16,athey-pnas16,shalit-icml17,tran-aaai19,kunzel-pnas19}.

HTE estimation methods can be broadly categorized into counterfactual prediction and effect estimation methods.
A counterfactual prediction model learns to predict each unit's potential outcomes when treated and not treated.
In turn, those predictions can be used to estimate a treatment effect~\citep{foster-statmed11,shalit-icml17,kunzel-pnas19}.
A heterogeneous effect estimation method learns a model that predicts the treatment effect directly~\citep{su-jmlr09,athey-pnas16,tran-aaai19}.
Effect estimation models can be learned through specialized objective functions~\citep{su-jmlr09,athey-pnas16} or in tandem with counterfactual prediction~\citep{louizos-neurips17,kunzel-pnas19}.
Typically, counterfactual prediction and effect estimation methods use all variables in the data for estimation relying on the strong ignorability assumption~\citep{montgomery-ajps18}.
However, as we discuss in Section~\ref{sec:structural_causal_models} and show in our experiments, including all variables can lead to errors in HTE estimation.

A third category of methods, which rely on structural causal models (SCMs)~\citep{pearl-book09}, provide criteria for the selection of variables that allow for unbiased causal effects. SCMs are graphical models that model explicitly the cause-effect relationships between variables.
The selection of variables through SCMs can help data-driven HTE estimation methods better estimate causal effects, as we show in our experiments.
However, the underlying causal mechanism of data and its SCM are rarely known in real-world settings. To circumvent this problem, one can use a causal structure learning algorithm~\citep{spirtes-book00} to learn a skeleton of the underlying causal model from the data itself.
This survey provides an overview of recently developed data-driven HTE estimation methods, and an overview of SCMs and how they can be used for HTE estimation.

The rest of the paper is organized as follows.
Section~\ref{sec:background} gives basic background knowledge.
Section~\ref{sec:survey} reviews counterfactual prediction and effect estimation methods.
Section~\ref{sec:structural_causal_models} provides an overview of data-driven HTE estimation with structural causal models.
Section~\ref{sec:experiments} empirically evaluates HTE estimators on several types of datasets.
Section~\ref{sec:conclusion} concludes the paper and discusses open problems.
\section{Experiments}\label{sec:experiments}

We evaluate HTE estimators on synthetic, semi-synthetic, and real-world experiments.
We test several prominent HTE estimators in our experiments.

\subsection{Experimental setup}\label{sec:setup}

In our experiments, we compare several HTE estimators.
We employ several methods from recent work:
Causal Trees~\citep{tran-aaai19}, Causal Forest~\citep{athey-annals19},
Meta-learners (S, T, XLearners)~\citep{kunzel-pnas19}, Doubly Robust learner (DRLearner)~\citep{kennedy-arxiv20}, Balancing Neural Network~\citep{johansson-icml16}, and TARNet~\citep{shalit-icml17}.
In addition, we also use a structure learning method to find the optimal adjustment set (O-set)~\citep{henckel-arxiv19}.
We opt to use a hybrid method, greedy fast causal inference (GFCI)~\citep{ogarrio-pgm16}, and abbreviate this method as full structure learning with O-set (FSL (O-set)).
For synthetic datasets, where the causal model is known, we use an Oracle method (Known Structure) that uses the causal model for variable selection using the \emph{optimal adjustment set (O-set)}~\citep{henckel-arxiv19} and a TLearner for estimation.
We use Logistic Regression and Linear Regression for propensity scores and estimating outcomes for meta-learners on synthetic datasets.
For semi-synthetic and real-world datasets, we use Random Forest for estimating outcomes.

\subsection{Datasets}
We utilize three types of datasets.
The first type of dataset is a synthetic dataset.
The second type of dataset is a semi-synthetic dataset, the Infant Health and Development Program (IHDP) dataset~\citep{hill-jcgs11}.
The final type of dataset are real-world datasets.
We use League of Legends, Cannabis, and SHARPS.

\subsubsection{Synthetic dataset}\label{sec:setup_synthetic}

We first evaluate HTE estimators on synthetic datasets generated by random structural causal models (SCMs).
As mentioned in Section~\ref{sec:hte_estimation_with_scms}, not all variables in the data should be included in HTE estimation.
A current gap in the literature is understanding how data-driven HTE estimators perform when variables should not be included due to the underlying structure, such as mediators or descendants of treatment.
Therefore, we seek to understand how estimators perform under differing structures in our synthetic experiments.
In this setting, we assume estimators do not have access to the underlying causal graph and instead use all features in the generated dataset, which mimics the real-world use of these estimators.

We generate SCMs by following a procedure similar to one for evaluating causal structure learning algorithms~\citep{heinze-stat18}.
SCMs and their corresponding datasets differ in the following characteristics: the number of variables \( d \); the probability of edge creation \( p_e \); the variance in the noise \( \sigma \); the strength of the noise term in non-source nodes \( \rho \).
In addition, we add parameters to control confounding \( \gamma \), mediation \( m \), and heterogeneity in treatment effects \( p_h, m_p \).
We initialize a set of parameters as follows:

\begin{enumerate}
    \item Number of variables: \( d \in [10, 20, 30] \)
    \item Edge probability parameter: \( p_e \in [0.1, 0.3, 0.5] \)
    \item Noise variance: \( \sigma \in [0.2, 0.4, 0.6] \)
    \item Magnitude in noise: \( \rho \in [0.1, 0.5, 0.9] \)
    \item Confounder: \( \gamma \in [\text{True, False}] \)
    \item Mediator chain length: \( m \in [0, 1, 2] \)
    \item Number of HTE inducing parents: \( p_h \in [0, 1, 2] \)
    \item HTE from mediating parents: \( m_p \in [\text{True, False}] \)
\end{enumerate}

Here, we describe how we generate synthetic datasets used in the experiments.
We first generate an adjacency matrix, \( \mathbf{A} \) over the \( d \) variables.
Assuming that the variables \( \{1, \dots, d \} \) are causally ordered.
For each pair of nodes \( i, j \), where \( i < j \) in the causal order, we sample from Bernoulli\(( p_e )\) to determine if there is an edge from \( i \) to \( j \).
With the adjacency matrix, we sample the corresponding coefficients, \( \mathbf{C} \), from Unif\(( -1, 1 )\).
Noise terms \( \epsilon_i \) are generated by sampling from a normal distribution with mean zero and variance \( \sigma \).
For a non-source node, $ i $ (nodes with parents), the noise term is multiplied by the magnitude in noise, $ \gamma $.

Next, we choose the treatment and outcome nodes based on two parameters, \( \gamma \) and \( m \). \( \gamma \) controls whether there is confounding in the generated dataset and \( m \) controls the length of the mediating chain (\( 0 \) means no mediator).
First, we select all pairs of nodes with a \emph{directed} \( m \) hop path.
Then, if \( \gamma \) is True, we filter all pairs that do not contain a \emph{backdoor path}.
Otherwise, we filter out pairs that contain a backdoor path.
Finally, we sample a pair of nodes from the filtered set to obtain the treatment and outcome nodes.
If no pairs satisfy the parameter setting (e.g., no confounding), the graph is resampled up to 100 times.
Otherwise, we discard the parameter setting.

A source node \( i \) (node with no parents) takes the value of the sampled noise terms, \( \epsilon_i \).
For a non-source node, \( i \), we consider several additive factors.
First is the base term from the coefficient matrix: \( C_{pa(i), i} X_{pa(i)} \), where \( C_{pa(i), i} \) are the coefficients from \( pa(i) \) to \( i \) and \( X_{pa(i)} \) are the values generated for \( pa(i) \).
Next are the interaction terms, which are determined by the mediating chain length \(m\), number of HTE inducing parents \( p_h \), and whether there is HTE from parents of mediators \( m_p \).
For mediators, we have: \( \mathbbm{1}_{m}(i) \mathbbm{1}_{m_p} X_{pa_m(i)} X_{pa_n(i)} \), where \( \mathbbm{1}_{m} \) is an indicator for if \( i \) is a mediator, \( \mathbbm{1}_{m_p} \) is an indicator for if there is HTE induced from mediator parents, and \( pa_m \) and \( pa_n \) represent mediating and non-mediating parents.
For outcomes, the interaction is defined similarly: \( \mathbbm{1}_{y}(i) X_{pa_m(i)} X_{pa_n(i)} \), where \( \mathbbm{1}_{y}(i) \) indicates whether the node is the outcome node.
The number of HTE parents determines how many non-mediator parents (\( pa_n \)) interact with the treatment.
A value of zero means there is no heterogeneity.
So the final data generating equation is:
\begin{align}
    X_i & = C_{pa(i), i} X_{pa(i)}
    \nonumber \\
    & + \mathbbm{1}_{m}(i) \mathbbm{1}_{m_p} X_{pa_m(i)} X_{pa_n(i)}
    \nonumber \\
    & + \mathbbm{1}_{y}(i) X_{pa_m(i)} X_{pa_n(i)} + \epsilon_i
\end{align}
For each parameter, we randomly sample from all other parameter values \( 500 \) times to generate SCMs.
For each of the $11,000$ SCMs, we create one dataset of size \( 10,000 \).

\subsubsection{IHDP dataset}
The covariates in the Infant Health and Development Program (IHDP) dataset come from a randomized experiment studying the effects of specialist home visits on future cognitive test scores~\citep{hill-jcgs11}.
Outcomes can be simulated using the covariates.
We simulate 1,000 instances of simulated outcomes and average over all results\footnote{Outcomes are generated using the JustCause package: https://github.com/inovex/justcause}.

\subsubsection{Real-world datasets}

We additionally explore three real-world datasets.
The first dataset is a League of Legends (LoL) dataset~\citep{he-fdg21}.
LoL is one of the most popular multiplayer online battle arena (MOBA) games.
Players compete in teams to capture the opposing team's base.
Riot Games, the developer of LoL, regularly updates the game's software through \emph{patches}.
Here, we use patches as the treatment variable and measure how it affects number of \emph{kills} on a team-level.

\emph{Cannabis} is a dataset consisting of tweets about the e-cigarette Juul and cannabis-related topics~\citep{adhikari-icwsm21}.
English tweets related to the e-cigarette Juul are collected from 2016 to 2018.
For the users who mention Juul, cannabis-related tweets are collected from 2014 to 2018.
In this dataset, stance towards Juul and cannabis in tweets are crowdsourced, and a model is trained for stance detection.
In this work, we are interested in how users' stance on Juul affects their stance on cannabis.
To do this, we find users who have Juul-related tweets occurring before a cannabis-related tweet and detect their stance on Juul and cannabis using the classifier from~\citep{adhikari-icwsm21}, where features are the word embeddings for tweets using BERT~\citep{devlin-nacl19}.
The treatment variable of interest is a user's stance on Juul, in favor (1) or against (0).
The outcome variable is their stance on cannabis, in favor (1) or against (0).
Since this dataset contains only one group, we bootstrap the estimators to compute standard deviation.

The third dataset is about a study of interventions that may improve performance among undergraduates at the University of Illinois at Urbana Champaign (UIUC): strengthening human adaptive reasoning and problem solving (SHARPS).
In this dataset, students were given a battery of cognitive, physical, and diet evaluations.
In subsequent weeks, they were given one of four interventions: no training (active control), fitness training, fitness and cognitive training, and fitness, cognitive, and mindfulness training.
Since most HTE estimators can only handle binary treatment variables, we measure the effect of fitness training compared to active training.
The outcome variable we are interested in is the change in the number of correct answers on a digit-symbol substitution task.
In order to compute standard deviation, we bootstrap estimators.

\subsection{Evaluation}

When the ground truth is known, a common metric is the \emph{precision in estimation of heterogeneous effects} (PEHE):
\begin{equation}
    PEHE = \sqrt{\frac{1}{n} \sum_i^n \Big( (Y_i(1) - Y_i(0)) - E[Y_i(1) - Y_i(0) \mid \mathbf{X}_i = \mathbf{x}_i] \Big)^2},
\end{equation}
where the expectation is the prediction of an HTE estimator.
This is the root mean squared error (RMSE) of the true effect and the predicted effect.

In practice, we may be more interested in how well an HTE estimator can predict whether an individual will benefit from treatment or not.
For synthetic datasets, we compute an \emph{accuracy of decisions}.
For an individual, $ i $, let the true individual effect be $ \tau_i $, and the predicted effect be $ \hat{\tau}_i $.
Let the indicator $ \mathbbm{1}\big[ sign(\tau_i) = sign(\hat{\tau}_i) \big] $ determine if the sign of the predicted effect matches the sign of the true effect.
The accuracy is defined as:
\begin{equation}
    Acc = \frac{1}{N} \sum_{i}^{N} \mathbbm{1}[sign(\tau_i) = sign(\hat{\tau}_i)].
\end{equation}

If the ground truth is unknown, several metrics for heuristically evaluating HTE estimators have been proposed~\citep{saito-cv19,schuler-model18}.
Therefore, we additionally report the $ \tau  $-risk for evaluation:
\begin{equation}
    \tau\text{-risk} = \frac{1}{n} = \sum_i^n \Big( (Y_i - \hat{m}(\mathbf{X}_i)) - (T_i - \hat{p}(\mathbf{X}_i)) \hat{\tau}(\mathbf{X_i}) \Big)^2,
\end{equation}

When comparing estimators across different parameter settings and SCMs, we also rank all methods based on their error for a particular SCM\@.
Then we average the rank of each method across SCMs and report on the average rank.
To study performance under specific settings, we also report on average ranks for subsets of parameter settings.

\input{05-01_synthetic.tex}

\input{05-02_semi.tex}

\input{05-03_real.tex}

\subsection{Synthetic dataset evaluation}\label{sec:synthetic_eval}

We first investigate how HTE estimation methods perform under different parameter settings (and thus SCMs).
Then we investigate the effect that causal feature selection has on HTE estimation.

\subsubsection{Multidimensional scaling}\label{sec:synthetic_mds}

\input{figures/mds_correlation.tex}

For each method and SCM, we compute the RMSE for the predicted effect.
This yields a (number of methods) $ \times $ (number of SCMs) matrix with RMSE values.
The correlation between pairs of SCM results is the similarity between two datasets.
Figure~\ref{fig:mds_correlation} shows a Multidimensional scaling (MDS) plot based on the correlation between SCM datasets.
The goal of the MDS plot is to visualize the similarity in performance of methods based on different SCM settings.
If two points are close together, then all learners in those experiments performed similarly.
The color and shape of a point represent the method that performed the best (lowest RMSE) for that SCM\@.

From this Figure, Causal Tree has the best performance in most settings, followed by Causal Forest and Balancing Neural Network.
Surprisingly, Causal Tree and Causal Forest are separated, meaning that the settings in which they perform the best have different RMSE results across all learners.
Balancing Neural Network occupies two regions close to Causal Tree and Causal Forest, with a majority near Causal Forest and TARNet.
This makes sense since TARNet is an extension of Balancing Neural Network.

\input{figures/pairwise_table.tex}

\subsubsection{Pairwise comparison}

Next, we investigate if there are methods that beat others consistently.
We compute the percentage of wins in Table~\ref{tab:pairwise} for which there is a 5\% decrease in RMSE\@.
For example, SLearner beats XLearner by a 5\% margin in only 0.34\% of scenarios, while XLearner beats SLearner by a 5\% margin in 50.12\% of scenarios.
From these results, there does not appear to be an overly dominant method.
For example, Causal Forest generally beats DRLearner by the margin but does not generally beat by Causal Tree or Balancing Neural Network.

From this Table, there is no clear dominant method, but there are methods which are dominated by others.
For example, SLearner generally does not perform well, especially its other meta-learning counterparts.
Causal Tree has an overall high percentage of wins but is also beaten by other methods at a relatively high percentage.
This suggests methods need to be carefully chosen for different settings and assumptions.

\subsubsection{Performance in different SCMs}

\input{figures/winner_table.tex}

In order to investigate how HTE estimators perform in different types of SCMs, we show the average rank under different parameter settings in Table~\ref{tab:mse_rankings} and the average RMSE in Table~\ref{tab:rmse_table}.
The Oracle method (Known Structure) is at the bottom of each table.
Bolded results indicate the best rank in the column (ignoring the Oracle).
In terms of ranking, structure learning (FSL) achieves the lowest average rank over most datasets.
Additionally, since some ranks are close, there is no clear winner among closely ranked methods, which supports the pairwise comparison between HTE estimation methods in Table~\ref{tab:pairwise}.

\input{figures/mse_table.tex}

We also show the average RMSE for each method in the same parameter settings in Table~\ref{tab:rmse_table}.
Here, Causal Tree and Balancing Neural Network achieve lower average RMSE.
DRLearner has high error and variance in predictions, which suggests it requires more fine-tuning of its estimators.

Looking at specific parameters, we see that SCMs with mediators generally result in higher errors than without mediators.
This makes sense since including mediators will bias the effect estimation.
Interestingly, a chain length of 1 is more challenging than a chain length of 2.
Increasing the number of HTE inducing parents also results in higher errors since 0 HTE parents imply average effect only, and more HTE inducing parents increases the complexity of the heterogeneous effect.
Causal Tree seems to be the least affected by parameter changes compared to other HTE estimators.

\subsection{Accuracy of decisions}

In practice, practitioners may be interested in how well an HTE estimator can assign a treatment decision to individuals, rather than the magnitude in the change in outcomes.
To investigate this, we explore how well an HTE estimator can decide whether it is beneficial to treat an individual or not.
To do this, we binarize the ground truth treatment effect into positive and negative effects.
For each HTE estimator, we binarize the predictions into positive and negative predictions.
Then, an accuracy is computed over multiple datasets, signifying how many examples are labeled correctly.

Table~\ref{tab:accuracy_table} shows the accuracy of decisions for each HTE estimator, using the same columns as in Tables~\ref{tab:mse_rankings} and~\ref{tab:rmse_table}.
In general, we see that Causal Tree achieves the highest accuracy of decisions, followed by Balancing Neural Network and FSL.
Meta-learners achieve very low accuracy when there are more invalid variables, such as when there are mediators.

\input{figures/accuracy_table.tex}

\subsection{Comparing evaluation metrics}

Since ground truth causal effects are generally not available in real-world datasets, several heuristic HTE evaluation metrics have been developed.
We consider inverse propensity weighting (IPW)~\citep{schuler-model18}, $\tau$-risk~\citep{schuler-model18,nie-rlearner17}, plugin validation (plugin $\tau$)~\citep{saito-cv19}, and counterfactual cross-validation (CFCV)~\citep{saito-cv19}.
Each metric uses estimators for propensity scores or outcomes to compute an error.
Each metric uses the full set of features in our experiments to compute a propensity score or outcome.
We use Logistic Regression for each metric that requires an estimate of the propensity score.
For $\tau$-risk, we use Linear Regression for outcome prediction.
For plugin validation and CFCV, we use Counterfactual Regression (CFR)~\citep{shalit-icml17} for the regression function.
In order to evaluate the reliability of these heuristic metrics, we compute the Spearman rank correlation of each metric compared to the ground truth ranking.
For each dataset, we rank each learner based on the heuristic error on the test set.
Then we compute a Spearman correlation in one of two ways.
The first way computes a correlation on the entire ranking for each dataset, and an average can be computed over all datasets.
The second way computes a correlation of the rank for each learner over all the datasets.

Table~\ref{tab:metric_dataset} shows the correlation of each heuristic metric ranking to the ground truth ranking over the synthetic datasets and the semi-synthetic (IHDP) datasets.
This table shows that most metrics are positively correlated with the ranking from the ground truth on synthetic datasets, except for plugin validation.
However, there is a large variance, and the average correlation is generally low, which suggests that having a hidden structural causal model with invalid variables may affect the metrics.
For the IHDP dataset, we see that all metrics except plugin validation positively correlate with the true ranking, with $ \tau  $-risk having the highest correlation.
Since the IHDP dataset does not contain mediators or descendants of treatment, this may improve the correlation of $ \tau  $-risk.

\begin{table*}[t]
    \centering

    \resizebox*{\linewidth}{!}{
        \begin{tabular}{l *{8}{M{1.7cm}}}
            & \multicolumn{7}{c}{Synthetic} \\
            \cmidrule(lr){2-8}

            & Overall & Has Mediator & No Mediator & Contains Confounding & No Confounding & 0 HTE Parents & 1 or 2 HTE Parents &  IHDP \\
            \midrule
            IPW &  0.10 $ \pm $ 0.55 & 0.03 $\pm$ 0.59 & 0.01 $\pm$ 0.43 & 0.02 $\pm$ 0.44 & 0.08 $\pm$ 0.45 & 0.11 $\pm$ 0.42 & 0.03 $\pm$ 0.44 & 0.10 $ \pm $ 0.43 \\
            $ \tau $-risk & 0.20 $\pm$ 0.47 & 0.01 $\pm$ 0.45 & 0.10 $\pm$ 0.43 & 0.13 $\pm$ 0.44 & 0.07 $\pm$ 0.43 & 0.13 $\pm$ 0.47 & 0.09 $\pm$ 0.43 & 0.53 $ \pm $ 0.32\\
            Plugin $\tau$ & -0.18 $\pm$ 0.55 & -0.23 $\pm$ 0.56 & -0.27 $\pm$ 0.44 & -0.27 $\pm$ 0.43 & -0.25 $\pm$ 0.44 & -0.24 $\pm$ 0.45 & -0.26 $\pm$ 0.43 & -0.34 $ \pm $ 0.38\\
            CFCV & 0.25 $\pm$ 0.59 & 0.14 $\pm$ 0.61 & 0.06 $\pm$ 0.47 & 0.03 $\pm$ 0.47 & 0.03 $\pm$ 0.45 & 0.15 $\pm$ 0.56 & 0.01 $\pm$ 0.47 & 0.20 $ \pm $ 0.44 \\
        \end{tabular}
    }
    \caption{Comparing the rank correlation of metrics with true performance.}\label{tab:metric_dataset}
    \vspace{1em}
\end{table*}
\section{Heterogeneous Treatment Effects}\label{sec:background}

In this section, we provide background on heterogeneous treatment effect estimation. We define the problem of heterogeneous treatment effect estimation and common assumptions to allow for estimation.
We introduce structural causal models~\citep{pearl-book09} and how they can be used to assess heterogeneity.

\subsection{Data model}\label{sec:data_model}

Before formalizing the problem and data model, we introduce the notation used in this survey.
We denote a space of random variables with a calligraphic letter (e.g., $ \mathcal{X} $).
We denote random variables with a capital letter (e.g., $Y$) and random vectors with a boldfaced capital letter (e.g., $\mathbf{X}$).
Subscripts signify the random variable associated with a specific unit (e.g., $Y_i$ or $\mathbf{X}_i$).
Lowercase letters denote an assignment of a random variable or vector (e.g., $Y_i = y$ or $\mathbf{X}_i = \mathbf{x}$).
Each unit $i$ has a vector of ``pre-treatment'' features (or covariates), $\mathbb{X}_i \in \mathcal{X}$, a binary treatment indicator, $T_i$, where $T_i=0$ and $T_i=1$ indicates $i$ is not treated and treated, respectively, and two outcomes when not treated and treated, $Y_i(T_i=0)$ and $Y_i(T_i=1)$, respectively.
For brevity, we write $Y_i(T_i=0)$ and $Y_i(T_i=1)$ as $Y_i(0)$ and $Y_i(1)$, respectively.
In this survey, we consider the treatment to be binary, as most work in estimating heterogeneous treatment effects only considers binary treatment, some recent work has considered how to estimate heterogeneous effects when the treatment is non-binary~\citep{tran-aaai19,schwab-arxiv19,kennedy-jrss17}
The underlying data distribution is represented by a tuple $(\mathbf{X}_i, T_i, Y_i(0), Y_i(1))^N$ of $N$ independent and identically distributed individuals.
The fundamental problem of causal inference is that for each individual, we can observe either the potential outcomes under treatment ($Y_i(1)$) or under control ($Y_i(0)$), but not both.
Therefore, the observed dataset as $\mathcal{D} = (\mathbf{x}_i, t_i, y_i)^N$ with $N$ individuals where $Y_i$ is defined as:
\begin{equation}
    \label{eq:potential_outcomes}
    y_i = t_i Y_i(1) + (1 - t_i) Y_i(0).
\end{equation}
We denote $ Y_i $ as the factual outcome and the counterfactual outcome, $ Y_i' $, as the outcome that would have occurred if the treatment was different.
We define the treatment group as all individual units $i$ such that $t_i = 1$ and the control group as all units $i$ such that $t_i = 0$.

\subsection{Heterogeneous Treatment Effect Estimation}\label{sec:hte}

The average treatment effect (ATE) over the population is defined as the expected difference in potential outcomes:
\begin{equation}
    \label{eq:ate}
    \text{ATE} := E[Y(1) - Y(0)].
\end{equation}
Since each unit reacts differently to treatment, we define the individual treatment effect (ITE) as the difference in the unit's potential outcomes, denoted by $\tau_i$:
\begin{equation}
    \label{eq:ite}
    \tau_i := Y_i(1) - Y_i(0).
\end{equation}
However, due to the fundamental problem of causal inference, $\tau_i$ can not be observed for any individual unit.
Instead, HTE estimation focuses on estimating the conditional average treatment effect (CATE):
\begin{equation}\label{eq:cate}
    \tau(\mathbf{x}) := \mathbb{E} \left[ Y_i(1) - Y_i(0) \mid \mathbf{X}_i = \mathbf{x} \right].
\end{equation}
The CATE describes the expected treatment effect for units in a subpopulation with features $\mathbf{x}$.
An HTE estimator, defined as $\hat{{\tau}}(\mathbf{X})$, seeks to estimate the CATE defined in equation~\eqref{eq:cate}.
The best estimator for CATE is also the best estimator for ITE in terms of the mean squared error~\citep{kunzel-pnas19}.
The goal of HTE estimation is to develop estimators of CATE which have the smallest root mean squared error for the ITE\@.
\begin{problem}\label{prob:hte_problem}
(Heterogeneous treatment effect estimation)
Given a dataset $\mathcal{D}=(\mathbf{X}_i, Y_i, T_i)$ of N instances with some true ITE, $\vect{\tau} = (\tau_i)^N$, estimate CATE $\hat{\tau}(\mathbf{X})$ such that the root-mean-squared error (RMSE) is minimized:
\begin{equation}\label{eq:cate_objective}
    \minl_{\hat{\tau}} \text{RMSE} = \frac{1}{N} \sum_{i=1}^N  (\tau_i - \hat{\tau}(\mathbf{X}_i)) ^ 2.
\end{equation}
\end{problem}
In practice, equation~\eqref{eq:cate_objective} cannot be estimated in real-world data.
HTE estimators are often evaluated on synthetic or semi-synthetic datasets~\citep{curth-neurips21}, but this may not capture performance in real-world scenarios.
Some work has been done on evaluating HTE estimators heuristically, \citep{schuler-model18,saito-cv19,schwab-perfect18}, which we describe in Section~\ref{sec:data_driven_background}.

\subsection{HTE Assumptions}\label{sec:hte_assumptions}

To estimate HTEs from observational data, some assumptions on the set of features and data given are needed.
We present some common assumptions used in heterogeneous treatment effect estimators.

\begin{assumption}\label{assumption:strong_ignorability}
    (Strong Ignorability) Given a binary treatment $T$, potential outcomes $Y(0)$, $Y(1)$, and a vector of features $\mathbf{X}$, strong ignorability states:
    \begin{align}
        \{ Y(0), Y(1) \} \perp T \mid \mathbf{X},
    \end{align}
\end{assumption}
This states the \emph{potential} outcomes are conditionally independent of the treatment conditioned on the features.
When strong ignorability holds, $\mathbf{X}$ is said to be ``admissible'', meaning treatment effects can be estimated without bias.
This assumption is also known as unconfoundedness or the assumption of no hidden confounding.

\begin{assumption}\label{assumption:overlap}
    (Overlap) For every unit $i$, $0 < P(T_i = 1) < 1$ with probability 1.
\end{assumption}
This assumption states that all units have a non-zero probability of assignment to each treatment condition.
Assumption~\ref{assumption:overlap} is also known as positivity or common support.

Another common assumption used is the Stable Unit Treatment Value Assumption (SUTVA) or the assumption of no interference~\citep{rubin-jasa80}.
\begin{assumption}\label{assumption:sutva}
    (Stable Unit Treatment Value Assumption (SUTVA)) Given any two units $i$ and $j$, the potential outcomes of $i$, $\{Y_i(0), Y_i(1) \}$ are unaffected by the treatment assignment of unit $j$:
    \begin{equation}
        \{Y_i(0), Y_i(1)\} \perp T_j, \forall j \neq i.
    \end{equation}

\end{assumption}
With the SUTVA assumption, individual outcomes are not affected by the treatment or outcomes of peers.
Though there has been some work focused on estimating HTEs without the SUTVA assumption~\citep{yuan-www21,tran-aaai22}, a majority of work in HTE estimation has focused on non-networked settings, which we focus on in this survey.

\subsection{Counterfactual prediction}\label{sec:counterfactuals}

The first category of methods predicts counterfactuals from data, the outcome when treated and control, $ \hat{Y}(1), \hat{Y}(0) $, respectively, and uses these predictions for effect estimation, $ \hat{Y}(1)-\hat{Y}(0)$.

\subsubsection{Single-model approach}

The simplest way to make counterfactual predictions is by using the \emph{single-model approach}.
The single-model approach estimates CATE using a single regression or classification model.
Given the data \( D = (\mathbf{X}_i, Y_i, T_i) \), a supervised method is trained to predict \( Y_i \) given \( \mathbf{X}_i, T_i \), \( \hat{\mu}(T_i, \mathbf{X}_i) = \hat{Y}_i \).
The CATE can be predicted as the difference between the prediction with \( T_i = 1 \) and \( T_i = 0 \).
\begin{equation}
    \hat{\tau}(\mathbf{x}) = \hat{\mu}( T = 1, \mathbf{X} = x) - \hat{\mu}(T = 0, \mathbf{X} = x).
\end{equation}
Hill used Bayesian Additive Regression Trees (BART) to predict treatment effects~\citep{hill-jcgs11}.
Athey and Imbens used regression trees as the single-model approach~\citep{athey-stat15}.
K{\"{u}}nzel et al.\ used a single-model approach with gradient-boosted trees and random forests~\citep{kunzel-pnas19} and name these single-model approaches S-Learner.

Since the single-model approach only requires regression or classification, any off-the-shelf supervised learning algorithm can be used.
However, the single-model approach may not be able to capture potential differences in treated and control groups and bias estimation.
Additionally, in some cases, it may not be possible for the treatment to be used in prediction (e.g., not used as a split in a regression tree).
This results in the estimator predicting the same outcome for different levels of treatment, which results in a zero for CATE estimation.

Besides methods using classical machine learning methods, several deep learning approaches for estimating HTEs have been proposed recently using the single-model approach.
Atan et al.\ propose a deep learning model that consists of two stages~\citep{atan-aaai18}.
In the first stage, Deep-Treat uses a bias-removing auto-encoder to learn a balanced representation, $ \Phi(\mathbf{x}) $, of treatment and control groups.
The balancing is learned through cross-entropy loss between the marginal and conditional treatment distributions given the hidden representation.
Deep-Treat trains a neural network to predict the outcome using the learned representation and the treatment in the second state.

\subsubsection{Two-model approaches}

In contrast to the single-model approach, another way to model heterogeneity is to model two potential outcomes separately, called the \emph{two-model approach}.
In the two-model approach, a separate model is built on treated and control groups.
Then, to estimate CATE, the difference between treatment group model prediction and control group model prediction is used to estimate the effect.

Formally, let \( \hat{\tau}_1(\mathbf{x}) \) and \( \hat{\tau}_2(\mathbf{x}) \) be supervised estimators for the outcome for the treated and control groups, respectively.
The CATE is estimated as:
\begin{align}
    \hat{\tau}(\mathbf{x}) & = E[Y \mid T = 1, \mathbf{X} = \mathbf{x}] - E[Y \mid T = 0, \mathbf{X} = \mathbf{x}] \\
    & = \hat{\tau}_1(\mathbf{x}) - \hat{\tau}_2(\mathbf{x}).
\end{align}
The main difference between the two-model and single-model approaches is that the two-model approach trains two separate estimators for treatment and control outcome prediction instead of one model for both.
The two model approach has been used with Linear Regression~\citep{cai-bio11}, Trees~\citep{athey-stat15,hill-jcgs11}, Random Forests~\citep{foster-statmed11}, and Gradient Boosted Trees~\citep{kunzel-pnas19}.
Figure~\refeq{fig:slearner_tlearner} contrasts how single-model and two-model approaches perform training and validation.

\begin{figure}[h]
    \centering
    \begin{subfigure}[b]{0.45\linewidth}
        \centering
        \includegraphics[width=\textwidth]{figures/ppt/SLearner.png}
        \caption{Single-model approach.}\label{fig:slearner}
    \end{subfigure}
    \hfill
    \begin{subfigure}[b]{0.45\linewidth}
        \centering
        \includegraphics[width=\textwidth]{figures/ppt/TLearner.png}
        \caption{Two-model approach.}\label{fig:tlearner}
    \end{subfigure}
    \caption{A comparison of training and prediction of effects using single-model and two-model approaches.
    The single-model approach uses one estimator for training, and for prediction, features are appended with a treatment indicator of 1 and a treatment indicator of 0 for two separate predictions based on the single model. Then the effect is estimated as the difference.
    In the two-model approach, training is done using two separate models for treated and control groups.
    For prediction, the features are input into the two separate models, and a difference is computed as the predicted effect.}\label{fig:slearner_tlearner}
\end{figure}

In addition, recent work has applied the two-model approach with deep learning estimators.
Johansson et al.\ proposed Balancing Neural Network (BNN), which aims to leverage representation learning and balanced representations for treatment and control groups~\citep{johansson-icml16}.
They learn a neural network architecture with two connecting parts to predict the outcome.
The first part learns the representation of the features, $ \Phi(\mathbf{x}) $.
In the second part, they append the treatment to the hidden representation and learn a representation to predict the outcome.
Separately, they constrain $ \Phi(\mathbf{x}) $ to be similar between treated, $ \Phi_{t=1}(\mathbf{x}) $, and control, $ \Phi_{t=0}(\mathbf{x})$.
The intuition is that the hidden representations should be ``balanced''.
The basic BNN architecture is shown in Figure~\ref{fig:bnn}.

Shalit et al.\ extend the idea of BNNs to the Treatment-Agnostic Representation Network (TARNet) and Counterfactual Regression (CFR)~\citep{shalit-icml17}.
They extend the discrepancy metric by introducing a more flexible family of algorithms in Integral Probability Metrics (IPM) for balancing distributions.
Instead of considering a one-headed network that predicts an outcome based on appending the treatment to the hidden representation, they consider two heads predicting outcome when treated and not treated separately.
The intuition is that the treatment indicator appended to $ \Phi(\mathbf{x}) $ may lose importance in the second part of the network, and so explicitly separating the final predictions can capture differences in treatment and control.

Liu et al.\ propose Split Covariate Representation Network (SCRNet)~\citep{liu-acml20}.
In this work, they propose to divide the covariate space into four categories: confounders, instrumental variables, adjustment variables, and irrelevant variables.
Their goal is to identify the category of each feature and use only relevant features for training a neural network.
The intuition is that irrelevant features are unnecessary for estimation and instrumental variables may bias effect estimations~\citep{pearl-arxiv12,myers-aje11,middleton-pa16}.
Instead, they identify confounders, which confound or are parents of the treatment in outcome in an SCM, and parents of outcome only which are referred to as adjustment variables.
Then, they train a neural network where confounders and adjustment variables are input separately and have their own hidden representations before being concatenated and separated into a two-head network, similar TARNet.
The main difference is the splitting of the initial input into confounding and adjustment variables.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.7\textwidth]{figures/ppt/bnn.png}
    \caption{A high-level architecture of the Balancing Neural Network (BNN) from~\citep{johansson-icml16}. The first part learns a representation of the features only, $ \Phi $.
    In the second part, the treatment is appended, and the neural network is optimized for prediction performance and minimizing the representations of treated and control representations.}\label{fig:bnn}
\end{figure}

\subsubsection{Multi-Model Approaches}

Related to the TARNet approach~\citep{shalit-icml17}, Schwab et al.\ propose Perfect Match~\citep{schwab-arxiv19}.
They extend TARNet in two ways.
First, they propose using $ k $ heads to predict multiple treatment outcomes.
Second, they propose a new way to optimize the neural network model based on nearest neighbor matching and the precision in estimation of heterogeneous effect (PEHE)~\citep{hill-jcgs11}, which measures the mean-squared error in HTE estimation, called NN-PEHE\@.
They propose computing an imputed treatment effect based on nearest neighbors.
Since the true causal effect is unknown for individual units, they substitute the true counterfactual outcome with the outcome from the nearest neighbor.
Their goal is to minimize NN-PEHE for training Perfect Match.
Shi et al.\ proposed another extension of TARNet, called Dragonnet~\citep{shi-neurips19}.
Instead of a two-head neural network, they introduce an additional head which predicts the propensity score from the hidden layer.
In addition, they introduce a modification of the objective function for training neural networks called targeted regularization.

\section{Heterogeneous Treatment Effects}\label{sec:background}

In this section, we provide background on heterogeneous treatment effect estimation, starting with notations from the perspective of the potential outcome framework~\citep{rubin1974}. Then we define the problem of heterogeneous treatment effect estimation and common assumptions to allow for estimation.
Finally, we introduce structural causal models~\citep{pearl-book09} and how they can be used to assess heterogeneity.

\subsection{Data model}\label{sec:data_model}

Before formalizing the problem and data model, we introduce the notations used in this survey.
We denote a space of random variables with a calligraphic letter (e.g., $ \mathcal{X} $).
We denote random variables with a capital letter (e.g., $Y$) and random vectors with a boldfaced capital letter (e.g., $\mathbf{X}$).
Subscripts signify the random variable for a specific unit (e.g., $Y_i$ or $\mathbf{X}_i$).
Lowercase letters denote an assignment of a random variable or vector (e.g., $Y_i = y$ or $\mathbf{X}_i = \mathbf{x}$).
Each unit $i$ has a vector of ``pre-treatment'' features or covariates, $\mathbb{X}_i \in \mathcal{X}$, a binary treatment indicator, $T_i$, where $T_i=0$ and $T_i=1$ indicates $i$ is not treated and treated, respectively, and two outcomes when not treated and treated, $Y_i(T_i=0)$ and $Y_i(T_i=1)$, respectively.
The underlying data distribution is represented by a tuple $(\mathbf{X}_i, T_i, Y_i(0), Y_i(1))^N$ of $N$ independent and identically distributed individuals.
The average treatment effect (ATE) over the population is defined as the expected difference in potential outcomes:
\begin{equation}
    \label{eq:ate}
    \text{ATE} := E[Y(1) - Y(0)].
\end{equation}
The fundamental problem of causal inference is that for each individual, we can observe either the potential outcomes under treatment ($Y_i(1)$) or control ($Y_i(0)$), but not both.
Therefore, we define the observed dataset as $\mathcal{D} = (\mathbf{x}_i, t_i, y_i)^N$ with $N$ individuals where $Y_i$ is defined as:
\begin{equation}
    \label{eq:potential_outcomes}
    Y_i = T_i Y_i(1) + (1 - T_i) Y_i(0).
\end{equation}
We denote $ Y_i $ as the factual outcome and the counterfactual outcome, $ Y_i' $, as the outcome that would have occurred if the treatment was different.

\subsection{Heterogeneous Treatment Effect Estimation}\label{sec:hte}

Since each unit reacts differently to treatment, we define the individual treatment effect (ITE) as the difference in the unit's potential outcomes, denoted by $\tau_i$:
\begin{equation}
    \label{eq:ite}
    \tau_i := Y_i(1) - Y_i(0).
\end{equation}
However, due to the fundamental problem of causal inference, $\tau_i$ is not observed for any individual.
The quantity of interest is often the conditional average treatment effect (CATE):
\begin{equation}
    \tau(\mathbf{x}) := \mathbb{E} \left[ Y_i(1) - Y_i(0) \mid \mathbf{X}_i = \mathbf{x} \right].
\end{equation}
The CATE describes the expected treatment effect for units in a subpopulation with features $\mathbf{x}$.
The best estimator for CATE is also the best estimator for ITE in terms of the mean squared error~\citep{kunzel-pnas19}, though in general, the ITE is not the same as CATE~\citep{vegetabile-arxiv21}.
The goal of HTE estimation is to develop estimators of CATE which have the smallest root mean squared error for the ITE.
\begin{problem}\label{prob:hte_problem}
    (Heterogeneous treatment effect estimation)
    Given a dataset $\mathcal{D}=(\mathbf{X}_i, Y_i, T_i)$ of N instances with some true ITE. $\tau_i$, estimate CATE $\hat{\tau}(\mathbf{X})$ such that the root mean squared error (RMSE) is minimized:
    \begin{equation}\label{eq:cate_objective}
         \minl_{\hat{\tau}} \text{RMSE} = \frac{1}{N} \sum_{i=1}^N  (\tau_i - \hat{\tau}(\mathbf{X}_i)) ^ 2.
    \end{equation}
\end{problem}

\subsection{HTE Assumptions}\label{sec:hte_assumptions}

In order to estimate HTEs from observational data, some assumptions on the set of features and data given are needed.
We present some common assumptions used in heterogeneous treatment effect estimators.

\begin{assumption}\label{assumption:strong_ignorability}
    (Strong Ignorability) Given a binary treatment $T$, potential outcomes $Y(0)$, $Y(1)$, and a vector of features or covariates $\mathbf{X}$, strong ignorability states:
    \begin{align}
        \{ Y(0), Y(1) \} \perp T \mid \mathbf{X},
    \end{align}
\end{assumption}
This states the \emph{potential} outcomes are conditionally independent of the treatment conditioned on the features.
When strong ignorability holds, $\mathbf{X}$ is said to be ``admissible'', meaning treatment effects can be estimated without bias.
This assumption is also known as unconfoundedness or the assumption of no hidden confounding.

\begin{assumption}\label{assumption:overlap}
    (Population Overlap) $0 < P(T_i = 1) < 1$ with probability 1.
\end{assumption}
This assumption states that all units have a non-zero probability of assignment to each treatment condition.
Assumption~\ref{assumption:overlap} is also known as positivity or common support.

Another common assumption used is the Stable Unit Treatment Value Assumption (SUTVA) or the assumption of no interference~\citep{rubin-jasa80}.
\begin{assumption}\label{assumption:sutva}
    (Stable Unit Treatment Value Assumption (SUTVA)) Given two individuals $i, j$, the potential outcomes of $i$, $\{Y_i(0), Y_i(1) \}$ are unaffected by the treatment assignment of other individuals $\forall j$:
    \begin{equation}
        \{Y_i(0), Y_i(1)\} \perp T_j, \forall j \neq i.
    \end{equation}
    SUTVA is also often referred to as the assumption of ``no interference''.
\end{assumption}
With the SUTVA assumption, individual outcomes are not affected by treatment or outcomes of peers.
Though there has been some work focused on estimating HTEs without the SUTVA assumption~\citep{yuan-www21,tran-aaai22}, a majority of work in HTE estimation has focused on non-networked settings, which we focus on in this survey.

\subsection{Data-driven HTE Estimation}\label{sec:data_driven_background}

Given an observed dataset $\mathcal{D}$, the goal of an \emph{HTE estimator} is to estimate CATE, $\hat{\tau}$, using the available features $\mathbf{X}$ in the data.
HTE estimators estimate CATE in a variety of ways.
Some HTE estimators rely on recursive partitioning using tree-based methods~\citep{tran-aaai19,athey-pnas16,athey-annals19}.
Other HTE estimators use predictions from machine learning methods to train a causal effect estimator~\citep{kunzel-pnas19,kennedy-arxiv20} or to predict counterfactual outcomes for estimation of effect~\citep{shalit-icml17,johansson-icml16}.
We describe these methods in Section~\ref{sec:survey}.

Since the true ITE is unknown in real-world scenarios, several heuristic HTE evaluation metrics have been developed for evaluating the performance of HTE estimators.
These metrics include inverse propensity weighting (IPW) validation~\citep{schuler-model18}, $ \tau $-risk~\citep{schuler-model18,nie-rlearner17}, nearest-neighbor imputation (NN-PEHE)~\citep{schwab-perfect18}, plugin validation~\citep{saito-cv19}, and counterfactual cross-validation (CFCV)~\citep{saito-cv19}.

IPW validation uses inverse propensity score weighting to a plugin estimate for the CATE, $ \tilde{\tau} $~\citep{schuler-model18}:
\begin{equation}
    \tilde{\tau}_{IPW} (\mathbf{X}_i, Y_i, T_i) = \frac{T_i}{\hat{e}(\mathbf{X}_i)}Y_i - \frac{1 - T_i}{1 - \hat{e}(\mathbf{X}_i)} Y_i.
\end{equation}
Here, $ \hat{e}(\mathbf{x}) $ is an estimate of the propensity score (i.e., $\mathbb{E}[T_i|\mathbf{X}_i]$)
This plugin is used to approximate the mean-squared error for the estimator $ \hat{\tau} $:
\begin{equation}
    \hat{R}_{IPW}(\hat{\tau}) = \frac{1}{n} \sum_{i=1}^n \Big(\tilde{\tau}_{IPW} (\mathbf{X}_i, Y_i, T_i) - \hat{\tau}(\mathbf{X}_i) \Big)^2.
\end{equation}

The $ \tau  $-risk takes as input an estimator $ \hat{\tau}  $ and approximates the objective defined in~\eqref{eq:cate_objective}.
To do this, two functions are defined: $ \hat{m}(\mathbf{X}_i) $ which is an outcome estimator for unit $i$ (i.e., $ \mathbb{E}[Y_i | \mathbf{X}_i] $) and $ \hat{e}(\mathbf{X}_i) $ which is an estimate of the propensity score.
Both $ \hat{m} $ and $ \hat{e} $ are learned through supervised techniques and optimized for predictive accuracy~\citep{nie-rlearner17}.
The $ \tau  $-risk is defined as:
\begin{equation}
    \tau\text{-risk} = \frac{1}{N} \sum_{i=1}^N \Big( \big(Y_i-\hat{m}(\mathbf{X}_i)\big) - \big(T_i - \hat{e}(\mathbf{X}_i)\big)\hat{\tau}(\mathbf{X}_i) \Big)^2.
\end{equation}

Nearest neighbor precision in estimation of heterogeneous effects (NN-PEHE) uses matching to approximate a treatment effect by substituting counterfactual outcomes from matched neighbors~\citep{schwab-perfect18}.
A unit $i$ is matched with one\footnote{Matching can be done with multiple units, and an aggregate can be computed.} control and one treated unit, denoted by NN$_0 (i)$ and NN$_1 (i)$.
Using the matched units, a substitute of the ground-truth treatment effect is computed as $ y(\text{NN}_1(i)) - y(\text{NN}_0(i)) $.
The NN-PEHE error, $ \hat{\epsilon}_{\text{NN-PEHE}} $ is computed as:
\begin{equation}
    \hat{\epsilon}_{\text{NN-PEHE}} = \frac{1}{N} \sum_{i=1}^n \bigg( \big(y(\text{NN}_1(i)) - y(\text{NN}_0(i)) \big) - \hat{\tau}(\mathbf{X}_i) \bigg)^2.
\end{equation}

Plugin validation uses arbitrary machine learning algorithms to estimate a plugin $ \tau $, similar to IPW validation.
To get a plugin estimate, two predictions are needed, $ f_1(\mathbf{x}) $ and $ f_0(\mathbf{{x}}) $, which are predictions of the outcome when treated and not treated, respectively, using only the features.
Then, the plugin validation metric is defined as:
\begin{equation}
    \hat{R}_{plugin}(\hat{\tau}) = \frac{1}{N} \sum_{i=1}^N \bigg( (f_1(\mathbf{X}_i) - f_0(\mathbf{X}_i)) - \hat{\tau}(\mathbf{X}_i) \bigg)^2.
\end{equation}

Counterfactual cross-validation (CFCV) uses the plugin estimate of the CATE in a doubly robust formulation~\citep{saito-cv19}:
\begin{equation}
    \tilde{\tau}_{DR}(\mathbf{X}_i, Y_i, T_i) = \frac{T - \hat{e}(\mathbf{X}_i)}{\hat{e}(\mathbf{X}_i) (1 - \hat{e}(\mathbf{X}_i))} (Y_i - f_T(\mathbf{X}_i)) + f_1(\mathbf{X}_i) - f_0(\mathbf{X}_i),
\end{equation}
where $ \hat{e} $ is an estimate of the propensity score and $ f_t(\mathbf{x}) $ is an arbitrary regression function for a specific treatment value $ t $ that predicts the outcome given the set of features.
Then, the plugin estimate $ \tilde{\tau}_{DR} $ can be used in the mean-squared error estimate:
\begin{equation}
    \hat{R}_{CFCV} = \frac{1}{N} \sum_{i=1}^N \bigg( \tilde{\tau}_{DR}(\mathbf{X}_i, Y_i, T_i) - \hat{\tau}(\mathbf{X}_i) \bigg)^2
\end{equation}

Each metric is an estimate of the error of an HTE estimator and performance based on these heuristic metrics has been shown to correlate with true performance.
Since these metrics are heuristic, we explore their reliability in synthetic and semi-synthetic experiments in Section~\ref{sec:experiments}.

\subsection{Structural Causal Models}\label{sec:scm_background}

\begin{figure}
    \centering
    \newcommand{0.25}{0.5}
        \begin{subfigure}[b]{\columnwidth}
            \centering
            \begin{tikzpicture}
                \node (t) at (-0.75, 0) [label=below:$T$,point];
                \node (x) at (0.0, 0.75) [label=left:\(X\),point];
                \node(y) at (0.75, 0.0) [label=below:$Y$,point];
                \node (ux) at (0.0, 1.5) [label=above:\(U_X\),point];
                \node (ut) at (-1.5, 0.75) [label=above:\(U_T\),point];
                \node (uy) at (1.5, 0.75) [label=above:\(U_Y\),point];

                \path (x) edge (t);
                \path (x) edge (y);
                \path (t) edge (y);
                \path (ux) edge (x);
                \path (uy) edge (y);
                \path (ut) edge (t);
            \end{tikzpicture}
        \end{subfigure}
    \caption{An example causal model.
    }\label{fig:example_scm}
\end{figure}

Now we introduce an alternative to data-driven estimation of causal effects through structural causal models (SCMs) and how they relate to heterogeneous effects.

A structural causal model (SCM), $M$, consists of two sets of variables, $\mathbf{U}$ and $\mathbf{V}$, and a set of structural equations, $\mathbf{F}$, describing how values are assigned to each \emph{endogenous} variable $V^{(i)} \in \mathbf{V}$ based on the values of $\mathbf{v}$ and $\mathbf{u}$~\citep{pearl-book09}: \( v^{(i)} = f^{(i)}(\mathbf{v}, \mathbf{u}) \).
A causal graph \( G \) is a directed acyclic graph that captures the causal relationships among the variables.
Each variable, $V^{(i)}$, is represented by a node in the causal graph, an edge from $V^{(i)}$ to $V^{(j)}$ signifies that $V^{(i)}$ is a cause of $V^{(j)}$.
We also refer to the edge relationship as a parent-child relationship, where the parent is a cause of the child.
The variables in $\mathbf{U}$ are considered \emph{exogenous}, and the causal graph may omit representing $\mathbf{U}$ explicitly.
Figure~\ref{fig:example_scm} shows an example of a causal graph with three exogenous variables, $ U_X, U_T, U_Y $ and three endogenous variables, $ X, T, Y $.
An example set of associated structural equations is shown in eq~\eqref{eq:see1}-\eqref{eq:see3}.
\begin{gather}
    X = U_X \label{eq:see1} \\
    T = Bern(X) + U_T \label{eq:see2}  \\
    Y = \alpha \cdot T + \beta \cdot X + U_Y \label{eq:see3}
\end{gather}
Here, each endogenous variable has a directed edge pointing to it from an exogenous variable (e.g., $ U_X \rightarrow X $ or $ U_T $ is a parent of $ X $).
This is represented in the structural equations as being part of the structural equation (e.g., $ U_X $ in eq~\eqref{eq:see1}).
It is assumed that every endogenous variable has a parent in $ \mathbf{U} $.
We denote the parents and children of a variable, $ X $, in the graph as $ pa(X) $ and $ ch(X) $, respectively.
Ancestors of $ X $, $ an(X) $, are variables that can reach $ X $ through a directed path, and descendants of $ X $, $ de(X) $, are variables that can be reached from $ X $ through directed paths out of $ X $.

\subsubsection{Identifying effects using SCMs}\label{sec:identifying_scms}

When the SCM underlying the dataset is known, its graph $ G $ supports the use of graphical criteria for variable selection for unbiased causal effect estimation, known as adjustment set selection~\citep{shpitser-uai06,shpitser-uai10,van-uai14,perkovic-uai15,henckel-arxiv19}.
These graphical criteria rely on the concept of $d$-separation.
\begin{definition}
    ($d$-separation) A set of $ \mathbf{Z} $ nodes is said to ``block'' a path $ p $ if either
    \begin{enumerate}
        \item[$(i)$] $ p $ contains at least one arrow-emitting node that is in $ \mathbf{Z} $ or
        \item[$(ii)$] $p$ contains at least one collision node that is outside $ \mathbf{Z} $ and has no descendant in $ \mathbf{Z} $.
    \end{enumerate}
\end{definition}
If $ \mathbf{Z} $ blocks all paths from a set of nodes $ \mathbf{X} $ to a set of nodes $ \mathbf{Y} $, it is said to ``d-separate'' $ \mathbf{X} $ and $ \mathbf{Y} $, and the conditional independence, $ \mathbf{X} \perp \mathbf{Y} \mid \mathbf{Z} $, holds~\citep{bareinboim-pnas16}.
One graphical criteria for adjustment set selection is the backdoor criterion~\citep{pearl-book09}.
The backdoor criterion aims to find a set that \emph{blocks} the paths from the treatment variable to the outcome variable in the graph, starting with an arrow into the treatment variable.
Essentially, the backdoor criterion blocks all the paths going backwards out of treatment.
For example, in Figure~\ref{fig:example_scm}, there is one backdoor path, namely $ T \leftarrow X \rightarrow Y $, and so the adjustment set $ \{X\} $ sufficiently satisfies the backdoor criterion.
Several adjustment criteria have been developed, such as the adjustment criterion~\citep{shpitser-uai06,van-uai14}, generalized backdoor criterion~\citep{maathuis-annals15}, generalized adjustment criterion~\citep{perkovic-uai15}, and optimal adjustment set (O-set)~\citep{henckel-arxiv19}.
These adjustment sets defined by the various adjustment criteria can be found through \emph{identification} algorithms such as the ID algorithm~\citep{shpitser-jmlr08}, the algorithmic framework by Van der Zander et al.~\citep{van-uai14} or the pruning procedure by Henckel et al.~\citep{henckel-arxiv19}.

When an adjustment set is discovered, unbiased causal effects can be estimated by conditioning on the adjustment set.
The adjustment set can be used instead of all the features for data-driven estimation.
We discuss data-driven HTE estimation with SCMs briefly in Section~\ref{sec:hte_estimation_with_scms}.

\subsubsection{C-specific effects}

Shpitser and Pearl provide a nonparametric estimand for HTEs for some feature assignment $C=c$, called $c$-specific effects~\citep{shpitser-uai06,pearl-smr17}.
One issue with estimating $c$-specific effects is that it requires knowledge of the SCM and which variables, C, exhibit heterogeneous effects.
When the SCM is known, the O-set~\citep{henckel-arxiv19} can be used on the learned structure since it implicitly includes all potential heterogeneity-inducing variables, while some identification algorithms may miss them. When the SCM is unknown, the identifiability of \( c \)-specific effects cannot be established.
Moreover, when the set of heterogeneity-inducing variables $\mathbf{C}$ is unknown, then all combinations of features and feature values would need to be checked for heterogeneity in effect, which is exponential in the number of valid features.
Circumventing this exponential search is one of the main reasons data-driven methods are popular in practice.

Shpitser and Pearl provide a nonparametric estimand for HTEs with respect to some feature assignment $C=c$, called $c$-specific effects~\citep{shpitser-uai06,pearl-smr17}.
There are three cases when $c$-specific effects are \emph{identifiable}, i.e., there exists an unbiased estimation for $\tau(\mathbf{c}) = E[Y(1) - Y(0) | \mathbf{C}=\mathbf{c}]$.
The first case is when $\mathbf{C}$ is observed and admissible which means that it satisfies the backdoor criterion in the causal graph~\citep{pearl-book09}.
A set of variables $\mathbf{Z}$ satisfies the backdoor criterion relative to $(T,Y)$ if no node in $\mathbf{Z}$ is a descendant of $T$ and $\mathbf{Z}$ blocks every path between $T$ and $Y$ that contains an arrow into $T$.
For example, $X$ is the only variable that is admissible in all SCMs in Figure~\ref{fig:fail_examples}, since it blocks the only backdoor path.

The second case is when $\mathbf{C}$ is part of an observed admissible set.
In this case, $\mathbf{C}$ itself is not admissible, but we can observe a set of features $\mathbf{S}$ such that $\mathbf{S} \cup \mathbf{C}$ is admissible.
In Figure~\ref{fig:marloes_scm}, $\{F\}$ could induce heterogeneity but is not admissible on its own. Adding $F$ to $\{X\}$ still satisfies the backdoor criterion making $\{X, F\}$ also admissible.
The third case for identifiable $c$-specific effects is when $\mathbf{C}$ is not part of an admissible set.
The first two cases which rely on \emph{covariate adjustment}~\citep{perkovic-jmlr18} and identifying heterogeneity-inducing variables are the basis of our work.
Estimating heterogeneous effects in the third case is beyond the scope of this paper and is reserved for future work.

One issue with estimating $c$-specific effects is that it requires knowledge of the SCM and which variables exhibit heterogeneous effects.
If the SCM is unknown, the identifiability of \( c \)-specific effects cannot be established.
If $\mathbf{C}$ is unknown, then all combinations of features and feature values would need to be checked for heterogeneity in effect, which is exponential in the number of valid features.
Circumventing this exponential search is one of the main reasons why data-driven methods are popular in practice.

\subsubsection*{Detecting latent heterogeneity}

In some cases, the observed covariates that induce heterogeneity, $C$, are not known, but there is unknown heterogeneity in effects present in the data~\citep{pearl-smr17}.
Pearl discusses how heterogeneity can be detected from statistical averages alone~\citep{pearl-smr17} when the effect of treatment on the treated (ETT) and the effect of treatment on the untreated (ETU) are identifiable.

\subsubsection{HTEs in SCMs}\label{sec:htes_scms_background}

In SCMs, heterogeneous effects are represented by an \emph{interaction} or \emph{effect modification} through ancestors and parents of a variable.
In a linear structural equation, an interaction occurs as a multiplication of two variables.
For heterogeneous effects, this interaction is specifically with the treatment variable.
Consider the causal model in Figure~\ref{fig:example_scm}, where the treatment variable of interest is $ T $ and the outcome of interest is $ Y $.
There is a potential interaction between $ X $ and $ T $.
Suppose $ Y $ is generated by the following structural equation:
\begin{equation}
    Y = \alpha \cdot T + \beta \cdot X + \gamma \cdot T \cdot X + U_Y.
\end{equation}
Here, there is a fixed effect \( \alpha \) of the treatment \( T \) and a heterogeneous effect \( \gamma \) due to the interaction between \( T \) and \( X \).
Since the effect of treatment varies with another variable, $ X $, there is a heterogeneous effect.
Heterogeneity through interaction can occur with any parent of the outcome or any parent of \emph{mediating} variables~\citep{vanderweele-epid07}.
Consider another example causal model, shown in Figure~\ref{fig:multi_interaction}.
Here, there is a \emph{mediating} variable, $ D $, a variable that occurs on a path from the treatment, $ T $, to the outcome, $ Y $.
Interaction can occur with $ X, C$, or $ F $.
For example, an interaction with $ X $ and $ F $ can be shown as:
\begin{equation}
    Y = \alpha \cdot D \cdot X + \beta \cdot D \cdot F + U_Y.
\end{equation}
Since the treatment is mediated through $ D $ and $ X $ and $ F $ interact with $ D $, there will be a heterogeneous effect.
Since $ D $ is a mediator, if $ C $ interacts with $ T $ in the structural equation, then there will also be a heterogeneous effect, such as:
\begin{gather}
    D = \alpha \cdot T + \beta \cdot T \cdot C + U_D, \\
    Y = \gamma \cdot D + U_Y.
\end{gather}

\begin{figure}
    \centering
    \begin{tikzpicture}
        \node (t) at (-1.5, 0) [label={below left}:$T$,point];
        \node (x) at (-1.5, 1.5) [label={above}:$X$,point];

        \node (c) at (0.0, 1.5) [label={above}:$C$,point];
        \node (d) at (0.0, 0) [label={below}:$D$,point];

        \node (f) at (1.5, 1.5) [label={above}:$F$,point];
        \node (y) at (1.5, 0.0) [label={below}:$Y$,point];

        \path (x) edge (t);

        \path (x) edge (y);
        \path(t)edge(d);
        \path(c)edge(d);
        \path (d) edge (y);

        \path (f) edge (y);
    \end{tikzpicture}
    \caption{Causal model with multiple possible interactions that cause heterogeneous treatment effects.}\label{fig:multi_interaction}
\end{figure}

\subsection{Semi-synthetic}

The IHDP dataset is a semi-synthetic dataset, where covariates are used to generate simulated outcome so that ground truth effects are available for evaluation.
We show two results in Table~\ref{tab:ihdp_rankings}.
The first column of Table~\ref{tab:ihdp_rankings} shows the average rankings on the IHDP dataset in terms of RMSE.
From this Table, we see that FSL, TLearner, and XLearner perform consistently well in this dataset in terms of RMSE compared to all other estimator, since their rank is much higher than the next best performing estimator, SLearner.
The second column shows the $ \tau  $-risk results.
We see that the $ \tau  $-risk rank is generally correlated with the RMSE rank of each learner.
A notable difference is that Causal Forest is ranked higher based on $ \tau  $-risk, but lower based on RMSE.

\input{figures/ihdp_table.tex}
\subsection{Data-driven HTE estimation using structure learning}\label{sec:structure_learning}

As mentioned in Section~\ref{sec:hte_estimation_with_scms}, when the structural causal model is known, an adjustment set can be used to select features for HTE estimation.
However, in many real-world scenarios, the underlying causal model is unknown, and therefore, there is \emph{structure uncertainty} of the causal mechanisms.
Therefore, ignoring the underlying causal models can lead to incorrect causal effect estimations.

To use the causal model under structure uncertainty, one approach is to use causal structure learning~\citep{heinze-stat18,vowels-arxiv21}.
Many traditional causal structure learning methods can be broadly classified into constraint-based~\citep{spirtes-book00,colombo-jstor12}, score-based~\citep{chickering-jmlr02,hauser-jmlr12,nandy-annals18}, a hybrid of the two~\citep{tsamardinos-ml06}, or based on asymmetries~\citep{shimizu-jmlr06}.
These methods find a partial graph or skeleton from the data, with some unoriented edges.
With the output from structure learning methods, an adjustment set can be applied to select variables for estimation~\citep{pearl-book09,henckel-arxiv19}, and these variables are used to train an HTE estimator.

Although causal structure learning can address some problems of data-driven HTE estimation, there are some drawbacks.
First, learning the full structure can be computationally expensive for high-dimensional datasets~\citep{chickering-jmlr04,heinze-stat18,vowels-arxiv21}.
Second, causal structure learning generally cannot learn the full causal model, so learned models will have undirected edges with an ambiguous cause-effect relationship.
In addition, adjustment criteria have been developed for total causal effect estimation without heterogeneity in mind.
Therefore, adjustment sets can potentially miss heterogeneity-inducing variables, leading to incorrect HTE estimation.
In our experiments, we use randomly generated SCMs and explore how data-driven HTE estimators perform under different types of SCMs.
\section{Data driven HTE estimation}\label{sec:survey}
\section{Counterfactual prediction and effect estimation}\label{sec:survey}

This section provides a survey of existing heterogeneous treatment effect estimation methods.
We categorize the methods into two main categories: methods that perform counterfactual prediction and methods that estimate treatment effect directly.

The first category consists of methods that perform counterfactual prediction.
Counterfactual prediction methods predict counterfactual outcomes when an individual is treated and not treated~\citep{johansson-icml16,shalit-icml17,kunzel-pnas19}.
Then an estimate of an individual treatment effect is the difference in predicted outcomes.
The second category consists of methods that perform effect estimation directly.
These methods use counterfactual predictions as part of a process to learn an effect estimator~\citep{athey-pnas16,tran-aaai19,athey-annals19} or use different objective functions to estimate effects~\citep{kunzel-pnas19,louizos-neurips17}.
Table~\ref{tab:method_overview} shows subcategories of each category, together with representative papers.

\begin{table}[h]
    \centering
    \caption{An overview of HTE estimator categories and subcategories.}\label{tab:method_overview}
    \begin{tabular}[t]{M{0.47\linewidth}  M{0.47\linewidth}}
        Counterfactual Prediction & Effect Estimation \\
        \midrule
        Single-Model Approaches \citep{hill-jcgs11,kunzel-pnas19,athey-stat15, atan-aaai18} & Transformed Outcome-based Approaches \citep{rosenbaum-bio83,athey-stat15} \\
        {Two-Model Approaches} \citep{foster-statmed11,cai-bio11,athey-stat15,kunzel-pnas19, johansson-icml16,shalit-icml17}  & {Meta-Learning Approaches} \citep{xie-soc12,powers-statmed18,kunzel-pnas19,kennedy-arxiv20} \\
        {Multi-Model Approaches} \citep{schwab-arxiv19,shi-neurips19} & {Tree-based Approaches} \citep{athey-pnas16,tran-aaai19,su-jmlr09} \\
        & {Forest-based Approaches} \citep{wager-jasa18,athey-annals19,powers-statmed18} \\
        & {Deep Learning Approaches} \citep{yoon-iclr18,louizos-neurips17} \\
    \end{tabular}
\end{table}

\input{03-01-counterfactual.tex}

\input{03-02-effect_estimation.tex}

\subsection{Evaluation of HTE estimators}\label{sec:data_driven_background}

Since the true ITE is unknown in real-world scenarios, several heuristic HTE evaluation metrics have been developed for evaluating the performance of HTE estimators.
These metrics include inverse propensity weighting (IPW) validation~\citep{schuler-model18}, $ \tau $-risk~\citep{schuler-model18,nie-rlearner17}, nearest-neighbor imputation (NN-PEHE)~\citep{schwab-perfect18}, plugin validation~\citep{saito-cv19}, and counterfactual cross-validation (CFCV)~\citep{saito-cv19}.

IPW validation uses inverse propensity score weighting to a plugin estimate for the CATE, $ \tilde{\tau} $~\citep{schuler-model18}:
\begin{equation}
    \tilde{\tau}_{IPW} (\mathbf{X}_i, Y_i, T_i) = \frac{T_i}{\hat{e}(\mathbf{X}_i)}Y_i - \frac{1 - T_i}{1 - \hat{e}(\mathbf{X}_i)} Y_i.
\end{equation}
Here, $ \hat{e}(\mathbf{X}_i) $ is an estimate of the propensity score, which is the probability of receiving the treatment of interest (i.e., $\mathbb{E}[T_i|\mathbf{X}_i]$)
This plugin is used to approximate the mean-squared error for the estimator $ \hat{\tau} $:
\begin{equation}
    \hat{R}_{IPW}(\hat{\tau}) = \frac{1}{n} \sum_{i=1}^n \Big(\tilde{\tau}_{IPW} (\mathbf{X}_i, Y_i, T_i) - \hat{\tau}(\mathbf{X}_i) \Big)^2.
\end{equation}

The second metric $ \tau  $-risk~\citep{nie-rlearner17} takes as input an estimator $ \hat{\tau}  $ and approximates the objective defined in~\eqref{eq:cate_objective}.
To do this, two functions are defined: $ \hat{m}(\mathbf{X}_i) $ which is an outcome estimator for unit $i$ (i.e., $ \mathbb{E}[Y_i | \mathbf{X}_i] $) and $ \hat{e}(\mathbf{X}_i) $.
Both $ \hat{m} $ and $ \hat{e} $ are learned through supervised techniques and optimized for predictive accuracy~\citep{nie-rlearner17}.
The $ \tau  $-risk is defined as:
\begin{equation}
    \tau\text{-risk} = \frac{1}{N} \sum_{i=1}^N \Big( \big(Y_i-\hat{m}(\mathbf{X}_i)\big) - \big(T_i - \hat{e}(\mathbf{X}_i)\big)\hat{\tau}(\mathbf{X}_i) \Big)^2.
\end{equation}

The third metric, neighbor precision in estimation of heterogeneous effects (NN-PEHE)~\citep{schwab-perfect18} uses matching to approximate a treatment effect by substituting counterfactual outcomes from matched neighbors.
A unit $i$ is matched with one\footnote{Matching can be done with multiple units, and an aggregate can be computed.} control and one treated unit, denoted by NN$_0 (i)$ and NN$_1 (i)$.
Using the matched units, a substitute of the ground-truth treatment effect is computed as $ y(\text{NN}_1(i)) - y(\text{NN}_0(i)) $.
The NN-PEHE error, $ \hat{\epsilon}_{\text{NN-PEHE}} $ is computed as:
\begin{equation}
    \hat{\epsilon}_{\text{NN-PEHE}} = \frac{1}{N} \sum_{i=1}^n \bigg( \big(y(\text{NN}_1(i)) - y(\text{NN}_0(i)) \big) - \hat{\tau}(\mathbf{X}_i) \bigg)^2.
\end{equation}

Plugin validation~\citep{saito-cv19}, the fourth metric, uses machine learning algorithms to estimate a plugin $ \tau $, similar to IPW validation.
To get a plugin estimate, two predictions are needed, $ f_1(\mathbf{x}) $ and $ f_0(\mathbf{{x}}) $, which are predictions of the outcome when treated and not treated, respectively, using only the features.
Then, the plugin validation metric is defined as:
\begin{equation}
    \hat{R}_{plugin}(\hat{\tau}) = \frac{1}{N} \sum_{i=1}^N \bigg( (f_1(\mathbf{X}_i) - f_0(\mathbf{X}_i)) - \hat{\tau}(\mathbf{X}_i) \bigg)^2.
\end{equation}

Counterfactual cross-validation (CFCV) uses the plugin estimate of the CATE in a doubly robust formulation~\citep{saito-cv19}:
\begin{equation}
    \tilde{\tau}_{DR}(\mathbf{X}_i, Y_i, T_i) = \frac{T - \hat{e}(\mathbf{X}_i)}{\hat{e}(\mathbf{X}_i) (1 - \hat{e}(\mathbf{X}_i))} (Y_i - f_T(\mathbf{X}_i)) + f_1(\mathbf{X}_i) - f_0(\mathbf{X}_i),
\end{equation}
where $ \hat{e} $ is an estimate of the propensity score and $ f_t(\mathbf{x}) $ is an arbitrary regression function for a specific treatment value $ t $ that predicts the outcome given the set of features.
Then, the plugin estimate $ \tilde{\tau}_{DR} $ can be used in the mean-squared error estimate:
\begin{equation}
    \hat{R}_{CFCV} = \frac{1}{N} \sum_{i=1}^N \bigg( \tilde{\tau}_{DR}(\mathbf{X}_i, Y_i, T_i) - \hat{\tau}(\mathbf{X}_i) \bigg)^2
\end{equation}

Each metric is an estimate of the error of an HTE estimator and performance based on these heuristic metrics has been shown to correlate with true performance.
Since these metrics are heuristic, we explore their reliability in synthetic and semi-synthetic experiments in Section~\ref{sec:experiments}.

\section{Heterogeneous Treatment Effects}\label{sec:background}

In this section, we provide background on heterogeneous treatment effect estimation, starting with notations from the perspective of the potential outcome framework~\citep{rubin1974}. Then we define the problem of heterogeneous treatment effect estimation and common assumptions to allow for estimation.
Finally, we introduce structural causal models~\citep{pearl-book09} and how they can be used to assess heterogeneity.

\subsection{Data model}\label{sec:data_model}

Before formalizing the problem and data model, we introduce the notations used in this survey.
We denote a space of random variables with a calligraphic letter (e.g., $ \mathcal{X} $).
We denote random variables with a capital letter (e.g., $Y$) and random vectors with a boldfaced capital letter (e.g., $\mathbf{X}$).
Subscripts signify the random variable for a specific unit (e.g., $Y_i$ or $\mathbf{X}_i$).
Lowercase letters denote an assignment of a random variable or vector (e.g., $Y_i = y$ or $\mathbf{X}_i = \mathbf{x}$).
Each unit $i$ has a vector of ``pre-treatment'' features or covariates, $\mathbb{X}_i \in \mathcal{X}$, a binary treatment indicator, $T_i$, where $T_i=0$ and $T_i=1$ indicates $i$ is not treated and treated, respectively, and two outcomes when not treated and treated, $Y_i(T_i=0)$ and $Y_i(T_i=1)$, respectively.
The underlying data distribution is represented by a tuple $(\mathbf{X}_i, T_i, Y_i(0), Y_i(1))^N$ of $N$ independent and identically distributed individuals.
The average treatment effect (ATE) over the population is defined as the expected difference in potential outcomes:
\begin{equation}
    \label{eq:ate}
    \text{ATE} := E[Y(1) - Y(0)].
\end{equation}
The fundamental problem of causal inference is that for each individual, we can observe either the potential outcomes under treatment ($Y_i(1)$) or control ($Y_i(0)$), but not both.
Therefore, we define the observed dataset as $\mathcal{D} = (\mathbf{x}_i, t_i, y_i)^N$ with $N$ individuals where $Y_i$ is defined as:
\begin{equation}
    \label{eq:potential_outcomes}
    Y_i = T_i Y_i(1) + (1 - T_i) Y_i(0).
\end{equation}
We denote $ Y_i $ as the factual outcome and the counterfactual outcome, $ Y_i' $, as the outcome that would have occurred if the treatment was different.

\subsection{Heterogeneous Treatment Effect Estimation}\label{sec:hte}

Since each unit reacts differently to treatment, we define the individual treatment effect (ITE) as the difference in the unit's potential outcomes, denoted by $\tau_i$:
\begin{equation}
    \label{eq:ite}
    \tau_i := Y_i(1) - Y_i(0).
\end{equation}
However, due to the fundamental problem of causal inference, $\tau_i$ is not observed for any individual.
The quantity of interest is often the conditional average treatment effect (CATE):
\begin{equation}
    \tau(\mathbf{x}) := \mathbb{E} \left[ Y_i(1) - Y_i(0) \mid \mathbf{X}_i = \mathbf{x} \right].
\end{equation}
The CATE describes the expected treatment effect for units in a subpopulation with features $\mathbf{x}$.
The best estimator for CATE is also the best estimator for ITE in terms of the mean squared error~\citep{kunzel-pnas19}, though in general, the ITE is not the same as CATE~\citep{vegetabile-arxiv21}.
The goal of HTE estimation is to develop estimators of CATE which have the smallest root mean squared error for the ITE.
\begin{problem}\label{prob:hte_problem}
    (Heterogeneous treatment effect estimation)
    Given a dataset $\mathcal{D}=(\mathbf{X}_i, Y_i, T_i)$ of N instances with some true ITE. $\tau_i$, estimate CATE $\hat{\tau}(\mathbf{X})$ such that the root mean squared error (RMSE) is minimized:
    \begin{equation}\label{eq:cate_objective}
         \minl_{\hat{\tau}} \text{RMSE} = \frac{1}{N} \sum_{i=1}^N  (\tau_i - \hat{\tau}(\mathbf{X}_i)) ^ 2.
    \end{equation}
\end{problem}

\subsection{HTE Assumptions}\label{sec:hte_assumptions}

In order to estimate HTEs from observational data, some assumptions on the set of features and data given are needed.
We present some common assumptions used in heterogeneous treatment effect estimators.

\begin{assumption}\label{assumption:strong_ignorability}
    (Strong Ignorability) Given a binary treatment $T$, potential outcomes $Y(0)$, $Y(1)$, and a vector of features or covariates $\mathbf{X}$, strong ignorability states:
    \begin{align}
        \{ Y(0), Y(1) \} \perp T \mid \mathbf{X},
    \end{align}
\end{assumption}
This states the \emph{potential} outcomes are conditionally independent of the treatment conditioned on the features.
When strong ignorability holds, $\mathbf{X}$ is said to be ``admissible'', meaning treatment effects can be estimated without bias.
This assumption is also known as unconfoundedness or the assumption of no hidden confounding.

\begin{assumption}\label{assumption:overlap}
    (Population Overlap) $0 < P(T_i = 1) < 1$ with probability 1.
\end{assumption}
This assumption states that all units have a non-zero probability of assignment to each treatment condition.
Assumption~\ref{assumption:overlap} is also known as positivity or common support.

Another common assumption used is the Stable Unit Treatment Value Assumption (SUTVA) or the assumption of no interference~\citep{rubin-jasa80}.
\begin{assumption}\label{assumption:sutva}
    (Stable Unit Treatment Value Assumption (SUTVA)) Given two individuals $i, j$, the potential outcomes of $i$, $\{Y_i(0), Y_i(1) \}$ are unaffected by the treatment assignment of other individuals $\forall j$:
    \begin{equation}
        \{Y_i(0), Y_i(1)\} \perp T_j, \forall j \neq i.
    \end{equation}
    SUTVA is also often referred to as the assumption of ``no interference''.
\end{assumption}
With the SUTVA assumption, individual outcomes are not affected by treatment or outcomes of peers.
Though there has been some work focused on estimating HTEs without the SUTVA assumption~\citep{yuan-www21,tran-aaai22}, a majority of work in HTE estimation has focused on non-networked settings, which we focus on in this survey.

\subsection{Data-driven HTE Estimation}\label{sec:data_driven_background}

Given an observed dataset $\mathcal{D}$, the goal of an \emph{HTE estimator} is to estimate CATE, $\hat{\tau}$, using the available features $\mathbf{X}$ in the data.
HTE estimators estimate CATE in a variety of ways.
Some HTE estimators rely on recursive partitioning using tree-based methods~\citep{tran-aaai19,athey-pnas16,athey-annals19}.
Other HTE estimators use predictions from machine learning methods to train a causal effect estimator~\citep{kunzel-pnas19,kennedy-arxiv20} or to predict counterfactual outcomes for estimation of effect~\citep{shalit-icml17,johansson-icml16}.
We describe these methods in Section~\ref{sec:survey}.

Since the true ITE is unknown in real-world scenarios, several heuristic HTE evaluation metrics have been developed for evaluating the performance of HTE estimators.
These metrics include inverse propensity weighting (IPW) validation~\citep{schuler-model18}, $ \tau $-risk~\citep{schuler-model18,nie-rlearner17}, nearest-neighbor imputation (NN-PEHE)~\citep{schwab-perfect18}, plugin validation~\citep{saito-cv19}, and counterfactual cross-validation (CFCV)~\citep{saito-cv19}.

IPW validation uses inverse propensity score weighting to a plugin estimate for the CATE, $ \tilde{\tau} $~\citep{schuler-model18}:
\begin{equation}
    \tilde{\tau}_{IPW} (\mathbf{X}_i, Y_i, T_i) = \frac{T_i}{\hat{e}(\mathbf{X}_i)}Y_i - \frac{1 - T_i}{1 - \hat{e}(\mathbf{X}_i)} Y_i.
\end{equation}
Here, $ \hat{e}(\mathbf{x}) $ is an estimate of the propensity score (i.e., $\mathbb{E}[T_i|\mathbf{X}_i]$)
This plugin is used to approximate the mean-squared error for the estimator $ \hat{\tau} $:
\begin{equation}
    \hat{R}_{IPW}(\hat{\tau}) = \frac{1}{n} \sum_{i=1}^n \Big(\tilde{\tau}_{IPW} (\mathbf{X}_i, Y_i, T_i) - \hat{\tau}(\mathbf{X}_i) \Big)^2.
\end{equation}

The $ \tau  $-risk takes as input an estimator $ \hat{\tau}  $ and approximates the objective defined in~\eqref{eq:cate_objective}.
To do this, two functions are defined: $ \hat{m}(\mathbf{X}_i) $ which is an outcome estimator for unit $i$ (i.e., $ \mathbb{E}[Y_i | \mathbf{X}_i] $) and $ \hat{e}(\mathbf{X}_i) $ which is an estimate of the propensity score.
Both $ \hat{m} $ and $ \hat{e} $ are learned through supervised techniques and optimized for predictive accuracy~\citep{nie-rlearner17}.
The $ \tau  $-risk is defined as:
\begin{equation}
    \tau\text{-risk} = \frac{1}{N} \sum_{i=1}^N \Big( \big(Y_i-\hat{m}(\mathbf{X}_i)\big) - \big(T_i - \hat{e}(\mathbf{X}_i)\big)\hat{\tau}(\mathbf{X}_i) \Big)^2.
\end{equation}

Nearest neighbor precision in estimation of heterogeneous effects (NN-PEHE) uses matching to approximate a treatment effect by substituting counterfactual outcomes from matched neighbors~\citep{schwab-perfect18}.
A unit $i$ is matched with one\footnote{Matching can be done with multiple units, and an aggregate can be computed.} control and one treated unit, denoted by NN$_0 (i)$ and NN$_1 (i)$.
Using the matched units, a substitute of the ground-truth treatment effect is computed as $ y(\text{NN}_1(i)) - y(\text{NN}_0(i)) $.
The NN-PEHE error, $ \hat{\epsilon}_{\text{NN-PEHE}} $ is computed as:
\begin{equation}
    \hat{\epsilon}_{\text{NN-PEHE}} = \frac{1}{N} \sum_{i=1}^n \bigg( \big(y(\text{NN}_1(i)) - y(\text{NN}_0(i)) \big) - \hat{\tau}(\mathbf{X}_i) \bigg)^2.
\end{equation}

Plugin validation uses arbitrary machine learning algorithms to estimate a plugin $ \tau $, similar to IPW validation.
To get a plugin estimate, two predictions are needed, $ f_1(\mathbf{x}) $ and $ f_0(\mathbf{{x}}) $, which are predictions of the outcome when treated and not treated, respectively, using only the features.
Then, the plugin validation metric is defined as:
\begin{equation}
    \hat{R}_{plugin}(\hat{\tau}) = \frac{1}{N} \sum_{i=1}^N \bigg( (f_1(\mathbf{X}_i) - f_0(\mathbf{X}_i)) - \hat{\tau}(\mathbf{X}_i) \bigg)^2.
\end{equation}

Counterfactual cross-validation (CFCV) uses the plugin estimate of the CATE in a doubly robust formulation~\citep{saito-cv19}:
\begin{equation}
    \tilde{\tau}_{DR}(\mathbf{X}_i, Y_i, T_i) = \frac{T - \hat{e}(\mathbf{X}_i)}{\hat{e}(\mathbf{X}_i) (1 - \hat{e}(\mathbf{X}_i))} (Y_i - f_T(\mathbf{X}_i)) + f_1(\mathbf{X}_i) - f_0(\mathbf{X}_i),
\end{equation}
where $ \hat{e} $ is an estimate of the propensity score and $ f_t(\mathbf{x}) $ is an arbitrary regression function for a specific treatment value $ t $ that predicts the outcome given the set of features.
Then, the plugin estimate $ \tilde{\tau}_{DR} $ can be used in the mean-squared error estimate:
\begin{equation}
    \hat{R}_{CFCV} = \frac{1}{N} \sum_{i=1}^N \bigg( \tilde{\tau}_{DR}(\mathbf{X}_i, Y_i, T_i) - \hat{\tau}(\mathbf{X}_i) \bigg)^2
\end{equation}

Each metric is an estimate of the error of an HTE estimator and performance based on these heuristic metrics has been shown to correlate with true performance.
Since these metrics are heuristic, we explore their reliability in synthetic and semi-synthetic experiments in Section~\ref{sec:experiments}.

\subsection{Structural Causal Models}\label{sec:scm_background}

\begin{figure}
    \centering
    \newcommand{0.25}{0.5}
        \begin{subfigure}[b]{\columnwidth}
            \centering
            \begin{tikzpicture}
                \node (t) at (-0.75, 0) [label=below:$T$,point];
                \node (x) at (0.0, 0.75) [label=left:\(X\),point];
                \node(y) at (0.75, 0.0) [label=below:$Y$,point];
                \node (ux) at (0.0, 1.5) [label=above:\(U_X\),point];
                \node (ut) at (-1.5, 0.75) [label=above:\(U_T\),point];
                \node (uy) at (1.5, 0.75) [label=above:\(U_Y\),point];

                \path (x) edge (t);
                \path (x) edge (y);
                \path (t) edge (y);
                \path (ux) edge (x);
                \path (uy) edge (y);
                \path (ut) edge (t);
            \end{tikzpicture}
        \end{subfigure}
    \caption{An example causal model.
    }\label{fig:example_scm}
\end{figure}

Now we introduce an alternative to data-driven estimation of causal effects through structural causal models (SCMs) and how they relate to heterogeneous effects.

A structural causal model (SCM), $M$, consists of two sets of variables, $\mathbf{U}$ and $\mathbf{V}$, and a set of structural equations, $\mathbf{F}$, describing how values are assigned to each \emph{endogenous} variable $V^{(i)} \in \mathbf{V}$ based on the values of $\mathbf{v}$ and $\mathbf{u}$~\citep{pearl-book09}: \( v^{(i)} = f^{(i)}(\mathbf{v}, \mathbf{u}) \).
A causal graph \( G \) is a directed acyclic graph that captures the causal relationships among the variables.
Each variable, $V^{(i)}$, is represented by a node in the causal graph, an edge from $V^{(i)}$ to $V^{(j)}$ signifies that $V^{(i)}$ is a cause of $V^{(j)}$.
We also refer to the edge relationship as a parent-child relationship, where the parent is a cause of the child.
The variables in $\mathbf{U}$ are considered \emph{exogenous}, and the causal graph may omit representing $\mathbf{U}$ explicitly.
Figure~\ref{fig:example_scm} shows an example of a causal graph with three exogenous variables, $ U_X, U_T, U_Y $ and three endogenous variables, $ X, T, Y $.
An example set of associated structural equations is shown in eq~\eqref{eq:see1}-\eqref{eq:see3}.
\begin{gather}
    X = U_X \label{eq:see1} \\
    T = Bern(X) + U_T \label{eq:see2}  \\
    Y = \alpha \cdot T + \beta \cdot X + U_Y \label{eq:see3}
\end{gather}
Here, each endogenous variable has a directed edge pointing to it from an exogenous variable (e.g., $ U_X \rightarrow X $ or $ U_T $ is a parent of $ X $).
This is represented in the structural equations as being part of the structural equation (e.g., $ U_X $ in eq~\eqref{eq:see1}).
It is assumed that every endogenous variable has a parent in $ \mathbf{U} $.
We denote the parents and children of a variable, $ X $, in the graph as $ pa(X) $ and $ ch(X) $, respectively.
Ancestors of $ X $, $ an(X) $, are variables that can reach $ X $ through a directed path, and descendants of $ X $, $ de(X) $, are variables that can be reached from $ X $ through directed paths out of $ X $.

\subsubsection{Identifying effects using SCMs}\label{sec:identifying_scms}

When the SCM underlying the dataset is known, its graph $ G $ supports the use of graphical criteria for variable selection for unbiased causal effect estimation, known as adjustment set selection~\citep{shpitser-uai06,shpitser-uai10,van-uai14,perkovic-uai15,henckel-arxiv19}.
These graphical criteria rely on the concept of $d$-separation.
\begin{definition}
    ($d$-separation) A set of $ \mathbf{Z} $ nodes is said to ``block'' a path $ p $ if either
    \begin{enumerate}
        \item[$(i)$] $ p $ contains at least one arrow-emitting node that is in $ \mathbf{Z} $ or
        \item[$(ii)$] $p$ contains at least one collision node that is outside $ \mathbf{Z} $ and has no descendant in $ \mathbf{Z} $.
    \end{enumerate}
\end{definition}
If $ \mathbf{Z} $ blocks all paths from a set of nodes $ \mathbf{X} $ to a set of nodes $ \mathbf{Y} $, it is said to ``d-separate'' $ \mathbf{X} $ and $ \mathbf{Y} $, and the conditional independence, $ \mathbf{X} \perp \mathbf{Y} \mid \mathbf{Z} $, holds~\citep{bareinboim-pnas16}.
One graphical criteria for adjustment set selection is the backdoor criterion~\citep{pearl-book09}.
The backdoor criterion aims to find a set that \emph{blocks} the paths from the treatment variable to the outcome variable in the graph, starting with an arrow into the treatment variable.
Essentially, the backdoor criterion blocks all the paths going backwards out of treatment.
For example, in Figure~\ref{fig:example_scm}, there is one backdoor path, namely $ T \leftarrow X \rightarrow Y $, and so the adjustment set $ \{X\} $ sufficiently satisfies the backdoor criterion.
Several adjustment criteria have been developed, such as the adjustment criterion~\citep{shpitser-uai06,van-uai14}, generalized backdoor criterion~\citep{maathuis-annals15}, generalized adjustment criterion~\citep{perkovic-uai15}, and optimal adjustment set (O-set)~\citep{henckel-arxiv19}.
These adjustment sets defined by the various adjustment criteria can be found through \emph{identification} algorithms such as the ID algorithm~\citep{shpitser-jmlr08}, the algorithmic framework by Van der Zander et al.~\citep{van-uai14} or the pruning procedure by Henckel et al.~\citep{henckel-arxiv19}.

When an adjustment set is discovered, unbiased causal effects can be estimated by conditioning on the adjustment set.
The adjustment set can be used instead of all the features for data-driven estimation.
We discuss data-driven HTE estimation with SCMs briefly in Section~\ref{sec:hte_estimation_with_scms}.

\subsubsection{C-specific effects}

Shpitser and Pearl provide a nonparametric estimand for HTEs for some feature assignment $C=c$, called $c$-specific effects~\citep{shpitser-uai06,pearl-smr17}.
One issue with estimating $c$-specific effects is that it requires knowledge of the SCM and which variables, C, exhibit heterogeneous effects.
When the SCM is known, the O-set~\citep{henckel-arxiv19} can be used on the learned structure since it implicitly includes all potential heterogeneity-inducing variables, while some identification algorithms may miss them. When the SCM is unknown, the identifiability of \( c \)-specific effects cannot be established.
Moreover, when the set of heterogeneity-inducing variables $\mathbf{C}$ is unknown, then all combinations of features and feature values would need to be checked for heterogeneity in effect, which is exponential in the number of valid features.
Circumventing this exponential search is one of the main reasons data-driven methods are popular in practice.

Shpitser and Pearl provide a nonparametric estimand for HTEs with respect to some feature assignment $C=c$, called $c$-specific effects~\citep{shpitser-uai06,pearl-smr17}.
There are three cases when $c$-specific effects are \emph{identifiable}, i.e., there exists an unbiased estimation for $\tau(\mathbf{c}) = E[Y(1) - Y(0) | \mathbf{C}=\mathbf{c}]$.
The first case is when $\mathbf{C}$ is observed and admissible which means that it satisfies the backdoor criterion in the causal graph~\citep{pearl-book09}.
A set of variables $\mathbf{Z}$ satisfies the backdoor criterion relative to $(T,Y)$ if no node in $\mathbf{Z}$ is a descendant of $T$ and $\mathbf{Z}$ blocks every path between $T$ and $Y$ that contains an arrow into $T$.
For example, $X$ is the only variable that is admissible in all SCMs in Figure~\ref{fig:fail_examples}, since it blocks the only backdoor path.

The second case is when $\mathbf{C}$ is part of an observed admissible set.
In this case, $\mathbf{C}$ itself is not admissible, but we can observe a set of features $\mathbf{S}$ such that $\mathbf{S} \cup \mathbf{C}$ is admissible.
In Figure~\ref{fig:marloes_scm}, $\{F\}$ could induce heterogeneity but is not admissible on its own. Adding $F$ to $\{X\}$ still satisfies the backdoor criterion making $\{X, F\}$ also admissible.
The third case for identifiable $c$-specific effects is when $\mathbf{C}$ is not part of an admissible set.
The first two cases which rely on \emph{covariate adjustment}~\citep{perkovic-jmlr18} and identifying heterogeneity-inducing variables are the basis of our work.
Estimating heterogeneous effects in the third case is beyond the scope of this paper and is reserved for future work.

One issue with estimating $c$-specific effects is that it requires knowledge of the SCM and which variables exhibit heterogeneous effects.
If the SCM is unknown, the identifiability of \( c \)-specific effects cannot be established.
If $\mathbf{C}$ is unknown, then all combinations of features and feature values would need to be checked for heterogeneity in effect, which is exponential in the number of valid features.
Circumventing this exponential search is one of the main reasons why data-driven methods are popular in practice.

\subsubsection*{Detecting latent heterogeneity}

In some cases, the observed covariates that induce heterogeneity, $C$, are not known, but there is unknown heterogeneity in effects present in the data~\citep{pearl-smr17}.
Pearl discusses how heterogeneity can be detected from statistical averages alone~\citep{pearl-smr17} when the effect of treatment on the treated (ETT) and the effect of treatment on the untreated (ETU) are identifiable.

\subsubsection{HTEs in SCMs}\label{sec:htes_scms_background}

In SCMs, heterogeneous effects are represented by an \emph{interaction} or \emph{effect modification} through ancestors and parents of a variable.
In a linear structural equation, an interaction occurs as a multiplication of two variables.
For heterogeneous effects, this interaction is specifically with the treatment variable.
Consider the causal model in Figure~\ref{fig:example_scm}, where the treatment variable of interest is $ T $ and the outcome of interest is $ Y $.
There is a potential interaction between $ X $ and $ T $.
Suppose $ Y $ is generated by the following structural equation:
\begin{equation}
    Y = \alpha \cdot T + \beta \cdot X + \gamma \cdot T \cdot X + U_Y.
\end{equation}
Here, there is a fixed effect \( \alpha \) of the treatment \( T \) and a heterogeneous effect \( \gamma \) due to the interaction between \( T \) and \( X \).
Since the effect of treatment varies with another variable, $ X $, there is a heterogeneous effect.
Heterogeneity through interaction can occur with any parent of the outcome or any parent of \emph{mediating} variables~\citep{vanderweele-epid07}.
Consider another example causal model, shown in Figure~\ref{fig:multi_interaction}.
Here, there is a \emph{mediating} variable, $ D $, a variable that occurs on a path from the treatment, $ T $, to the outcome, $ Y $.
Interaction can occur with $ X, C$, or $ F $.
For example, an interaction with $ X $ and $ F $ can be shown as:
\begin{equation}
    Y = \alpha \cdot D \cdot X + \beta \cdot D \cdot F + U_Y.
\end{equation}
Since the treatment is mediated through $ D $ and $ X $ and $ F $ interact with $ D $, there will be a heterogeneous effect.
Since $ D $ is a mediator, if $ C $ interacts with $ T $ in the structural equation, then there will also be a heterogeneous effect, such as:
\begin{gather}
    D = \alpha \cdot T + \beta \cdot T \cdot C + U_D, \\
    Y = \gamma \cdot D + U_Y.
\end{gather}

\begin{figure}
    \centering
    \begin{tikzpicture}
        \node (t) at (-1.5, 0) [label={below left}:$T$,point];
        \node (x) at (-1.5, 1.5) [label={above}:$X$,point];

        \node (c) at (0.0, 1.5) [label={above}:$C$,point];
        \node (d) at (0.0, 0) [label={below}:$D$,point];

        \node (f) at (1.5, 1.5) [label={above}:$F$,point];
        \node (y) at (1.5, 0.0) [label={below}:$Y$,point];

        \path (x) edge (t);

        \path (x) edge (y);
        \path(t)edge(d);
        \path(c)edge(d);
        \path (d) edge (y);

        \path (f) edge (y);
    \end{tikzpicture}
    \caption{Causal model with multiple possible interactions that cause heterogeneous treatment effects.}\label{fig:multi_interaction}
\end{figure}