text:
\section{Introduction}  \label{sec:intro}
	The application of advanced statistical learning tools in causal inference has gained popularity in recent years, partly due to the fact that large datasets are becoming available at relatively lower costs (thanks to e.g.~electronic health records, social network data etc.). One of the increasingly common objectives of causal inference in many disciplines is to draw inferences about individual-level treatment effects, as opposed to inferring treatment effects on average across the entire population. The importance of inferring individual-level treatment effects lies in the fact that treatment effects are very often heterogeneous across units of analysis. Two such examples arise in precision medicine and personalized advertisement, where the ultimate goal is to make decisions at the level of an individual patient or user. For instance, patients with high cholesterol levels respond differently to statin prescriptions based on their demographics. This level of analysis requires causal inference methods that can accurately predict the impact of treatment, as well as quantify its uncertainty, at a fine resolution. To this end, popular statistical learning algorithms that exhibit excellent predictive performance, such as tree ensembles, kernel methods and neural networks, can be exploited also, with due adjustments, in causal settings. This paper is therefore motivated by the growing number of non-parametric regression methods for modelling \emph{Individual Treatment Effects} (ITE) using large datasets to answer individual-level questions, and reviews some of the most popular ones laying down the underlying assumptions and pros/cons for each of them.

	Regardless of the precise causal questions of interest, the fundamental challenge in causal inference is that the quantity of interest  --- namely the effect of treatment --- is an unknown model's parameter. Moreover, treatment effect and selection into treatment are invariably entangled. As a result, advanced statistical learning tools, capable of exploiting large datasets in supervised learning settings, cannot be directly applied. Instead, the latent treatment effect is inferred by reconstructing counterfactual statements either through sampling (randomization in the administration of treatment) and/or by adjusting for covariates that affect both the outcome of interest and treatment assignment, and by doing so ``confound'' the treatment effect (confounding).

	Randomized Controlled Trials (RCTs), where treatment is randomly allocated, mar\-ginally on confounding factors (such as medical or socio-economic characteristics), are considered to be the gold standard for causal inference. However, fully randomized studies are costly, difficult to access, and sometimes suffer from problems such as non-compliance and other issues of missingness not at random that might invalidate the randomization mechanism, or external validity of the results.

	In contrast, data of observational nature, where treatment administration is not randomized, are more easily accessible and abundantly present in many applied fields. However, observational studies present several drawbacks, largely attributable to three complex phenomena. The first is \textit{selection bias}, which manifests when the treatment allocation mechanism is not under the researcher's control, but determined by observable or unobservable factors. This constitutes a potential source of confounding that needs to be controlled for, as it generates structural differences between the treated and the control groups. The second phenomenon is known as \textit{partial overlap}, which occurs when there are regions in the space of relevant covariates where only treated, or only control, units are present. As a result, units in these particular regions lack an appropriate comparator with similar demographics. These two issues are closely related, as partial overlap may be a direct consequence of selection bias. The third phenomenon relates to the fact that treatment allocations and their corresponding outcomes may not be independent across individuals.

	In this paper we review some of the most recent non-parametric regression methods for estimating heterogeneous effects arising from a binary treatment assignment, at a individual level, and introduce a taxonomy that classifies these methods under the same unified framework. More in details, we provide an overview of the implied assumptions and compare methods' performance with the help of two simulated studies. In addition, we illustrate practical application of some of the methods on a real-world dataset to investigate the effect of participation to school meal programs on students' health indicators, which is likely to display heterogeneity patterns depending on students' demographics. We focus on methods which circumvent the potential issues stemming from lack of complete (and controlled) randomization by making a set of assumptions which are common in the literature. The first assumption is \emph{unconfoundedness}, which ensures that there is no unobserved confounding factor driving selection into treatment. Unconfoudedness might be, in some cases, a strained assumption to make, but represents less of a threat in settings where a large number of relevant covariates is available. The second assumption is \emph{common support}, which states that each unit, identified by a given set of demographics, has non-zero probability of being observed in each of the treatment groups. In other words, common support ensures that there is no deterministic component in treatment administration. Finally, the last assumption is known as \emph{Stable Unit Treatment Value Assumption} (SUTVA), and states that the response to treatment of one unit is not affected by other units' assignment to treatment, thus ensuring that there is no interference.

	The rest of the work is organized as follows. Section \ref{sec:ProblemSetup} formulates the problem of estimating \emph{Conditional Average Treatment Effect} (CATE) by using the potential outcomes representation by \cite{rubin1978}. Section \ref{sec:cate} reviews the most recent methods of estimating CATE and the problem of model selection. Section \ref{sec:simu} presents results from two different simulated studies that we conducted to compare performance of some of the methods. Section \ref{sec:realdata} provides an example of a real-world social sciences application. Section \ref{sec:discussion} concludes with a discussion.

	\section{Problem setup}\DIFaddbegin \label{sec:ProblemSetup}
	\DIFaddend

	We will follow the Neyman-Rubin Causal Model, outlined in \cite{rubin1978} and \cite{imbens_rubin_2015}, which conceives causal inference as a missing data problem. For each unit of analysis $i \in \{ 1, ..., N \}$, given a binary treatment assignment $Z_i \in \{ 0,1 \}$, where $Z_i = 1$ indicates exposure to the treatment and $Z_i = 0$ indicates no exposure, the framework defines the quantities $\big( Y_i^{(0)}, Y_i^{(1)} \big)$ as potential outcomes. $Y_i^{(1)}$ corresponds to unit $i$'s outcome under exposure to the treatment, while $Y_i^{(0)}$ corresponds to its outcome under no exposure, and only one of them is actually observed. We will consider, throughout this work, a setting where the outcome variable of interest is continuous, i.e.,~$\big( Y_i^{(0)}, Y_i^{(1)} \big) \in \mathbb{R}^2$, but most of the ideas and methods presented can be generalized. It is worth pointing out that the Neyman-Rubin framework is not the only viable approach to causal inference. For example, \cite{dawid_2000} develops a critique of potential outcomes and proposes a different framework based on Bayesian decision analysis.

	Given the framework outlined above, in case both potential outcomes were observed, one would have access also the \emph{Individual Treatment Effects} (ITE), defined as $Y_i^{(1)} - Y_i^{(0)}$. However, as described earlier, the fundamental problem of causal inference is that, for each individual $i$, only one of the two potential outcomes $\big( Y_i^{(0)}, Y_i^{(1)} \big) \in \mathbb{R}^2$ is observed, corresponding to the realized treatment assignment: $Y_i = Z_i Y_i^{(1)} + (1 - Z_i) Y_i^{(0)}$. This implies that $Y_i^{(1)} - Y_i^{(0)}$ cannot be estimated as in a usual regression problem.

	Given a dataset of either observational or randomized nature $\mathcal{D}_i = \{ \bm{X}_i, Z_i, Y_i \} $, with $i \in \{ 1, ..., N \}$, where $\bm{X}_i \in \mathcal{X}$ denotes a $d$-dimensional set of observed covariates for the individual $i$ which are potential source of confounding to be controlled for, the aim is to estimate Conditional Average Treatment Effect (CATE), defined as
	\begin{equation} \label{CATE}
	\tau(\bm{x}_i) = \mathbb{E} \Big[ Y_i^{(1)} - Y_i^{(0)} \mid \bm{X}_i = \bm{x}_i \Big]~.
	\end{equation}

	\begin{figure}[t]
		\centering
		\begin{minipage}{.5\textwidth}
			\centering
			\includegraphics[scale=0.52]{POs.pdf}
		\end{minipage}%
		\begin{minipage}{.4\textwidth}
			\centering
			\includegraphics[scale=0.52]{ITE.pdf}
		\end{minipage}
		\caption{\small Simulated example with one single covariate $X$. Left panel: observed outcomes for the treated (blue dots) and control (red dots) groups with underlying conditional mean function (dashed lines), and unobservable counterfactual outcomes (grey dots). Right panel: unobservable true ITE (grey dots) and corresponding conditional mean function (dashed line).}
		\label{example}
	\end{figure}

	\noindent The two quantities $\mu_Z (\bm{x}_i) = \mathbb{E} \big[ Y^{(Z_i)} \mid \bm{X}_i = \bm{x}_i \big]$ in (\ref{CATE}) are the conditional average potential outcomes. The intuition behind the estimation of $\tau(\bm{x}_i)$ is the following. In case both potential outcomes were observed, then $Y_i^{(1)} - Y_i^{(0)}$ (ITE) would be modelled as the response variable in a regression framework where $\bm{X}_i$ are the $d$ regressors, and where the aim is to estimate the conditional mean of the outcome, namely $\tau(\bm{x}_i) = \mathbb{E} [ Y_i^{(1)} - Y_i^{(0)} \mid \bm{X}_i = \bm{x}_i ]$. Figure \ref{example} provides a graphical representation of a simple single-covariate example, where coloured dots show observed values $Y_i = Y^{(Z_i)}$ of the response, and grey their corresponding (unobservable) counterfactuals $Y_i^{(1-Z_i)}$. Notice that the example is purely illustrative and serves as visual aid to introduce the reader to the key concepts in the Rubin-Neyman framework.

	An additional quantity of interest under the Neyman-Rubin framework is the \emph{propensity score}, which is defined, for each unit of analysis $i$, as the probability of being selected into treatment, given a set of observed covariates that we denote by $\tilde{\bm{X}}_i = \tilde{\bm{x}}_i$, to stress the fact that it might be different from the set of covariates used for the estimation of $\tau(\bm{x}_i)$; more formally:
	\begin{equation*}
	\pi(\bm{x}_i) = \mathbb{P} ( Z_i = 1 \mid \tilde{\bm{X}}_i = \tilde{\bm{x}}_i ) ~ .
	\end{equation*}
	Thus, for each unit, the binary treatment assignment $Z_i$ can be seen as the outcome of a Bernoulli experiment where $Z_i \sim \text{Bern} \big( \pi(\bm{x}_i) \big)$. In the simple example of Figure \ref{example}, the propensity score is generated as a monotone function of the covariate $X$. This is why treated units are more frequently observed for higher values of $X$, while control units are more frequent for lower values of $X$. The propensity score distribution in this case is also slightly skewed to the right; as a consequence, there are more units in the control group compared to the treated one.

	Under the potential outcomes framework, the untestable unconfoundedness assumption is equivalent to assuming the conditional independence
	\begin{equation} \label{CIA}
	\big( Y_i^{(0)}, Y_i^{(1)} \big) \perp \!\!\! \perp Z_i \mid \bm{X}_i  ~ ,
	\end{equation}
	which ensures that selection is determined by observed demographics. A well-known result concerning the propensity score, derived by \cite{rosenbaum_1983}, is that if (\ref{CIA}) holds, then
	\begin{equation} \label{CIAps}
	\big(Y_i^{(0)}, Y_i^{(1)} \big) \perp \!\!\! \perp Z_i \mid \pi(\bm{X}_i)~.
	\end{equation}
	Conditional independence in (\ref{CIAps}) represents a different way of expressing unconfoundedness, by conditioning only on the propensity score rather than on the full set of covariates. In the context where the number of covariates is high, $\pi(\bm{X}_i)$ represents a 1-dimensional summary of a $d$-dimensional covariate set $\bm{X}_i \in \mathcal{X}$ which achieves conditional independence between outcome and treatment assignment. There are many different ways of making effective use of propensity score for deriving population-level or individual-level causal estimators. In some approaches propensity score assumes a central role, such as in matching and blocking methods \citep{rosenbaum_1984, rosenbaum_1985}, compared to regression-based methodologies that we focus on in this work. \cite{imbens2004} offers a nice and thorough overview of methods based on the propensity score, and also of other methods (e.g.~regression and matching).

	The assumption of common support, under the Rubin-Neyman framework, implies that propensity score is such that $0 < \pi(\bm{x}_i) < 1$, for each $i$, which equivalently means that treatment assignment is not deterministic, and that an individual with covariate values $\bm{X}_i = \bm{x}_i$ can be potentially observed in either treatment group with non-zero probability.

	Under unconfoundedness and common support assumptions, CATE can be estimated from observed data in $\mathcal{D}_i = \{ \bm{X}_i, Z_i, Y_i \}\;i \in \{ 1, ..., N \} $. In fact, under unconfoudedness, the conditional average potential outcomes are such that
	\begin{equation}  \label{mu_ide}
	\begin{split}
	\mu_{Z} (\bm{x}_i) & = \mathbb{E} \left[ Y^{(Z_i)} \mid \bm{X}_i = \bm{x}_i \right] = \mathbb{E} \left[ Y^{(Z_i)} \mid Z_i = z_i, \bm{X}_i = \bm{x}_i \right] \\
	& = \mathbb{E} \left[ Y_i \mid Z_i = z_i, \bm{X}_i = \bm{x}_i \right]~,
	\end{split}
	\end{equation}
	for $z_i \in \{ 0, 1 \}$, where the second equality generates from the conditional independence between $\big( Y_i^{(0)}, Y_i^{(1)} \big)$ and $Z_i$ given $\bm{X}_i$, while the last one from the identity $Y_i = Z_i Y_i^{(1)} + (1 - Z_i) Y_i^{(0)}$. Hence, as a straightforward implication of (\ref{mu_ide}), one can derive an estimator for CATE as
	\begin{equation} \label{CATEide}
	\begin{split}
	\tau(\bm{x}_i) & = \mathbb{E} \left[ Y_i^{(1)} - Y_i^{(0)} \mid \bm{X}_i = \bm{x}_i \right]
	\\
	& = \mathbb{E} \left[ Y_i^{(1)} \mid \bm{X}_i = \bm{x}_i \right] - \mathbb{E} \left[ Y_i^{(0)} \mid \bm{X}_i = \bm{x}_i \right]
	\\
	& = \mathbb{E} \left[ Y_i^{(1)}  \mid Z_i = 1 , \bm{X}_i = \bm{x}_i \right] - \mathbb{E} \left[ Y_i^{(0)} \mid Z_i = 0 , \bm{X}_i = \bm{x}_i \right]
	\\
	& = \mathbb{E} \left[ Y_i \mid Z_i = 1 , \bm{X}_i = \bm{x}_i \right] - \mathbb{E} \left[ Y_i \mid Z_i = 0 , \bm{X}_i = \bm{x}_i \right]~,
	\end{split}
	\end{equation}
	where, as in (\ref{mu_ide}), the third equality is given by unconfoundedness, and the last one by the identity $Y_i = Z_i Y_i^{(1)} + (1 - Z_i) Y_i^{(0)}$. Common support assumption is then needed to guarantee that the two conditional average potential outcomes $\mu_{Z} (\bm{x}_i) = \mathbb{E} \left[ Y^{(Z_i)} \mid \bm{X}_i = \bm{x}_i \right]$ exist for each values $z_i$ and $\bm{x}_i$ in their supports, and thus can be estimated through the observed quantities in the conditional expectations $\mathbb{E} \left[ Y_i \mid Z_i = z_i, \bm{X}_i = \bm{x}_i \right]$. To clarify, suppose that common support does not hold for $\bm{X}_i = x^*$ and that $\pi (x^*) = 0$ (without loss of generality), then the conditional average potential outcome $\mu_1 (x^*) = \mathbb{E} \big[ Y^{(1)} \mid \bm{X}_i = x^* \big]$ does not theoretically exist, and it would not make sense to attempt to even estimate it.

	\subsection{Non-Parametric regression framework for CATE estimation}

	Throughout this work, we will review non-parametric regression approaches where the outcome surface $Y_i$ is modelled as a function of $( \bm{X}_i$, $Z_i )$ and some unobservable error term $\varepsilon_i$~: $Y_i = g (\bm{X}_i, Z_i, \varepsilon_i)$. More specifically, the reviewed methods generally assume that the error term is additive and normally distributed, which leads to the following setup:
	\begin{equation} \label{onestage}
	Y_i = f(\bm{X}_i, Z_i) + \varepsilon_i~, \qquad \text{where} \quad \varepsilon_i \sim \mathcal{N} (0,\sigma^2)~,
	\end{equation}
	where $f(\bm{X}_i, Z_i) = \mathbb{E} \big[ Y_i \mid \bm{X}_i, Z_i \big]$ is learnt from the data. The strength of non-parametric regression models is their robustness to misspecification of the functional form of $f(\cdot)$ (e.g.~tree-based methods model $f(\cdot)$ as piece-wise constant, splines as piece-wise polynomial, etc.). The covariates $\bm{X}_i \in \mathcal{X}$ represent a potential source of confounding, in that they might affect both the treatment and the outcome, and are thus controlled for. In a setting where the number of available covariates is high, one might want to resort to regularization in the estimation of $f(\cdot)$; this is more likely to be of interest in observational studies in socio-economic disciplines rather than clinical research, where closely controlling for which covariates are included in the model is preferable. However, as explained in both \cite{hahn_2018} and \cite{hahn_2020}, regularization should be handled carefully in this context; we will return to this point in Section \ref{sec:taulearner}.

	In the regression setup illustrated in (\ref{onestage}), some of the methods reviewed in the next section reserve a specific role for the propensity score. For the remaining methods which do not explicitly use the propensity score, in the simulated studies that we conduct in Section \ref{sec:simu}, we follow the suggestion of \cite{hahn_2020} and incorporate it via the following two-stage regression framework:
	\begin{align} \label{twostage}
	\begin{split}
	\pi(\bm{X}_i) = & ~ \mathbb{P} \big( Z_i = 1 \mid \bm{X}_i \big)   \\
	Y_i = & ~ f \Big( \big[ \bm{X}_i ~~ \pi(\bm{X}_i) \big] , Z_i \Big) + \varepsilon_i ~ .
	\end{split}
	\end{align}
	The first stage in (\ref{twostage}) involves estimating the propensity score, while the second embeds it as an extra covariate in the covariate set. Any probabilistic classifier is suitable for use in the first stage regression (e.g.~logistic regression, probit Bayesian Additive Regression Trees, Gaussian Process Classifier, etc.). As explained in \cite{hahn_2020}, and as we will describe later in Section \ref{sec:taulearner}, the addition of $\pi(\bm{X}_i)$ to the covariate set in (\ref{twostage}) represents an effective way to tackle bias deriving from \emph{targeted selection}. Targeted selection arises when individuals are selected into treatment based on the prediction of otherwise adverse outcome, i.e.,~when $\pi(\bm{X}_i)$ is a strictly monotone function of $\mathbb{E} [ Y_i^{(0)} \mid \bm{X}_i = \bm{x}_i ]$, and is common in many observational studies (e.g.~medical or socio-economic studies).

	CATE estimators can be derived from the representations in (\ref{onestage}) and (\ref{twostage}). There are a few different approaches for deriving an estimator for CATE from (\ref{onestage}) and (\ref{twostage}), that will be analyzed in the next section.

	\section{Estimating CATE} \label{sec:cate}
	Given the framework outlined in the previous section, various meta-algorithms designed to derive a CATE estimator have been proposed in the literature. These meta-algorithms are often referred to as ``Meta-Learners'', in that they are built on top of ``base-learners'', which are common machine learning algorithms (e.g.~tree ensembles, regularized regressions, neural networks etc.). We attempt to build a taxonomy of these ``Meta-Learners'' approaches in Section \ref{sec:metalearners}, while in Section \ref{sec:modelsel} we present a review of the metrics used for model selection when estimating CATE, which is a substantially hard, arguably impossible, problem.

	\subsection{Meta-Learners} \label{sec:metalearners}

	As mentioned in the earlier sections, we will partly build on top of the review of \cite{kunzel_2017}, and we will illustrate the most recent contributions stemming from both statistics and computer science literature. A concise summary of the presented ``Meta-Learners'', together with the relevant references, can be found in Table \ref{metalearners}.

	\begin{table}
		\caption{Summary of meta-learners described in this paper}
		\footnotesize
		\centering
		\begin{tabular}{>{\small\arraybackslash}m{3.3cm} | >{\small\arraybackslash}m{4cm} | >{\small\arraybackslash}m{6.1cm}}
			&   \textbf{References}   &    \textbf{CATE estimator}
			\\
			\midrule

			\textbf{S-Learner} &  \cite{hill_2011, foster_2011}  &   $\tau (\bm{x}_i) = f(\bm{x}_i, 1) - f(\bm{x}_i, 0)$
			\\
			\midrule

			\textbf{T-Learner} &  \cite{athey_2016, lu_sadiq_2018}, \cite{powers_2018}  &   $\tau (\bm{x}_i) = f_1(\bm{x}_i) - f_0(\bm{x}_i)$
			\\
			\midrule

			\textbf{X-Learner} &  \cite{kunzel_2017}  &  $ \tau(\bm{x}_i) = \pi(\bm{x}_i) \tau_0(\bm{x}_i) + \big(1 - \pi(\bm{x}_i) \big) \tau_1(\bm{x}_i) $
			\\
			\midrule

			\textbf{R-Learner} &  \cite{wager_2019}  &  $  \tau(\bm{x}_i) = \argmin_{\tau} \Big \{  L_n \big( \tau(\cdot) \big)  + \Lambda_n \big(  \tau(\cdot) \big)  \Big \}  $
			\\
			\midrule
			\textbf{Multitask-Learner} &  \cite{vanderschaar_2017, vanderschaar_2018}  &   $\tau(\bm{x}_i) = \mathbf{f^{\top}} (\bm{x}_i) \mathbf{e} $
			\\
			\midrule
			$\tau$\textbf{-Learner} &  \cite{hahn_2020}  & $\tau (\bm{x}_i)$ as explicit model parameter
			\label{metalearners}
		\end{tabular}
	\end{table}

	\subsubsection{S-Learners}

	``Single-Learners", shortened to S-Learners, have been implicitly proposed in two early contributions \citep{hill_2011, foster_2011} and derive an estimator for CATE by using treatment assignment as ``just another covariate" in the covariate space $\mathcal{X}$, which means that CATE is estimated as
	\begin{equation}
	\tau (\bm{x}_i) = f(\bm{x}_i, 1) - f(\bm{x}_i, 0)~.
	\end{equation}
	An S-Learner fits a single surface $f(\cdot)$ through a base-learner and derives CATE estimates by taking the difference between the two conditional average potential outcomes, which are represented by the fitted $\widehat{f} (\cdot)$ with $Z_i = 1$ and $Z_i = 0$ respectively. The underlying assumption is that the group-specific conditional average potential outcomes stem from the same model, with conditional mean function $f (\cdot)$ and error term $\varepsilon_i$. Regression trees are popular base-learners employed in the context of S-Learners. For instance, \cite{hill_2011} advocates the use of Bayesian Additive Regression Trees (BART), while \cite{foster_2011} of random forests.

	The left panel plot of Figure \ref{STBART} shows a S-Learner BART fit for the conditional mean $\widehat{f} (\cdot)$ of the single-covariate simulated example already encountered in Figure \ref{example}. Notice that the dashed line representing $\widehat{f} (\cdot)$ has a unique color (grey) to emphasize the fact that S-Learner fits a unique surface.

	Since a S-Learner fits a single regression, it is quite restrictive in the sense that it does not incorporate the way in which $f(\cdot)$ changes with respect to $Z$; and it does not take into account the fact that the distribution of the covariates $\bm{X}_i \in \mathcal{X}$ may vary in $Z$, as a result of selection bias. This is obviously more problematic when working with observational data, while it is less so with randomized studies where selection bias is less likely to be present. \cite{vanderschaar_2018} and \cite{hahn_2020} have both identified that the main drawback of S-Learners is their lack of ability in adapting to variable levels of sparsity and smoothness across the two treatment groups, since they impose the same regularizing conditions for both treated and control groups. A S-Learner will then perform poorly in a situation where the outcome surface complexity is very different across the two groups. On the contrary, S-Learner will perform well when CATE is of small magnitude, as $\tau (\bm{x}_i)$ surface complexity does not heavily depend on $Z_i$, or, in other words, does not vary much across treatment groups. For example, consider the case of a S-Learner employing a tree ensemble base-learner, such as BART. Since a tree ensemble method like BART picks splitting variables at each node in each tree randomly, it might not even choose $Z$ as splitting variable in some of the trees in the ensemble, so that $Z$ will possibly be included in most of the trees fitting the response $Y$, but not necessarily in all of them. The exclusion of $Z$ from the splitting rules of a tree is more likely to happen as the number of covariates $\bm{X}_i$ grows larger, in that the model has a larger set of splitting variables to pick from. This intuitively explains why S-Learners turn out to be appropriate in situations where the magnitude of ground-truth CATE is rather small, relative to the variation in outcome attributable to the covariates. It may happen in real world applications, e.g.~in clinical studies, that the researcher has a prior idea regarding the magnitude of the treatment effect. As it becomes clearer in the following sections, this is a non-negligible piece of information when choosing a method to perform the analysis.

	\begin{figure}[t]
	\centering
	\begin{minipage}{.5\textwidth}
		\centering
		\includegraphics[scale=0.52]{SBART.pdf}
	\end{minipage}%
	\begin{minipage}{.5\textwidth}
		\centering
		\includegraphics[scale=0.52]{TBART.pdf}
	\end{minipage}
	\caption{\small Simulated one-covariate data from Section \ref{sec:ProblemSetup}. Left panel: conditional mean fit from a S-Learner BART (dashed grey line). Right panel: conditional mean fit from a T-Learner BART (blue and red dashed lines).}
	\label{STBART}
\end{figure}

	\subsubsection{T-Learners}

	``Two-Learners'', shortened to T-Learners, derive an estimator for CATE by fitting two separate surfaces for the treated and control groups and computing their difference:
	\begin{equation}
	\tau (\bm{x}_i) = f_1(\bm{x}_i) - f_0(\bm{x}_i)~.
	\end{equation}
	Versions of T-Learners can be found in many contributions in the literature. For instance, \cite{athey_2016}, \cite{lu_sadiq_2018} and \cite{powers_2018} offer some examples employing decision trees, random forests and gradient boosted trees as base-learners respectively.
	In contrast to S-Learners, T-Learners separate the two treatment groups when modelling response variable $Y$, and assume that group-specific conditional average potential outcomes are derived from separate regression models, characterized by different conditional means $f_1(\cdot)$ and $f_0(\cdot)$ and independent error terms $\varepsilon_{1i}$ and $\varepsilon_{0i}$. This allows to preserve distributional differences across the two groups that might originate from selection bias, and to take into account different degrees of sparsity and smoothness that vary with $Z$, when regressing $Y$ against $X$. On the other hand, a shortcoming of T-Learners is that, as a result of splitting the sample in two, they do not allow sharing common underlying information between the groups when estimating the two surfaces. This is particularly not ideal in a scenario where individuals in the two groups share the same distributional characteristics, such as in a randomized study.

	The right panel plot in Figure \ref{STBART} displays a T-Learner BART fit for the conditional means of the two treatment groups $f_1 (\cdot)$ and $f_0 (\cdot)$, on the same one-covariate simulated example of Figure \ref{example}. Notice that fitted $\widehat{f}_1 (\cdot)$ and $\widehat{f}_0 (\cdot)$ are differentiated by colors (blue and red dashed lines respectively), emphasizing the fact that T-Learner fits two separate surfaces.

	T-Learners work particularly well when complexity of the response surface is very different across treatment groups, and so CATE itself turns out to be a rather complex function. In addition, as formally derived by \cite{vanderschaar_2018}, T-Learners are expected to do well as sample size $N$ goes to infinity, and more observations per group are available to estimate $f_1 (\cdot)$ and $f_0 (\cdot)$. However, this is not usually the case with real world data. For example, in presence of imbalanced designs, where one group is larger than the other, splitting the sample in two would leave few observations for the estimation of $f_{Z}$ in the smaller group. On the contrary, T-Learners tend to perform quite poorly in settings where CATE surface is relatively simple and heterogeneity patterns are not so sophisticated, i.e.,~situations where S-Learner usually performs better. Hence if subject-matter knowledge suggests that treatment impact is likely to be of significant magnitude, a T-Learner might be the preferred choice.

	\subsubsection{X-Learners}

	X-Learners have been proposed by \cite{kunzel_2017} as an extension of T-Learners, and derive a CATE estimator in three steps. In the first step, conditional average potential outcomes are fitted as in a T-Learner approach, that is by using two separate regression models for the conditional means  $f_1(\bm{x}_i)$ and $f_0(\bm{x}_i)$. This carries the same underlying assumptions, benefits and drawbacks of T-Learners, illustrated in the previous section. Then in the second step, ``imputed treatment effects" are computed for each group separately; these are defined as the differences between the group-specific observed outcome $Y_i^Z$, and the estimated unobservable conditional average potential outcome $\widehat{Y}_i^{(Z)}$ derived in the first step, more formally:
	\begin{equation}  \label{impTE}
	\begin{split}
	\tilde{D}^1_i & = Y^1_i - \widehat{Y}_i^{(0)} \qquad \text{if} \quad Z_i = 1 \\
	\tilde{D}^0_i & = \widehat{Y}_i^{(1)} - Y^0_i \qquad \text{if} \quad Z_i = 0  ~.
	\end{split}
	\end{equation}
	The second step thus attempts to recover the unobservable differences $ D_i =Y_i^{(1)} - Y_i^{(0)}$ (ITE) by replacing the unobservable potential outcomes with the relative conditional average potential outcome estimates $\widehat{Y}_i^{(Z)}$, separately for the treated and control group.
	In the last step, $\tilde{D}^1_i$ and $\tilde{D}^0_i$ are used as response variables in two separate regressions, employing the chosen base-learner (linear regression, random forest, BART, etc.), to obtain estimates of $\hat{\tau}_1(\bm{x}_i)$ and $\hat{\tau}_0(\bm{x}_i)$, using covariates $\bm{X}_i$ as regressors. These two independent regressions can be depicted as
	\begin{equation} \label{lastXstep}
	\begin{split}
	\tilde{D}^1_i & = \tau_1(\bm{X}_i) + \eta_{1i}  \qquad \text{if} \quad Z_i = 1 \\
	\tilde{D}^0_i & = \tau_0(\bm{X}_i) + \eta_{0i}  \qquad \text{if} \quad Z_i = 0   ~ ,
	\end{split}
	\end{equation}
	where the two estimated quantities $\hat{\tau}_1(\bm{x}_i)$ and $\hat{\tau}_0(\bm{x}_i)$ from (\ref{lastXstep}) are group-specific CATE estimates. The final CATE estimate is then obtained through a weighted average of the two group-specific CATE estimates,
	\begin{equation} \label{Xweight}
	\widehat{\tau}(\bm{x}_i) = g(\bm{x}_i) \widehat{\tau}_0(\bm{x}_i) + \big(1 - g(\bm{x}_i) \big) \widehat{\tau}_1(\bm{x}_i)  ~ ,
	\end{equation}
	where $g(\bm{x}_i) \in [0,1]$ is a given weighting function. The authors propose to set $g(\cdot)$ equal to a propensity score estimate $g(\bm{x}_i) = \widehat{\pi}(\bm{x}_i)$, but it can take different forms (e.g.~$g(\bm{x}_i) = 1$ or $g(\bm{x}_i) = 0$).

	The intuition behind the X-Learner approach is the following. If both potential outcomes were observable for each individual $i$, then $\mathbb{E} \big[ Y_i^{(1)} - Y_i^{(0)} \mid \bm{X}_i = \bm{x}_i \big]$ (CATE) would be estimated by regressing $ Y_i^{(1)} - Y_i^{(0)} $ (ITE) directly on the covariates $\bm{X}_i$. However, since the difference $ Y_i^{(1)} - Y_i^{(0)} $ is not observable, one can estimate it by using the observed potential outcome $Y_i^Z$ and replacing the counterfactual unobservable outcome $\widehat{Y}_i^{(1-Z)}$ with an estimate, derived via T-Learner in the first step.

	The procedure laid out by the X-Learner in the first two steps is comparable to that of a T-Learner as both methods imply dividing the sample into the two groups. The main difference lies then in the last additional step, where the X-Learner attempts to re-balance CATE estimates through weighting via the estimated propensity score. X-Learners are in fact shown to perform well in unbalanced settings, where one of the treatment groups (typically the control group) is larger in size compared to the other one, and where T-Learner might face issues. A final note about the X-Learner balancing step is that it requires careful specification of the propensity score model, as poor PS estimates would defy the re-balancing purpose. Propensity score estimation is no different from a standard prediction problem and boils down to performing a probabilistic binary outcome regression of $Z$ on $\bm{X}$. We will see that also other Meta-Learners envisage an important balancing role for propensity score, and thus must deal with its estimation.

	\begin{figure}[t]
		\centering
		\begin{minipage}{.5\textwidth}
			\centering
			\includegraphics[scale=0.52]{XBART1.pdf}
		\end{minipage}%
		\begin{minipage}{.5\textwidth}
			\centering
			\includegraphics[scale=0.52]{XBART2.pdf}
		\end{minipage}
		\caption{\small X-Learner BART applied to the simulated one-covariate example. Left panel: unobservable ITE (grey dots) and imputed treatment effects $D^1$ and $D^0$ (blue and red triangles), estimated as in (\ref{impTE}) using T-Learner BART. Right panel: group-specific CATE estimates (blue and red dashed lines) obtained from the two regressions in (\ref{lastXstep}), and final weighted CATE estimates (green dashed line) obtained from the re-balancing step in (\ref{Xweight}).}
		\label{XBART}
	\end{figure}

	Figure \ref{XBART} offers a simple example of a X-Learner, with BART as base-learner, applied to the one-covariate simulated data encountered in Figure \ref{example} and Figure \ref{STBART}. X-Learner's first step essentially derives, via T-Learner, the same BART estimates seen in the right panel plot of Figure \ref{STBART}. The output of the second step, namely the imputed treatment effects $\tilde{D}^1_i$ and $\tilde{D}^0_i$, are depicted in the left panel plot of Figure \ref{XBART} (red and blue triangles), together with the true ITE (grey dots). The graph on the right instead shows the estimated group-specific CATE $\hat{\tau}_1(\bm{x}_i)$ (blue dashed line) and $\hat{\tau}_0(\bm{x}_i)$ (red dashed line), derived from the two regressions in (\ref{lastXstep}), and the final CATE estimate $\widehat{\tau}(\bm{x}_i)$ (green dashed line), obtained from the weighting step in (\ref{Xweight}). Propensity score estimates employed for the weighting step were retrieved via probit version of BART. Notice that the final CATE estimates $\widehat{\tau}(\bm{x}_i)$ lie in between the two fitted group-specific $\hat{\tau}_1(\bm{x}_i)$ and $\hat{\tau}_0(\bm{x}_i)$.

	\subsubsection{R-Learners}

	R-Learner was proposed by \cite{wager_2019} as a two-stage meta-algorithm, which aims at minimizing a loss function, specifically defined on CATE, through parameter tuning. The derivation of the two-step procedure stems from \cite{robinson_1988} decomposition of the outcome model in (\ref{onestage}). Define the following two quantities:
	\begin{align}  \label{robinson}
	\begin{split}
	Y_i & =  \mu(\bm{X}_i) + \tau(\bm{X}_i) Z_i + \varepsilon_i   \\
	m(\bm{X}_i) & = \mathbb{E} (Y_i \mid \bm{X}_i) = \mu(\bm{X}_i) + \tau(\bm{X}_i) \pi(\bm{X}_i)
	\end{split}
	\end{align}
	as the outcome model and the \emph{conditional mean outcome} model, respectively. Under this setup, unconfoudedness assumption implies that the error term is such that $\mathbb{E} \big[ \varepsilon_i \mid \bm{X}_i, Z_i \big] = 0 $. Notice that under this parametrization $\tau(\bm{X}_i)$ (CATE) enters explicitly in the outcome regression model. By combining the two quantities above, \cite{robinson_1988} noticed that the outcome model can be re-written as:
	\begin{equation} \label{transout}
	Y_i - m(\bm{X}_i) = \Big(Z_i - \pi(\bm{X}_i) \Big) \tau (\bm{X}_i) + \varepsilon_i ~ .
	\end{equation}
	Starting from this decomposition, \cite{wager_2019} derive a loss function that can be used for parameter tuning in the estimation of CATE; the optimal CATE estimates are defined as the minimizer of the following loss function:
	\begin{equation} \label{tauargmin}
	\tau(\bm{X}_i) = \argmin_{\tau} \Biggl\{ \mathbb{E} \Bigg[ \Big( \big(  Y_i - m(\bm{X}_i)  \big)   -  \big(  Z_i - \pi(\bm{X}_i)  \big)  \tau(\bm{X}_i)  \Big)^2 \Bigg]  \Biggr\} ~ .
	\end{equation}
	The idea is that a base-learner that relies heavily on parameter tuning (e.g.~random forest or gradient boosted trees) can be tuned on the modified parametrization of the outcome model in (\ref{transout}), which includes a version of the outcome net of the impact of the covariates $m(\bm{X}_i)$ and propensity score balancing, instead of being tuned on the raw outcome $Y_i$ as one would do in an S- or T-Learner procedure.
	Since we cannot observe directly the quantities in (\ref{tauargmin}) for the minimization problem, the R-Learner overcomes the issue through the following two-step approach, i.e., by:
	\begin{enumerate}
		\item Fitting $\widehat{m}(\bm{x}_i)$ and $\widehat{\pi} (\bm{x}_i)$, by minimizing usual prediction error.

		\item Plugging in estimates from the first step to estimate $\hat{\tau} (\bm{x}_i)$, by minimizing the regularized sample equivalent of (\ref{tauargmin}), via cross-validation parameters tuning, that is:
		\begin{equation}   \label{rlearner}
		\begin{split}
		\widehat{\tau}(\bm{X}_i) & = \argmin_{\tau} \Big \{  \widehat{L}_n \big( \widehat{\tau}(\bm{X}_i) \big)  + \Lambda_n \big(  \widehat{\tau}(\bm{X}_i) \big)  \Big \} ~ , \quad \text{where}  \\
		\widehat{L}_n \big( \widehat{\tau}(\bm{X}_i) \big) & = \frac{1}{N} \sum_{i = 1}^{N}  \Bigg(  \Big(  Y_i - \widehat{m}^{-i} (\bm{X}_i)  \Big)  -  \Big( Z_i - \widehat{\pi}^{-i} (\bm{X}_i)  \Big)  \widehat{\tau}(\bm{X}_i)  \Bigg)^2 ~ ,
		\end{split}
		\end{equation}

	\end{enumerate}

	\noindent and where $\Lambda_n \big(  \hat{\tau}(\cdot) \big)$ is a term representing regularization (e.g.~Ridge or LASSO penalization, splines smoothness penalization, etc.), and the super-script $(-i)$ refers to the $i$-th observation being held-out from the estimation subsample, and used for leave-one out cross validation (LOO-CV).

	The advantage of R-Learners lies in the fact that they take a two-stage procedure where the first step produces an estimate of propensity score $\widehat{\pi} (\bm{x}_i)$, and by doing so mitigates confounding generated by the non-randomized selection mechanism. While in the second step, exploiting the parametrization in (\ref{robinson}), they make use of a loss function through which parameters can be optimized to estimate $\widehat{\tau}(\bm{x}_i)$ directly. R-Learners naturally work better with types of base-learners that rely heavily on parameter tuning, such as ensemble methods (e.g.~gradient boosting, random forests, etc.).

	\subsubsection{Multitask-Learners}

	The idea of multitask-learning, or multi-output learning, for causal inference was introduced by \cite{vanderschaar_2017} and \cite{vanderschaar_2018}, in the context of Gaussian Processes. The multitask perspective on CATE estimation consists in viewing the two potential outcomes $Y_i^{(0)}$ and $Y_i^{(1)}$ as output of a function $f: \mathcal{X} \rightarrow \mathbb{R}^2$, with $d$-dimensional input space and $2$-dimensional output space, where the output space is indexed by $Z_i$, that acts as ``task identifier". CATE estimator is defined as the difference between the elements of the $2$-dimensional output of $f (\cdot)$, i.e.,
	\begin{equation} \label{multitau}
	\widehat{\tau}(\bm{x}_i) = \widehat{f}_1(\bm{x}_i) - \widehat{f}_0(\bm{x}_i) = \mathbf{\hat{f}^{\top}} (\bm{x}_i) \bm{\xi} ~ , \quad ~ \text{where} \quad \bm{\xi^{\top}} = [ -1 ~~ 1 ] ~ .
	\end{equation}

	Equation (\ref{multitau}) displays a very similar formulation to a T-Learner; and as in the T-Learner procedure, the sample is divided into the two subgroups for the estimation. However, the advantage of viewing CATE estimation as a multitask problem is that, instead of estimating the two potential outcomes independently as one would do in a T-Learner or X-Learner, they are estimated ``jointly", through the specification of some hyperparameters that trigger a joint optimization for the two ``tasks": learning $f_{Z=1}$ and $f_{Z=0}$. Hence, this approach fits separate conditional mean functions (as in a T-Learner), but at the same time attempts to recover common underlying patterns between the two groups (as in an S-Learner) that would be otherwise lost due to the sample split.

	In the case of \cite{vanderschaar_2017}, multitask learning is induced through the specification of a particular structure in the stationary kernel function of the Gaussian Process prior. The method is labelled as \emph{Causal Multitask Gaussian Process} (CMGP). \cite{vanderschaar_2018} then proposed a similar method where multitask learning is induced via a non-stationary version of the GP kernel function (\emph{Non-Stationary Gaussian Process} - NSGP). Both causal multitask GPs allow to compute posterior credible interval for CATE as natural result of their Bayesian approach.

	\cite{vanderschaar_2017} and \cite{vanderschaar_2018} are two example of Multitask-Learners employing Gaussian Process regression as a base-learner, but there are different ways of inducing multitask learning using other types of base-learners (e.g.~tree ensembles). Due to their similarity with T-Learners in deriving a CATE estimator, we expect Multitask-Learners to perform well when complexity of the response surfaces $f_1, f_0$ varies across groups, and CATE itself turns out to be a rather complex function.

	\subsubsection{$\tau$-Learners}  \label{sec:taulearner}

	The last type of Meta-Learner presented in this work was developed by \cite{hahn_2020}, under the name of ``Bayesian Causal Forest". The authors tackle the problem of CATE estimation with a Bayesian approach, where they exploit the same parametrization seen in the context of R-Learners. Particularly, they noticed that the parametrization
	\begin{equation} \label{bcf}
	Y_i = \mu(\bm{X}_i) + \tau (\bm{X}_i) Z_i + \varepsilon_i ~ ,
	\end{equation}
	can be viewed as a Bayesian regression framework where the \emph{prognostic score}, defined as the impact of the covariates $\bm{X}_i \in \mathcal{X}$ on the outcome $Y_i$ in absence of the treatment, $\mu(\bm{x}_i) = \mathbb{E} \big[ Y_i \mid Z_i = 0, \bm{X}_i = \bm{x}_i \big]$, plays the role of the intercept, while $\tau (\bm{x}_i)$ the role of the slope.
	In this perspective, CATE is an explicit parameter of the model and thus can be treated in a Bayesian fashion through the specification of a prior distribution $p \big( \tau(\cdot) \big)$, which can be shaped to convey prior knowledge and more targeted regularization that can capture even simple patterns of heterogeneity. Bayesian Causal Forest of \cite{hahn_2020} is composed by a pair of separate and independent BART priors placed on $\mu (\cdot)$ and $\tau (\cdot)$ respectively, but the parametrization in (\ref{bcf}) can be exploited using other Bayesian regression methods (e.g.~Gaussian Process, Dirichlet Process regression, etc.).

	In addition to the parametrization shown in (\ref{bcf}), \cite{hahn_2020} make use of the two-stage procedure seen in (\ref{twostage}), Section \ref{sec:ProblemSetup}. The two-stage approach is motivated by the presence of a particular type of confounding, which the authors in \cite{hahn_2018} and \cite{hahn_2020} call \emph{Regularization Induced Confounding} (RIC). The intuition behind RIC is the following: regularization applied directly on the two curves $f_1, f_0$ may have unintended consequences on the induced regularization on $\tau (\cdot)$, leading to biased estimates of CATE. RIC is shown to have a stronger effect when there is strong confounding, such as in presence of \emph{targeted selection}, that is when individuals are selected into treatment based on the prediction of otherwise adverse outcome. Targeted selection implies a potential strictly monotone relationship between the propensity score $\pi (\bm{x}_i)$ and the prognostic score $\mu(\bm{x}_i) = \mathbb{E} \big[ Y_i \mid \bm{X}_i, Z_i = 0 \big]$, and is rather common in studies of observational nature. The proposed way to tackle confounding from targeted selection is precisely to use the two-stage representation illustrated in (\ref{twostage}), where a probabilistic estimate of the propensity score $\widehat{\pi} (\bm{x}_i)$, obtained from the first stage regression, is added to the covariates for the estimation of $\mu(\bm{x}_i) = \mathbb{E} \big[ Y_i \mid \bm{X}_i, Z_i = 0 \big]$ in the second stage, to account for their potential relationship.

	We name the above approach $\tau$-Learner, as it involves an explicit parametrization in terms of $\tau (\bm{x}_i)$ and a direct Bayesian approach to CATE estimation. \cite{hahn_2020} specifically make use of BART for estimation of $\mu(\bm{x}_i)$ and $\tau(\bm{x}_i)$, but any other Bayesian method could potentially work. As a further advantage, the direct Bayesian approach allows to obtain valid posterior credible intervals for CATE.

	\subsection{Model Selection} \label{sec:modelsel}

	Model selection is a challenging problem in causal inference. The main reason is that researchers cannot access counterfactual outcomes $Y_i^{(1-Z_i)}$ and thus observe the difference $Y_i^{(1)} - Y_i^{(0)}$, $\forall i \in \{ 1,..., N \}$, which distinguishes it from other classical model selection problems. The aim here is  to select a model $\mathcal{M}^* \in \{ \mathcal{M}_1, ..., \mathcal{M}_d \}$, which minimizes a loss function of the estimated CATE $\hat{\tau}$. We will review some of the methods that have been proposed, partly following the review in \cite{tibshirani_2018}.

	Standard likelihood-based criteria such as AIC and BIC are doomed to fail in a structural missing data context such as the one depicted. Other popular metrics in the ML literature might be used, such as prediction error in estimating the two conditional average potential outcomes $\mu_{z_i} (\bm{x}_i) = \mathbb{E} [ Y_i \mid \bm{X}_i = \bm{x}_i, Z_i = z_i ]$, that are the expression of the two surfaces $f_1 (\cdot)$ and $f_o (\cdot)$, that is
	\begin{equation}  \label{prederr}
	\widehat{\mathbb{E}}_{\mu} [\mathbb{\ell}(\widehat{\mu}_z, \mu_z)] = \frac{1}{N} \sum_{i = 1}^{N} \Big( \widehat{\mu}_{z_i} (\bm{x}_i) - y_i^{(z_i)} \Big)^2 ~ .
	\end{equation}
	A slightly modified version of (\ref{prederr}) takes into account the observational nature of the data at hand and uses inverse propensity score estimates for weighting the prediction error in estimating each $\mu_{z}$, that is,
	\begin{equation} \label{muIPTW}
	\widehat{\mathbb{E}}_{\mu,IPTW} [\mathbb{\ell}(\widehat{\mu}_z, \mu_z)] = \frac{1}{N} \sum_{i = 1}^{N} \frac{\big( \widehat{\mu}_{z_i} (\bm{x}_i) - y_i^{(z_i)} \big)^2}{\widehat{\pi}_{z_i} (\bm{x}_i)} ~ .
	\end{equation}

	The main reason why these metrics work poorly is that they do not express a direct loss function on CATE. CATE squared-loss function, of the type $\mathbb{\ell}( \hat{\tau} , \tau ) = \mathbb{E} [(\hat{\tau} - \tau)^2] $, was defined as \textit{Precision in Estimating Heterogeneous Treatment Effects} (PEHE) by \cite{hill_2011} and takes the following form:
	\begin{equation}  \label{PEHE}
	\mathbb{E} \big[ (\hat{\tau}_i - \tau_i)^2 \mid \bm{X}_i = \bm{x}_i \big] ~ .
	\end{equation}
	PEHE would be the ideal loss function to use, but cannot be computed since one of the two potential outcomes is unobservable. Attempts have been made in the literature to render the estimation of PEHE feasible, through plug-in estimates of $\tau_i$, but understandably none of them has been successful and commonly used so far. These approaches generally propose to use a Meta-Learner to estimate CATE on a portion of training sample $N_{train}$, then to derive CATE estimates also on a portion of held-out validation sample $N_{val}$ (where $N=N_{train} + N_{val}$). The training set CATE estimates, denoted by $\tilde{\tau} (\bm{x}_i)$, will act as the ground-truth CATE, while the validation set estimates, denoted by $\hat{\tau} (\bm{x}_i)$, are instead used to validate the model. Given these two quantities, an estimator of PEHE is constructed as
	\begin{equation} \label{taurisk}
	\widehat{\mathbb{E}}_{\tau} [\mathbb{\ell}(\hat{\tau}, \tau)] = \frac{1}{N_{val}} \sum_{i = 1}^{N_{val}} \big( \hat{\tau} (\bm{x}_i) - \tilde{\tau} (\bm{x}_i) \big)^2 ~ .
	\end{equation}
	A slightly different version of (\ref{taurisk}) replaces the training set estimate of CATE $\tilde{\tau} (\bm{x}_i)$ with a version of the observed outcome $y_i$ weighted by the inverse propensity score, to correct for selection bias typically found in observational studies, and reads:
	\begin{equation} \label{tauIPTW}
	\widehat{\mathbb{E}}_{\tau, IPTW} [\mathbb{\ell}(\hat{\tau}, \tau)] = \frac{1}{N_{val}} \sum_{i = 1}^{N_{val}} \left( \hat{\tau} (\bm{x}_i) - \frac{ (2z_i - 1) y_i}{\tilde{\pi}_{z_i}(\bm{x}_i)}  \right)^2 ~ ,
	\end{equation}
	where $\tilde{\pi}_{z_i}(\bm{x}_i)$ is a propensity score estimate obtained on the training set $N_{train}$. The advantage of this version is that instead of employing the training set CATE estimate $\tilde{\tau} (\bm{x}_i) $ as ground-truth, it replaces it with $\frac{ (2z_i - 1) y_i}{\tilde{\pi}_{z_i}(\bm{x}_i)} $, which contains observed $y_i$ and $z_i$, and an estimate of the propensity score that can be obtained through supervised approaches.

	However, the evident setback of these plug-in methods is that they exhibit high model-dependent bias, which derives from the fact that the plug-in estimates of the ground-truth CATE $\tilde{\tau} (\bm{x}_i) $ (or those of the propensity score) obtained on a training subset of the data, clearly depend on the method used for the estimation. In this way $\widehat{\mathbb{E}}_{\tau} [\mathbb{\ell}(\hat{\tau}, \tau)]$ actually measures the error of $\widehat{\tau} (\bm{x}_i)$ in estimating the training set CATE $\widetilde{\tau} (\bm{x}_i)$, rather than in estimating the ground-truth CATE $\tau(\bm{x}_i)$.

	Finally, it is worth noticing that the R-Learner presented in the previous section develops an actual loss function for $\tau(\cdot)$ in (\ref{rlearner}). However this can be exploited exclusively in the context of R-Learners (e.g.~comparing a random forest R-Learner with a BART R-Learner), as it implies assuming \cite{robinson_1988} parametrization as in (\ref{robinson}).

	\begin{table}
		\caption{Tested models}
		\label{testedmodels}
		\small
		\centering
		\begin{tabular}{lcl}
			\toprule

			\textbf{Family}     &  \textbf{Name}     &  \textbf{Description}  \\
			\midrule

			\multirow{3}{3cm}{Linear Models}  &  \multirow{1}{2.7cm}{\centering S-OLS}  &   Linear regression as S-Learner  \\
			&  T-OLS  &   Linear regression as T-Learner \\
			&  R-LASSO  &   LASSO regression as R-Learner \\
			\midrule

			Naive Non-Parametrics	&  $k$NN  &   $k$-Nearest Neighbors \\

			\midrule

			\multirow{5}{3.2cm}{Tree-Based Methods}  &  \multirow{1}{2cm}{\centering CF}  &   Causal Forest \\
			&  S-BART  &   BART as S-Learner \\
			&  T-BART  &   BART as T-Learner \\
			&  BCF  &   Bayesian Causal Forest \\
			&  R-BOOST  &   Gradient Boosting as R-Learner \\

			\midrule

			\multirow{2}{3.2cm}{Gaussian Processes}  &  \multirow{1}{2cm}{\centering CMGP}  &   Causal Multi-task Gaussian Process \\
			&  NSGP  &   Non-Stationary Gaussian Process \\

			\bottomrule
		\end{tabular}
	\end{table}

	\section{Simulation studies} \label{sec:simu}

	In this section we report and comment on results from two different semi-simulated studies, carried out to compare performance of some of the models presented above in estimating CATE. A third supplemental semi-simulated study can be found in the Appendix Section \ref{sec:appendix}. A semi-simulated study consists in simulating only the outcome surface $Y_i$ in the tuple $\mathcal{D}_i = \{ \bm{X}_i, Z_i, Y_i \} $, starting from real-world $\bm{X}_i$ and $Z_i$. In the case of observational semi-simulations, \cite{hill_2011} introduced a practical way of recreating an observational study from a randomized one. This is essentially done by leaving out a non-random portion of the treated group, so that treatment assignment is no longer random. Recreating an observational study from a purely randomized one has the main advantage of allowing control over the selection mechanism and that the common support assumption holds for the treated group.

	We provide results based on the analysis of two real world randomized controlled trials, after transforming them into observational studies. The first semi-simulated setup employs the IHDP dataset, firstly introduced by \cite{hill_2011} and popular in both the computer science and statistics literature on CATE estimation. The second and third setups instead employ the ACTG-175 dataset, available in the \textsf{R} package \textsf{speff2trial}, and differ in the way the outcome and ITE are generated. Both datasets are publicly available. As mentioned above, we present here below results from the IHDP data simulation and one of the two setups involving the ACTG-175 data, while we leave the other ACTG-175 setup in the Appendix Section \ref{sec:appendix}.

	The models tested include: linear regressions with S-Learner and T-Learner approach respectively (S-OLS and T-OLS); a penalized LASSO regression, through which we induce variable selection, with R-Learner approach (R-LASSO); a simple non parametric $k$-Nearest Neighbors ($k$NN) as T-Learner; a Causal Forest (CF), proposed by \cite{athey_2019} as a modified version of Random Forests, where after partitioning the sample space on the covariates $X$, the last split in each tree in the ensemble is made on $Z$; Bayesian Additive Regression Trees as S-Learner and T-Learner respectively (S-BART and T-BART); a Bayesian Causal Forest (BCF), i.e.,~a $\tau$-Learner; Gradient Boosted trees with R-Learner procedure (R-BOOST); and finally the two causal multitask versions of Gaussian Processes, one with stationary kernel (CMGP) and the other with non-stationary kernel (NSGP), developed by \cite{vanderschaar_2017} and \cite{vanderschaar_2018} respectively. A summary of the tested models is provided in Table \ref{testedmodels} below.

	For each of the two datasets analyzed, in order to provide a comparison of the methods presented above, we computed $\sqrt{\mbox{PEHE}}$ estimates for each of the $B=1000$ Monte Carlo simulations, and we averaged it over all the simulations. Estimates of $\mbox{PEHE}$ were obtained through its sample equivalent, namely:
	\begin{equation}
	\mbox{P}\widehat{\mbox{EH}}\mbox{E}_{\tau} = \frac{1}{N} \sum_{i=1}^{N} \Big(  \tau(\bm{x}_i) - \widehat{\tau} (\bm{x}_i) \Big)^2 ~ ,
	\end{equation}
	where $\widehat{\tau} (\bm{x}_i)$ is the CATE estimate obtained under the given method, while $\tau(\bm{x}_i)$ is the ground-truth CATE, which is always unknown in the real world, but in this case is simulated. Data are randomly split in 70\% train set used to train the models, and 30\% test set to evaluate the model on unseen data. $\sqrt{\mbox{PEHE}}$ is reported for both train and test data.

	\DIFaddend \begin{figure}[t]
		\centering
		\begin{minipage}{.5\textwidth}
			\centering
			\includegraphics[width=1\linewidth]{Joy_Train.pdf}
		\end{minipage}%
		\begin{minipage}{.5\textwidth}
			\centering
			\includegraphics[width=1\linewidth]{Joy_Test.pdf}
		\end{minipage}
		\caption{\small $\sqrt{\mbox{PEHE}}$ distribution on train set (left) and test set (right), IHDP data.}
		\label{fig:IHDPdistr}
	\end{figure}

	\subsection{IHDP data}

	The first semi-simulated setup makes use of the IHDP dataset, popular in the literature for CATE estimation and  used for the first time in \cite{hill_2011}. It includes data taken from the Infant Health and Development Program (IHDP), a randomized controlled trial carried out in 1985, aimed at improving health and cognitive status of premature infants with low weight at birth, through follow-ups and parent support groups. The dataset includes 25 covariates, 6 continuous and 19 binary. The data are transformed into observational by leaving out a non-random portion of the treated individuals, namely those with non-white mothers. This leaves 139 observations in the treated group and 608 in the control group, for a total of 747 observations.

	ITE is derived as the difference between the simulated potential outcomes, which are generated as:
	\begin{equation} \label{hillPO}
	\begin{split}
	Y^{(0)} & \sim \mathcal{N} \Big( \exp \big( (X+W) \beta_B \big), 1 \Big) ~ ,  \\
	Y^{(1)} & \sim \mathcal{N} \big( X \beta_B - \omega^b_B, 1 \big)  ~ ,
	\end{split}
	\end{equation}
	where $W$ is an offset matrix of same dimension as $X$ with every entry equal to 0.5, and $\beta_B$ is a 25-dimensional vector of regression coefficients $\big( 0, 0.1, 0.2, 0.3, 0.4 \big)$, sampled with probabilities $ \big( 0.6, 0.1, 0.1, 0.1, 0.1 \big)$ (``Response surface B'' in \cite{hill_2011}).
	For each simulation $b \in \{ 1, ..., 1000 \}$, $\omega^b_B$ is an offset chosen to guarantee that the population-level ``Average Treatment effect on the Treated'' (ATT) is equal to $ATT = \mathbb{E} \big[ Y^{(1)} - Y^{(0)} \mid Z = 1 \big] = 4$ .

	\begin{table}
		\caption{\small IHDP and ACTG-175 data. $\sqrt{\mbox{PEHE}_{\tau}}$ estimates $\pm$ 95\% confidence intervals for each tested model, on the train and test sets respectively.}
		\label{simuresults}
		\footnotesize
		\centering
		\begin{tabular}{ccc c cc}
			\toprule

			&   \multicolumn{2}{c}{\textbf{IHDP}} & & \multicolumn{2}{c}{\textbf{ACTG-175}}  \\

			\cmidrule(r){2-3}  \cmidrule(r){5-6}

			& Train  &  Test &   & Train  &  Test  \\

			\midrule

			S-OLS  		&	4.78 $\pm$ 0.32		&  4.78  $\pm$ 0.32  &   &  0.299 $\pm$ 0.003	&   0.300 $\pm$ 0.003  \\
			T-OLS  		&	1.94 $\pm$ 0.13	&  2.02 $\pm$ 0.13      & 	&  0.591 $\pm$ 0.007	&  0.601 $\pm$ 0.008  \\
			R-LASSO  		&	1.99 $\pm$ 0.15 	&   2.05 $\pm$ 0.16      &	 &  0.309 $\pm$ 0.006 	&    0.310 $\pm$ 0.006	\\

			\midrule

			$k$NN 		&	2.67 $\pm$ 0.19 	&   2.98 $\pm$ 0.21      &   &  0.903 $\pm$ 0.004	&    0.902 $\pm$ 0.005	\\

			\midrule

			CF  		&	3.71 $\pm$ 0.25	&  3.68 $\pm$ 0.25      &    &  0.292 $\pm$ 0.004	&   0.291 $\pm$ 0.004	\\
			S-BART  		&	1.93 $\pm$ 0.13 	&  2.06 $\pm$ 0.14      &	&  0.298 $\pm$ 0.003 	&  	0.299 $\pm$ 0.003	\\
			T-BART  		&	0.94 $\pm$ 0.04 	&  	1.41 $\pm$  0.09  &	 &  0.530 $\pm$ 0.005 	&  	0.530 $\pm$ 0.006	\\
			BCF  		&	0.88 $\pm$ 0.05		&     1.34 $\pm$ 0.10  &   &  \textbf{0.281 $\pm$ 0.003}	&  \textbf{0.282 $\pm$ 0.003}	\\
			R-BOOST  		&	2.28 $\pm$ 0.14 	&  	2.48 $\pm$ 0.16     &	&  0.408 $\pm$ 0.006 	&   0.405 $\pm$ 0.006	\\

			\midrule

			CMGP  		&   0.64 $\pm$ 0.04 	&	0.76 $\pm$ 0.06      &	   &  0.354 $\pm$ 0.009 	&	0.346 $\pm$ 0.008	\\
			NSGP		&  	\textbf{0.55 $\pm$ 0.02} 	&	\textbf{0.71 $\pm$ 0.05}      &  	&  0.293 $\pm$ 0.003 	&	0.292 $\pm$ 0.003	\\

			\bottomrule
		\end{tabular}
	\end{table}

	Simulation results are reported in Table \ref{simuresults} and Figure \ref{fig:IHDPdistr}. The best models appear to be the multitask Gaussian Processes (CMGP and NSGP) of \cite{vanderschaar_2017} and \cite{vanderschaar_2018}. This does not come as unexpected for two reasons related to the way the potential outcomes are simulated. First, the simulated setting tends to accommodate Meta-Learners fitting separate surfaces $f_1, f_0$ (T-Learners, X-Learners, Multitask-Learners), since the two simulated potential outcomes display very different degrees of complexity. Indeed, T-Learners in Table \ref{simuresults} show better performance than their S-Learner counterparts (e.g.~T-BART vs S-BART). The second reason is that potential outcomes generated in (\ref{hillPO}) are rather smooth, and Gaussian Process regression tends to adapt well to non-linear continuous surfaces. We report also, in Figure \ref{fig:IHDPdistr}, the distribution of $\sqrt{\mbox{PEHE}_{\tau}}$ over the $B=1000$ iterations on both train and test data, for each of the models. We also notice that tree-based methods are relatively more prone to overfitting.

	An important remark about the data generating process described by (\ref{hillPO}) is that it does not really induce strong confounding, since it is easy for a non-parametric base-learner to distinguish the two underlying polynomials $\mathbb{E} [ Y^{(Z_i)} \mid \bm{X}_i = \bm{x}_i ]$. Besides, the fact that noise around the conditional mean is independently simulated for the two potential outcomes produces extra noise around CATE that often results in unrealistically large variation in the simulated ITE. In the ACTG-175 simulated example illustrated in Section 4.2, we will follow the parametrization found in \cite{wager_2019} and \cite{hahn_2020} in the data generating process of the outcome surface, in order to induce stronger confounding (which is believed to be common in observational studies) and avoid creating unnecessarily high noise around CATE.

	\subsection{ACTG-175 data}

	The second semi-simulated setup presented here is created using ACTG 175 clinical trial data, freely available in the \textsc{R} package \textsf{speff2trial}. The data come from a randomized placebo-controlled trial aimed at comparing monotherapy versus a combination of therapies in HIV-1-infected subjects with CD4 cell counts between 200 and 500 (\cite{doi:10.1056/NEJM199610103351501} for details). As in the case of IHDP data, an observational study is recreated by throwing away a non-random subset of patients, namely those not showing symptomatic HIV infection. The final dataset consists of 813 observations and 12 variables (3 continuous and 9 binary). We simulate continuous outcome and treatment effects by trying to resemble sign and magnitude of the estimated treatment effect fitted with a simple linear model, and also the scale of the original outcome (CD4 cell counts, on a log scale). Variables included in the dataset are shown in Table \ref{AIDSvar}.

	\begin{table}
		\caption{ACTG 175 data variables}
		\centering \footnotesize
		\begin{tabularx}{10.5cm}{l | X}

			\textbf{Variable}  &  \textbf{Description} \\
			\midrule

			\textit{age} &  Numeric	\\
			\textit{wtkg} &	Numeric	\\
			\textit{hemo} &	 Binary (hemophilia = 1)	\\
			\textit{homo} &	 Binary (homosexual = 1)	\\
			\textit{drugs} & Binary (intravenous drug use = 1)	\\
			\textit{oprior} & Binary (non-zidovudine antiretroviral therapy prior to initiation of study treatment = 1) \\
			\textit{z30} &	Binary (zidovudine use in the 30 days prior to treatment initiation = 1)	\\
			\textit{preanti} &	Numeric (number of days of previously received antiretroviral therapy)	\\
			\textit{race} &	 Binary	\\
			\textit{gender} &	Binary 	\\
			\textit{str2} &	 Binary: antiretroviral history (0 = naive, 1 = experienced)   \\
			\textit{karnof\_hi} &	Binary: Karnofsky score (0 = $<100$, 1 = $100$)
			\label{AIDSvar}
		\end{tabularx}
	\end{table}

	Unlike the case of IHDP data, response surface $Y_i$ is not generated via simulation of the two potential outcomes. Instead, we generate continuous outcome $Y_i$ according to the parametrization
	\begin{equation}
	Y_i = \mu(\bm{X}_i) + \tau(\bm{X}_i) Z_i + \varepsilon_i ~ ,
	\end{equation}
	which allows to specify CATE directly, instead of starting from the simulation of potential outcomes, and to induce more confounding through the prognostic score function $\mu(\bm{x}_i)$. The prognostic score $\mu(\bm{x}_i)$ and CATE $\tau(\bm{x}_i)$ are generated as:
	\begin{equation}  \label{setup1}
	\begin{split}
	\mu(\bm{x}_i) = & ~ 8 - 0.07hemo - 0.002 | wtkg - 1 | + 0.06gender - \frac{0.1}{age + 2} \\
	& + 0.007karnof_{hi} - 0.1z30  - 0.05race
	\\[8pt]
	\tau(\bm{x}_i) = & ~ 0.1 + 0.1age \cdot (karnof_{hi} + 2) ~ .
	\end{split}
	\end{equation}
	Then noise is added by simulating normally distributed errors $\varepsilon_i \sim \mathcal{N} (0, \sigma^2)$, with standard deviation equal to $\sigma =   (\mu_{max} - \mu_{min}) / 2$, where $\mu_{max}$ is the maximum value of the generated prognostic score across the sample, that is $\mu_{max} = \big \{ \mu(x_j) : \mu(x_j) > \mu(\bm{x}_i), ~ \forall i \in \{ 1,...,N \} \backslash \{j\} \big \}$, while $\mu_{min}$ is the minimum. Notice that, unlike the case of IHDP data, the error term is not simulated independently for the two groups, which avoids imposing too much noise around CATE. This translates into generally lower and less variable estimated $\mbox{PEHE}$ for all models compared to the case of IHDP data, as shown both in Table \ref{simuresults} and Figure \ref{fig:AIDSdistr}.

	In this simulated setting CATE is of rather simple and known form. Hence, we expect S-Learner to outperform their T-Learner counterparts. Indeed, this is exactly what Table \ref{simuresults} shows. S-Learners all display lower $\sqrt{\mbox{PEHE}}$ than T- and Multitask-Learners, with the only exception of NSGP, which instead performs quite well. BCF turns out to be the best method in this scenario, thanks to its ability to adjust to stronger confounding and apply more regularization when learning CATE, which allows to capture small and weakly heterogeneous treatment effect patterns. Figure \ref{fig:AIDSdistr} depicts again the distribution of $\sqrt{\mbox{PEHE}}$ over the $B=1000$ simulations, for both train and test data, for all the tested models.

	We employ ACTG-175 data also in a third semi-simulated setup featuring more complex $\mu(\bm{x}_i)$ and $\tau(\bm{x}_i)$ surfaces, compared to the ones in (\ref{setup1}). Description and results of this third example are provided in the Appendix Section \ref{sec:appendix}.

	\begin{figure}[t]
	\centering
	\begin{minipage}{.5\textwidth}
		\centering
		\includegraphics[width=1\linewidth]{AIDS_Joy_Train.pdf}
	\end{minipage}%
	\begin{minipage}{.5\textwidth}
		\centering
		\includegraphics[width=1\linewidth]{AIDS_Joy_Test.pdf}
	\end{minipage}
	\caption{\small $\sqrt{\mbox{PEHE}}$ distribution on train set (left) and test set (right), ACTG-175 data.}
	\label{fig:AIDSdistr}
\end{figure}

	\section{The effect of school meal programs on health indicators} \label{sec:realdata}
	\DIFaddend

	Social, economic and political sciences are among the disciplines that rarely benefit from RCTs studies, and often have to rely on observational or non-fully randomized data to address causal questions. In this section we provide an example of a real world social sciences application of the models reviewed in previous sections, to demonstrate their usefulness in applied research. To this end, we make use of the \textsf{nhanes\textunderscore bmi} dataset found in the \textsf{ATE} \textsf{R} package. The dataset consists in a subset of the 2007-2008 National Health and Nutrition Examination Survey (NHANES), containing the variables shown in Table \ref{NHANES}. The data are employed in this context to investigate whether participation to school meal programs is associated to improved children's health indicators. There is reason to believe that the impact of such programs is heterogeneous across individuals, in that demographics such as age, gender or ethnicity are likely to play a role in how effective participation is (e.g.~younger kids might benefit more than older ones, etc.). This supports the use of methods for individual treatment effects estimation.

	Contrary to a simulated study, in a real world application it is impossible to recognize the type of heterogeneity patterns one encounters. For this reason, it is usually good practice to deploy more than just a single model to estimate CATE. In this specific context, we derive and compare CATE estimates from both Bayesian Causal Forest (BCF) and non-stationary Multitask Gaussian Process (NSGP), as from the simulation studies they emerged to be the most flexible methods when it comes to adapt to different types of surfaces.

	We find that CATE estimates obtained from the two models are highly correlated, with Pearson's and Spearman's coefficients equal to 0.93 and 0.92 respectively, but BCF exhibits higher variance in the estimates (sample standard deviation of 0.09) compared to NSGP (sample standard deviation of 0.03). This means that BCF estimates capture much more heterogeneity. The left panel plot of Figure \ref{figNHAMES} depicts the densities of CATE estimates derived from each model, and NSGP's estimates appear to be far more concentrated around the mean and generally indicate larger impact. The main source of the difference in dispersion between the two obtained estimates is likely to be child's age covariate. As illustrated in the right panel plot in Figure \ref{figNHAMES}, BCF attributes higher degree of treatment heterogeneity to child's age compared to NSGP, although they both indicate slightly non-linear monotone relationship between the two quantities.

	The impossibility of directly measuring performance of the methods prevents from carrying out proper model selection, as outlined in details in Section \ref{sec:modelsel}. However, it appears that both models support the existence of treatment heterogeneity, although with different patterns. In particular, BCF detects a stronger moderating effect of child's age on treatment effect than NSGP, while to other covariates (such as gender and ethnicity) the two methods attribute very similar moderating effects.

	\begin{table}
	\caption{\small NHANES data variables}
	\centering \footnotesize
	\begin{tabularx}{9.7cm}{l | X}

		\textbf{Variable}  &  \textbf{Description} \\
		\midrule

		\textit{BMI} &  Numeric. Body mass index (outcome variable)	\\
		\textit{school\textunderscore meal} &	Binary (treatment indicator)	\\
		\textit{age} &	 Numeric (child's age)	\\
		\textit{childSex} &	 Binary (male = 1)	\\
		\textit{afam} & Binary (African American = 1)	\\
		\textit{hisam} & Binary (Hispanic = 1) \\
		\textit{povlev\textunderscore 200} & Binary (family above 200\% federal poverty lvl = 1)	\\
		\textit{sup\textunderscore nutr} &	Binary (supplementary nutrition program = 1)	\\
		\textit{stamp\textunderscore prog} & Binary (food stamp program = 1)	\\
		\textit{food\textunderscore sec} &	Binary (food security in household = 1)	\\
		\textit{ins} &	 Binary (any insurance = 1)   \\
		\textit{refsex} &	Binary (adult respondent gender is male = 1)  \\
		\textit{refage} &	Numeric (adult respondent's age)

		\label{NHANES}
	\end{tabularx}
\end{table}

	\section{Conclusion} \label{sec:discussion}
	\DIFaddend

	In the previous sections, we covered the most recent developments on the estimation of heterogeneous treatment effect in the context of non-randomized studies. We use here the phrase ``non-randomized" studies as the methods illustrated are well suited for observational studies, but also for RCTs that might display threats to randomization such as imperfect compliance and non-random missingness (obviously also pure RCTs, sample size permitting).

	Our review and simulation studies lead to a few general observations. A first simple but practical suggestion in applied setting where CATE is to be estimated, is to implement more than just one single method, to provide robustness to the results, given that it is impossible to discern which particular scenario the researcher faces (weak vs strong heterogeneity, simple vs complex CATE patterns, etc.). Secondly, it is clear from the semi-simulated exercises that non-parametric regression methods adapt well to non-linearities, and in general to complexity of the CATE function, whereas parametric ones are naturally less robust to functional form misspecification, unless heterogeneity patterns are very weak (reverting to homogeneity). In addition, Bayesian approaches to non-parametric modelling sacrifice some computational cost for less cumbersome parameter tuning and nicer posterior properties (i.e., posterior credible interval).

	A final remark concerns $\tau$-learners, and BCF in particular. Bayesian Causal Forest of \cite{hahn_2020} appears to show nice adaptability and robustness to misspecification regarding the CATE surface (i.e., smoothness, sparsity, heterogeneity vs homogeneity of treatement effect) compared to other Meta-Learners examined. This ability stems directly from the explicit inclusion of $\tau(\bm{x}_i)$ as a ``parameter" in a regression framework, which is designed specifically to detect moderating impact of the covariates on treatment effects, and allows to place a direct prior distribution in a Bayesian fashion. From an empirical perspective, this implies some degree of flexibility in the choice of $p \big( \tau(\bm{x}_i) \big)$. In fact, the prior on $\tau(\bm{x}_i)$ can translate into a more or less targeted regularization based on prior knowledge on the complexity of CATE surface. In a case where CATE is believed to be complex one might, for example, modify default options to have deeper trees in the estimation of $\tau$. On the contrary, in case little or no a priori knowledge is possessed, a more agnostic prior, which prevents overfitting via stronger regularization and leaves some degree of freedom in CATE estimation, can be chosen. This is something extraneous to the other methods found in the literature that, although the promising performances (the ``causal'' Gaussian Processes, CMGP and NSGP, in particular), do not allow to treat CATE as a model parameter and to specify of a prior distribution on it, thus not allowing to convey a priori knowledge, which is frequently of interest in causal analysis.

\DIFaddend \begin{figure}[t]
	\centering
	\begin{minipage}{.5\textwidth}
		\centering
		\includegraphics[width=1\linewidth]{DENS_TAU.pdf}
	\end{minipage}%
	\begin{minipage}{.5\textwidth}
		\centering
		\includegraphics[width=1\linewidth]{AGE_TAU.pdf}
	\end{minipage}
	\caption{\small Left panel: kernel density of CATE estimates from the two models (dashed lines are the sample means). Right panel: CATE estimates from the two models as a function of child's age.}
	\label{figNHAMES}
\end{figure}

	\nocite{*}
	\bibliographystyle{chicago}

redpajama_set_name: RedPajamaArXiv


text:
Microsoft to Enter Formula One?
Posted in F1 Design, F1 History, FIA rules, Formula 1, Mclaren, Microsoft, Motor Racing on July 13th, 2006
One of my commenters has pointed out a BBC report that suggests the FIA has agreed that Microsoft supply all engine control units in F1. The full thread of the resulting conversation can be read on my post, Here's to the FIA.
What can I say?
Apart from all the MS jokes that will inevitably result, the news raises some interesting questions. Further digging has revealed that Microsoft are in partnership with a McLaren subsidiary called McLaren Electronic Systems. If true, this would surely give McLaren a potentially unfair advantage over the other teams.
Given Microsoft's complete lack of a track record in ECU systems, one must presume that the real work would be done by McLaren's subsidiary. And how tempting will it be for them to ensure that the software works best with Mercedes engines? In fact, that would be quite natural, as the engine most readily available for testing the system would be their own. Is this why Fernando Alonso is so sure that McLaren are going to be the team to beat in the future?
It has been reported today that the World Motor Sport Council has approved the changes to the 2006 and 2007 F1 regulations. The concentration has been on the reduction of final qualifying from 20 minutes to 15, and the elimination of third cars from Friday practice. Presumably any mention of standardized engine control units will wait until the engine homologation issue has been signed, sealed and delivered.
The whole business reminds me of the FIA's introduction of standardized fuel delivery systems several years ago. The company awarded the contract for these produced a system that has given constant trouble with sticking fuel nozzles, yet the teams were forced to soldier on with them. Will we see cars stopping inexplicably in the middle of a race while monitor screens suddenly turn blue and advise that the engine has committed an illegal operation? Will the teams have to phone MS to activate the software every time they change an engine? The mind boggles.
If the report turns out to be accurate, I can see an interesting struggle developing between the FIA and the constructors. With McLaren Electronic Systems involved, it might well become a war resulting in a split similar to CART's. Can the FIA really be serious?
5 Responses to "Microsoft to Enter Formula One?"
[...] Update on Microsoft to Enter Formula One?: [...]
By Who Owns Formula One? » Formula 1 Motor Racing News and Views on July 20th, 2006 at 1:22 pm
[...] The potential storm over the FIA's choice of manufacturer for a standardized engine control unit, as mentioned in my post Microsoft to Enter Formula One?, is brewing nicely. [...]
By Microsoft, McLaren and the FIA » Formula 1 Motor Racing News and Views on July 24th, 2006 at 1:49 pm
I may not be old enough to remember 'the good old days' but it seems to me that the FIA is doing whatever they can to make this sport more like A1GP only with more powerful cars and more egotistical drivers. I understand that they are trying to keep costs down for the sake of everyone and i also believe that some changes are necessary to make the races more competitive. I would like to see a series in which the engines are regulated but the areo rules are very lax. Give me your Hammerhead noses, fan cars, active suspensions, wheel inlets, maybe even different engine displacements based on overall car weight. Massive amounts of money will always be thrown at problems in F1, instead of a 1% increase gained from money spent, why not have that millions go to a crazy new design that will revolutionize the sport.
By Dan M on March 14th, 2007 at 1:51 pm
That's a very interesting idea, Dan. When we consider the innovations springing from F1 over the last forty years or so, it seems a pity that such inventiveness should have the lid put on it through regulations designed to keep costs down. In the early days of GP racing (before the FIA was even a gleam in anyone's eye), there was no accepted norm for engines and some pretty weird solutions were tried – so, if we free up the aerodynamics, we could see some amazing things resulting!
By Clive on March 14th, 2007 at 1:57 pm
[...] that the FIA has been able to ban traction control in F1 from 2008, thanks to the introduction of the standardized Microsoft/McLaren ECU. Anything that removes driver aids from the cars and puts more emphasis on the skill of the driver [...]
By Formula 1 Latest » Decisions From the FIA on March 31st, 2007 at 1:51 pm
Art NYC
Beauty on a Budget
The Money Blog

redpajama_set_name: RedPajamaCommonCrawl


text:
Manhattan season 2 episode 11 gravity falls full free
Watch manhattan online full episodes all seasons yidio. With sean astin, jason biggs, greg cipes, rob paulsen. Manhattan season two recaps atomic heritage foundation. The second season focuses on the challenges of designing and assembling the gadget for the trinity test, the security situation and espionage, and personal relationships.
Take back the falls s2 e20 online watch online anytime. When donnie upgrades metalheads artificial intelligence, he must decide if he can trust a sentient metalhead. Manhattan returns to wgn america on tuesday october for its second season and it hasnt made one of those leaps that might cause critics to yell into the void that this is a classic drama. If you like manhattan 20142015 you might like granite flats 202015 or taken 2002. So, i travel a lot, and get bored alot, so i was wondering if their was anywhere i could download full seasons of gravity falls for like a phone. Gravity falls reaction series season 2 episode 11 youtube. The gadget dangles from its tower in the pouring rain. Watch gravity falls season 2 episode 11 online not what he seems.
After agents powers and trigger take stan into custody, dipper and mabel begin to wonder how much they actually know about their grunkle stan. In the opening scene, we flash forward to july 16, 1945, the morning of the trinity test. Gravity falls show summary gravity falls is just west of weird, slightly east of eerie and always north of normal in a fresh new comedy series that follows city kids and twins dipper and mabel pines who are sent to spend the summer with their eccentric great uncle stan aka grunkle stan who runs a tourist trap called the mystery shack in a. Scott called it an absolute delight funny and charming, fastmoving and full of surprises. Simmons joins the cast as of season 2, episode 12, a tale of two stans, and plays. Aug 03, 2014 watch full gravity falls season 2 episode 2. This will always be the last on anybodys list of gravity falls. Meanwhile lil gideon takes advantage of dippers and mabels height discrepancy. A new era is ushered in with the first test of an atomic weapon when the series concludes. Northwest mansion mystery is counted 11 instead of 10. Gravity falls season 2 episode 11 watch online the full episode. The 50 best things to watch on disney plus right now the new. Gravity falls s02e11 not what he seems video dailymotion. It takes place during the manhattan project during the year 1943.
Here is a list of our recaps of each episode of season two of wgn americas television show manhattan. He is known for his roles as kevin girardi in the television series joan of arcadia, ethan haas in the class, sean walker in the nbc series the event, dipper pines in gravity falls, and. Watch gravity falls episodes online season 2 2019 tv guide. Simmons joins the cast as of season 2, episode 12, a tale of two stans, and plays stans twin brother, stanford fillbrick ford pines. Good diva, bad diva promo today, total divas season 4. A new era is ushered in with the first test of an atomic weapon. The first season of the american animated television series gravity falls consisted of 20 episodes on the disney channel, and aired from june 15, 2012 to august 2, 20. Just east of eerie gravity falls is an american animated television series produced by disney television animation for disney channel. Good diva, bad diva online free megavideo, total divas season 4. Oct 21, 2015 here is a list of our recaps of each episode of season two of wgn americas television show manhattan.
It takes place during the time when countries had their scientists racing to be the first to build an atomic bomb. What was that about disney it was a full 30 minute episode, what up with that. Watch manhattan season 1 episode 1 online tv fanatic. A development overseas reinvigorates the hunt for a spy on the hill. The series was created by alex hirsch, a former writer for the marvelous misadventures of flapjack and fish hooks. Jul 27, 2014 july 2 1943 61 countries at war, over 40 million casualties. July 2 1943 61 countries at war, over 40 million casualties.
With rachel brosnahan, michael chernus, christopher denham, alexia fast. Disneys emmy awardwinning animated series gravity falls reaches its final showdown in the highly anticipated onehour finale episode on monday, february 15 7. Good diva, bad diva sneak peek, total divas season 4. And wow, the tenth episode of the season and shit is going. In this fantasy series, a young woman is drawn to a small maine town and. Here is the annalysis for gravity falls season 2 episode 11 subscribe and like the video. Click here and start watching the full episode in seconds. Oct 09, 2015 manhattan returns to wgn america on tuesday october for its second season and it hasnt made one of those leaps that might cause critics to yell into the void that this is a classic drama. Gravity falls season 2 episode 20 weirdmageddon 3 take back the falls duration. The second and final season of the american animated television series gravity falls began on. After agents powers and trigger take stan into custody, dipper and mabel begin to wonder how.
Rotten tomatoes reports a 92% approval rating, with an average rating of 8. The only reason why i didnt give this episode a full 10 is because of the opening. Tvguide has every full episode so you can stayuptodate and watch your favorite show gravity falls anytime, anywhere. Feb 22, 2015 the only reason why i didnt give this episode a full 10 is because of the opening. Gravity falls tv show season 2 episodes list next episode. Good diva, bad diva full episode long, total divas season 4. Twins dipper and mabel pines are sent to spend the summer with their. May 05, 2020 a guide listing the titles and air dates for episodes of the tv series gravity falls. Mayor night from the story gravity falls season 3 by cartoonfuel hannah j. Oct 15, 2015 season two of manhattan begins with a bright, mysterious light. Season 2 may 25, 2018 the final season for the prematurely cancelled manhattan is an allaround strong period piece, albeit one that sometimes mistakes heavy. A man in is the desert playing golf during a terrible storm, looking at his notes and listening to the news.
Disney xd to air finale of emmywinning series gravity falls. Season two of manhattan begins with a bright, mysterious light. Season 2 guide for gravity falls tv series see the episodes list with schedule and episode summary. Watch gravity falls season 2 full episodes watch cartoon online. Welcome to the new series of gravity falls season 2. In the processes of investigating gravity falls, stans brother met up with bill the illuminati incarnate that helped him design a machine that would give a person the knowledge of the universe, but bill had an agenda to destroy the world with the information overload. Toby determined, a journalist for gravity falls gossiper, is voiced by gregg turkington, and will forte plays tyler cutebiker.
Not what he seems is the 11th episode of the second season of gravity falls and the 31st. Running with scissors season 3, episode 11 in which the bigcarrie affair hits rock bottom. Teenage mutant ninja turtles metalhead rewired tv episode. One of the hallmarks of the entire series was mabels constant. A guide listing the titles and air dates for episodes of the tv series gravity falls. The episode also serves as the midseason finale of season 2.
Watch gravity falls season 1 episode 11 online little dipper. Dipper is selfconscious that mabel is taller than him and seeks a magic solution. The free pizza poster from headhunters is seen on the agents record of stan. Get your team aligned with all the tools you need on one secure, reliable video platform. Season two renewal for wgn america series by trevor kimball, october 14, 2014 though its viewership has been quite small, manahattan. Filmed in new mexico by lionsgate television, skydance television and tribune studios. Feb 25, 2015 here is the annalysis for gravity falls season 2 episode 11 subscribe and like the video. Manhattan is about making the first atomic bomb, like the movie fat man and little boy 1989. Gravity falls s02 e11 not what he seems video dailymotion. Anyway, tuned in march 9th for the eleventh episode of season two, not what he seems. Gravity falls season 2 subtitles english 67 subtitles. Jason ritter born february 17, 1980 is an american actor.
The series was created by alex hirsch, a former writer for the marvelous misadventures of. Nathan fillion, stana katic, susan sullivan, molly quinn, ruben santiagohudson, tamala jones, jon huertas, seamus dever. Now that ive warmed to your characters and highminded story about the creators of the atom bomb toiling in the new mexico desert amid discussions of good and evil and physics and adultery, please pick up the pace. Bigs wife, natasha, catches carrie in their apartment, chases her out, and on the way falls. In the season 2 premiere, charlie is forced into a new leadership role in the race to build the worlds first. Watch manhattan season 1 episode 1 online and lets take a trip to new mexico. Gravity falls season 2 episode 2 video dailymotion. Buy, rent gravity falls, season 2 episode 20, is available to watch and stream on disney. Manhattan season 2 episode 10 manhattan season 2 episode 10 manhattan season 2 episode 10 manhattan season 2 episode 10 manhattan season 2 episode 10. Gravity falls prepares for the end by focusing on dipper and mabels bond.
1370 33 374 593 25 521 723 409 1537 557 712 1348 169 1286 1460 920 611 1367 1565 352 1122 680 795 719 787 463 1385 374 767 1034 1025 1066 492 222 1434 1492 429 626 5 433 1091


text:
Adams, Abraham (1753)
Adams, Henry (1531)
Adams, Jackie Devon (1947)
Adams, Jeremy (1604)
Adams, Jerry Hudnall (1898)
Adams, John (1500)
Adams Sr, John (1666)
Adams, Mary Ann (1791)
Adams, Rose (1536)
Airs / Ayers, John (1721)
Akin, Isabella America (1863)
Alexander, Andrew (Drew) (1860)
Alexander, Joseph C. (1891)
Alexander, Martha Melvina (1866)
Alexander, Thomas Houston (1809)
Alexander, William Lee (Pa) (1900)
Allen, Rainey Edgar (1919)
Anderson, Sarah ' Sally ' Iola (1886)
Andrews, Catherine M. (1652)
Antonucci, Faustino (1739)
Antonucci, Fillippo (1819)
Antonucci, Giacomo Antonio 'John George' (1863)
Antonucci, Pasquale (1795)
Antonucci / George, Elena 'Lee' Francesca (1924)
Antonucci / George, Florence 'Flo' (1919)
Antonucci / George, John 'Jack' (1915)
Antonucci / George, Juan 'Jean' Marie (1914)
Antonucci / George, Kathryn 'Kay' J. (1914)
Antonucci / George, Mary Jean (1907)
Antonucci / George, Susan 'Suzie' (1906)
Armstrong, Annie Mary (1907)
Armstrong, Etta Belle (1901)
Armstrong, George Lee (1873)
Armstrong, Ida Ramoth (1873)
Armstrong, Jefferson Benjamin (1846)
Armstrong, Lula Green (1876)
Arnold, Annie Elizabeth (1881)
Arnold, Clara 'Rudy' (1899)
Arnold, Elsie (1915)
Arnold, James Albert (1858)
Arnold, John (1807)
Arnold, John Clyde (1905)
Arnold, Lois (1890)
Arnold Sr., Luther Lee (1894)
Arnold Jr., Luther Lee (1927)
Arnold, Margaret (1851)
Arnold, Martha Mae (1921)
Arnold, Robert 'Ruben' (1886)
Arnold, Wilma Irene (1916)
Aten, Adrian (1695)
Aten, Hendrick 'Henry' M. (1723)
Aten, John (1765)
Atkinson, Hezakiah (1853)
Aultman, James Franklin (1945)
Aultman, Woodrow Wilson (1912)
Austin, William ' Buck' Lewis (1921)
Auten, Alexander Colley (1839)
Auten, Almira Matilda (1852)
Auten, Annie Ethel (1873)
Auten, Emma Lovina (1856)
Auten, George Russell (1903)
Auten, George W. (1875)
Auten, Henry (1789)
Auten, Henry Clinton (1839)
Auten, Jacob Cotner (1813)
Auten, Jacob Wesley (1844)
Auten, James Finley (1859)
Auten, John Ellison (1840)
Auten, Josephine Elizabeth (1908)
Auten, Mary Eliza (1840)
Auten, Mary Elizabeth (1885)
Auten, Myrtle Dell (1870)
Auten, Ralph Wesley (1906)
Auten, Rebecca Alice (1854)
Auten, Robert Love (1879)
Auten, William Allen (1846)
Aycock, Amy (1897)
Aycock, Leonard (1905)
Aycock, Lula Bama (1893)
Aycock, Simon (1863)
Aycock, Zeffie C.B. (1898)
Ayers, Abigail Saphronia (1820)
Ayers, Allie Mae (1902)
Ayers, Asa William (1846)
Ayers, Benjamin Franklin (1842)
Ayers, Capt. John 'Jack' David (1791)
Ayers, Elbert L. (1856)
Ayers, Elizabeth 'Lizzie' Frances (1897)
Ayers, Emily (1878)
Ayers, Ishmael A. (1850)
Ayers, James Thomas (1868)
Ayers, James William (1871)
Ayers, Jeremiah (1883)
Ayers, Jesse McKinley (1896)
Ayers, John (1763)
Ayers, John David (1876)
Ayers, John Jefferson (1907)
Ayers, John Wood (1933)
Ayers, Louraine (1870)
Ayers, Martha Charity (1875)
Ayers, Martha Jane (1870)
Ayers, Penellope Sarah (1821)
Ayers, Queen Ann (1849)
Ayers, Rebecca Jane (1872)
Ayers, Samuel (1855)
Ayers, Telitha 'Annie' Ann Frances (1824)
Ayers, Thomas W. (1841)
Ayers, Wealthy Elizabeth (1844)
Ayers, Wilburn Asa (1878)
Ayers, William (1798)
Ayers, William Asa (1895)
Ayers, William Henry (1872)
Ayers, William Nathaniel (1842)
James Franklin Aultman
Woodrow Wilson Aultman
11/10/1912 - 1/6/2005
Mary Lucille Davis
About James Franklin Aultman
James Franklin Aultman was born on August 24, 1945 in Calhoun County, Florida and died on May 29, 2014 in Calhoun County, Florida. He lived to be 68 years old.
Called "Jimmy" . He married Cynthia L.
Mr. James Franklin (Jimmy) Aultman, age 68, of Blountstown, FL passed away early Thursday morning, May, 29, 2014 in Panama City, FL.
Jimmy was born on August 24, 1945 in Marianna, FL and had lived in Blountstown all of his life. Jimmy retired with the Florida Department of Corrections with 24 years of service. He was a member of Dixie Masonic Lodge #109 F&AM in Blountstown, FL, a member of the York Rite Bodies in Tallahassee, FL, a member of the Shaddai Shrine Temple in Panama City, FL. Jimmy enjoyed doing carpentry and electrical work and was an avid fisherman. He was a member of the Florida Peace Officers Association.
He was an excellent pistol shooter and a member of the Single Action Pistol Society, where he participated in shooting cowboy action matches. Jimmy was a member of the First Baptist Church in Blountstown, FL.
Survivors include, Wife, Cynthia Aultman of Blountstown, FL Daughter, Cindy Golden and her husband, Grant of Blountstown, FL Son, Derek Aultman of Live Oak, FL Sister, Margaret (Aultman) Davis of Easley, SC
7 grandchildren, Brit Golden, Lane Golden, Kaylin Bontrager, Kelsey Bontrager, Alex Aultman, Kaydan Aultman and Alicia Aultman. 3 great-grandchildren, Mason Wiehl, Konnor Parish and Karter Parish
Funeral services will be held Sunday, June 1, 2014 at 3:00 pm (CDT) at the First Baptist Church in Blountstown with Reverend Tim Rhoads officiating. Interment will follow in Nettle Ridge Cemetery in Blountstown with Masonic Rites. The family will receive friends Sunday, June 1, 2014 from 2:00 pm (CDT) until service time at 3:00 pm (CDT) at the First Baptist Church. All arrangements are under the direction of Marlon Peavy at Peavy Funeral Home in Blountstown, FL.
Death Certificate 5/29/2014 Blountstown FL (Calhoun County) 68 yrs old
Name:	James Franklin Aultman, Event Place:	Blountstown, Calhoun, Florida, Birth Date:	Aug.24 , 1945, Death Date:	May 29 , 2014, Cemetery:	Nettle Ridge Cemetery
Click on any of the News, Events, or Discoveries buttons above to see historical things that happened during James Aultman's life. These are only some of the major events that affected the life and times of James, his family, and friends. For example, James is 9 years old when The nationally televised McCarthy hearings lead to Senator Joseph McCarthy's downfall; Segregation by race in schools is declared unconstitutional by the Supreme Court; Congress adds the words "under God" to the Pledge of Allegiance and requires "In God We Trust" to appear on all American currency.
1 1946 The university of Pennsylvania develops the ENIAC computer, containing 18,000 vacuum tubes
2 1947 Researches at Bell Labs invent the first transistor
6 1951 The first color television is introduced in the U.S.
7 1952 The first sex-change surgery was performed to change George Jorgensen into Christine Jorgensen.
65 2010 CERN scientist trap anti-matter for the first time.
5 1950 NEWS HEADLINES: On June 25, the Korean War begins and the United Nations officially declaring war on North Korea two days later.
66 2011 March 15, 2011 - Largest recorded earthquake hits Japan and everyone watches the horrifying scenes as survivors post videos taken with their cell phones on YouTube.
9 1954 The nationally televised McCarthy hearings lead to Senator Joseph McCarthy's downfall; Segregation by race in schools is declared unconstitutional by the Supreme Court; Congress adds the words "under God" to the Pledge of Allegiance and requires "In God We Trust" to appear on all American currency.
65 2010 April 20, 2010 - A British Petroleum deep water oil rig explodes in the Gulf of Mexico, causing the largest oil spill in the history of the United States, killing eleven workers, and devastating the envirnoment. It also severely damaged the fishing and tourism industries of gulf states.
66 2011 March 15, 2011 - The largest recorded earthquake hits Japan followed by a huge tsunami. Three nucular power plants are danaged and remain in critical status for the next six weeks. Manufacturing plants in Japan can no longer produce and impacts the world economy.
James Aultman
The Children of James Franklin Aultman


text:
Seven Hills Charter Public School Feasibility Study
Hope High School Auditorium Study
Moses Brown Squash Courts
Lincoln School STEAM Hub for Girls
Community Preparatory School Expansion
Lincoln School Master Plan
Head of School Residence
Salvation Army Child Care Center
Moses Brown School Facilities Audit
Community Preparatory School Feasibility Studies
LLB Architects is well-acquainted with the Lincoln School campus, having completed a limited feasibility study for the School in 2005. In 2016, LLB led a master planning process during which we asked: What does a leading school for girls look like? For Lincoln, it is a place "where tradition meets innovation." The school's location is one of its best advantages but also one of its greatest challenges.
Our team reviewed Lincoln School's Strategic Plan and previous Master Plans, existing building plans, spaces, and infrastructure. We conducted a comprehensive facility conditions assessment to provide order-of-magnitude cost analysis on the state of Lincoln's buildings and grounds, combined with quantitative data extracted from newly created 3D models. We held multiple meetings with the steering committee, stakeholders, faculty departments, Buildings and Grounds, and the steering committee to identify and itemize program requirements and issues. In addition, we engaged many Lincoln School constituents through the use of online surveying to rapidly identify shared opinions of Lincoln School's strengths and weaknesses.
The findings presented in this planning study were a culmination of the team's multifaceted research, interviews, analysis, creative thinking, and consensus-building.


More and more people are choosing blinds over traditional curtains for more stylish and modern windows. With hundreds of blinds available, our partners even offer next day delivery to Somerset.	We have a wide selection of window coverings available to accomodate almost any need from traditional wooden blinds through to fancy new velux blinds. Off the peg or made to measure!
If you are buying made to measure blinds in Somerset click on one of the images above. And if you don't see the type of product you need, then just click on the closest to find your way to the blinds that are exactly right for you in Somerset.
If you look after your blinds, they should last a long long time. Just remember to follow the manuacturer's guidelines and you will get the maximum life. Our partners can offer most of the window solutions available on the market, and deliver them direct to your door in Somerset.
This site is for blinds in Somerset. We offer a wide variety of wooden blinds, vertical blinds, roller blinds, venetian blinds, blackout blinds, panel blinds, velux blinds, roman blinds, and all can be made to measure for Somerset.


November 29, 2017 in Movies, News.
November 28, 2017 in Movies, News.
November 24, 2017 in Movies, News.
November 22, 2017 in Movies, News.
November 21, 2017 in News, Television.
November 20, 2017 in Movies, News.
November 18, 2017 in Movies, News.
November 17, 2017 in Movies, News.
November 16, 2017 in Movies, News.
November 15, 2017 in Movies, News.