The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
93,349 characters
Median-based Splitting Rules for Causal Trees and Forests
\newtheorem{assump}{Assumption}
\maketitle
\begin{abstract}
Heavy-tailed and skewed outcomes are common in the randomized experiments and observational studies used to estimate heterogeneous treatment effects, yet the mean-squared-error criterion that guides splitting in honest causal trees is sensitive to the extreme values they generate.
Building on the causal forest framework \citep{athey_recursive_2016, wager_estimation_2018}, we introduce the Median Squared Deviation (MSD) criterion, which replaces the leafwise difference in means in the honest splitting objective with the Hodges--Lehmann location estimator while leaving honest leaf estimation and forest inference unchanged. Two further median-based rules, the Median Absolute Deviation (MAD) and the Least Median of Squares (LMS), serve as robust baselines. We evaluate the criteria in a simulation study covering precision, bias, and confidence interval coverage. MSD restricts its robustness to split selection and lowers the error of conditional average treatment effect estimates under heavy-tailed and skewed outcomes. Further, we re-visit two empirical applications: the first analyzes the electoral effects of a Mexican conditional cash transfer program on precinct-level observations, while the second application studies antiretroviral treatments in HIV-positive adults.
\end{abstract}
\begin{keywords}
Heterogeneous treatment effects, causal forests, Hodges--Lehmann estimator, robust estimation, heavy-tailed outcomes
\end{keywords}
\begin{jelclass}
C14, C21, C52
\end{jelclass}
\section{Introduction}
\label{sec:Introduction}
Extreme and heavy-tailed outcomes pose a recurring challenge for estimating treatment effects in randomized trials and observational studies in economic applications \citep{li_robust_2021, li_robust_2023, Athey2023, ghosh_robustness_2026}. A natural and established remedy is to base inference on the ranks of the outcomes rather than their magnitudes \citep{lehmann_nonparametrics_1975} within the potential outcomes framework \citep{imbens_causal_2015}. Building on this idea, \cite{Rosenbaum1993} combines the location estimator of \cite{hodges_estimates_1963} with rank statistics to estimate a constant additive treatment effect. \cite{ghosh_robustness_2026} establish an asymptotic theory for this rank-based estimator.
Although the constant treatment effect assumption serves as a convenient and widely used starting point for causal inference \citep{rosenbaum_covariance_2002, Athey2023}, many applications in economics, social sciences, and public health exhibit treatment effect heterogeneity \citep{kunzel_meta-learners_2019, yao_survey_2021, kennedy_towards_2023, cinelli_challenges_2025}. Motivated by the goal of combining nonparametric heterogeneous treatment effect (HTE) estimation with robustness to heavy-tailed or skewed outcomes, this paper extends the causal forest framework \citep{athey_recursive_2016, wager_estimation_2018} by using the Hodges--Lehmann estimator as the basis for median-based splitting rules in the underlying causal trees. The construction retains the honest sample splitting of \citet{athey_recursive_2016} and \citet{wager_estimation_2018} and we use pointwise confidence intervals based on the bootstrap of little bags (BLB) \citep{sexton_standard_2009, athey_generalized_2019} to quantify uncertainty around the conditional average treatment effect (CATE) estimates.
Nonparametric estimation of HTEs has become a central theme in modern causal inference, motivated by the need to understand how treatment effects vary across subpopulations without imposing restrictive parametric assumptions on the outcome model. In the potential outcomes framework \citep{imbens_causal_2015}, where each unit has unobserved counterfactual responses under treatment and control, a broad class of flexible methods exists, from regression trees and forests to meta-learners and Bayesian nonparametric models \citep{wager_estimation_2018, kunzel_meta-learners_2019, hahn_bayesian_2020}, that adaptively model treatment effects as functions of covariates without specifying a fixed functional form.
These approaches are useful in settings where the true relationship between covariates, outcomes, and treatment effects is complex and high-dimensional. While much of the HTE literature focuses on smooth or well-behaved outcome distributions, heavy-tailed or irregular potential outcomes are prevalent in real-world applications.
\cite{li_robust_2021, li_robust_2023} analyze electronic health record data to estimate and compare HTEs between two alternative therapies for hypertension with a right-skewed outcome variable, flexibly adjusting for a set of demographic and clinical covariates.
\cite{Athey2023} study two applications using skewed and heavy-tailed house price data and medical expenditures for pneumonia patients as outcome variables. In this paper, we revisit the conditional cash transfer experiment of \cite{de_la_o_conditional_2013}, also studied by \cite{ghosh_robustness_2026}. This randomized experiment evaluates the electoral impact of Mexico's conditional cash transfer program on skewed political support outcomes. Additionally, we analyze the ACTG 175 antiretroviral trial of \cite{hammer_trial_1996}, which was also used as an empirical application by \cite{leqi_median_2022}. For both applications, we examine potential treatment effect heterogeneity while accounting for their skewed outcomes.
Causal trees and their aggregation into causal forests provide a flexible, nonparametric framework for estimating the HTE by adapting recursive partitioning to target CATE rather than prediction error. \cite{athey_recursive_2016} and \cite{wager_estimation_2018} formalize an honest estimation approach by separating tree construction from leaf-level effect estimation to reduce adaptive bias and enable asymptotically valid inference, while forest aggregation and subsampling stabilize estimates in high-dimensional settings.
Subsequent work connects causal forests to semiparametric efficiency theory through orthogonalized and doubly robust score constructions \citep{chernozhukov_doubledebiased_2018, kennedy_towards_2023}, thereby mitigating bias from nuisance function estimation and improving robustness to model misspecification. Nonetheless, practical performance can deteriorate under heavy-tailed or skewed outcomes, where classical mean-squared-error splitting criteria become sensitive to extreme values \citep{galimberti_notes_2011}.
Within random forest algorithms based on the framework of \citet{Breiman2001}, splitting rules of the underlying trees typically rely on MSE minimization for regression tasks or impurity measures such as the Gini index for classification. For predictive tasks, several robust alternatives have been proposed.
\citet{Hastie2009} review absolute and Huber loss functions as splitting criteria to reduce sensitivity to extreme outcomes, and \citet{roy_robustness_2012} propose median-based strategies encompassing tree aggregation, splitting criterion choice, and outcome transformation to robustify random forest predictions. \citet{ishwaran2015effect} study weighted splitting criteria of which MSE and Gini index splitting emerge as special cases, while \citet{L1_splitting} compare log-rank tests and integrated absolute differences as splitting criteria in survival forests.
\citet{LiMartin2017} provide a unifying framework for forest-type regression that supports a wide range of robust loss functions, nesting classical random forests and quantile regression forests as special cases. Translating these robustness ideas to the causal forest setting is non-trivial, however, because the splitting criteria must be adapted to the honest sample-splitting framework of \citet{athey_recursive_2016} and \citet{wager_estimation_2018}. In this framework, the data used for split selection are separated from the data used for leafwise treatment-effect estimation. Beyond the standard MSE-based splitting rules for causal trees, \citet{athey_recursive_2016} discuss several alternative partitioning estimators, including transformed-outcome trees, fit-based trees, which choose splits by improvements in outcome fit, and t-statistic trees, which prioritize treatment-effect heterogeneity across candidate leaves. Further extensions include \citet{Lechner_Modified_CF}, who augment variance-based splitting by a propensity-score penalty, and \citet{chernozhukov_riesznet_2022}, who generalize the splitting mechanism via their ForestRiesz estimator.
Our approach differs in construction from the weighted least-absolute-deviation criterion of \citet{li_robust_2021}. We retain the honest, mean-based splitting objective of \citet{athey_recursive_2016} and \citet{wager_estimation_2018} and replace its difference-in-means anchor with the Hodges--Lehmann location estimator. To our knowledge, anchoring the honest splitting objective on the Hodges--Lehmann estimator has not previously been studied. Our primary contribution is the resulting Median Squared Deviation (MSD) criterion, anchored on this estimator within the causal forest framework. We implement a computationally convenient version of this criterion and evaluate its performance in simulations and applications.
We further examine two related heuristic variants, the Median Absolute Deviation (MAD) and the Least Median of Squares (LMS), which we treat as practical alternatives without the same formal guarantees.
The paper is structured as follows. Section~\ref{sec:PO_framework} reviews the potential outcome framework. Section~\ref{sec:causal_trees} introduces causal trees and the honest sample-splitting construction. Section~\ref{sec:MedianSplitting} introduces the Hodges--Lehmann estimator as a consistent estimator of the leafwise location shift and implements it in new splitting rules, primarily the Median Squared Deviation (MSD) criterion and, as heuristic variants, the Median Absolute Deviation (MAD) and the Least Median of Squares (LMS). Section~\ref{sec:simulations} presents a simulation study comparing the different causal tree splitting rules, embedded in the resulting causal forest estimators, in terms of ATE and CATE precision and confidence interval coverage. Section~\ref{sec:empappl} illustrates the methods on the Progresa conditional cash transfer experiment of \cite{de_la_o_conditional_2013} and the ACTG 175 antiretroviral trial of \cite{hammer_trial_1996}. Section~\ref{sec:conclusion} concludes.
\section{Potential Outcomes Framework}
\label{sec:PO_framework}
To formalize the HTEs motivated in Section~\ref{sec:Introduction}, we work within the potential outcomes framework, a standard approach for causal inference in economics and the social sciences, particularly for observational data and randomized experiments \citep{dominici_controlled_2021, hernan_miguel_a_what_nodate}.
The general idea was introduced by Jerzy Neyman in his 1923 work on randomized agricultural experiments \citep{splawa-neyman_application_1990}. It was developed further in a series of papers \citep{rubin_estimating_1974, rubin_inference_1976, rubin_bayesian_1978, rubin_randomization_1980} and addresses what \cite{holland_statistics_1986} later termed the fundamental problem of causal inference, that only one of the two potential outcomes is observed for each unit.
The resulting Rubin Causal Model \citep{rubin_rubin_2011, hernan_miguel_a_what_nodate} rests mainly on the three assumption of stable unit treatment values, unconfoundedness, and common support, as reviewed below.
We consider independent and identically distributed sampled data $(X_i, Y_i, D_i)$ for $i = 1, ..., N$ with the pair of potential outcomes $(Y_i(0), Y_i(1))$ for each unit $i$ based on a binary treatment indicator $D_i \in \{0,1\}$ where \(D_i=1\) denotes treatment and \(D_i=0\) denotes control, and a covariate vector $X_i$ of dimension $K$, where $K$ is the number of covariates. The individual treatment effect (ITE),
\begin{equation}
\label{eq:ITE}
\tau_i \coloneqq Y_i(1)-Y_i(0),
\end{equation}
is based on the existence of a potential outcome pair for each observation $i$. The notation of potential outcomes relies on the Stable Unit Treatment Value Assumption (SUTVA), which rules out interference between units and multiple or hidden versions of treatment. Assumption~\ref{assump:SUTVA} formalizes the link between the potential outcomes and the observed outcome.
\begin{assump}{Stable unit treatment value assumption (SUTVA).}
\label{assump:SUTVA}
\begin{align*}
\text{If } D_i = d \text{, then } Y_i(d) = Y_i^{obs} \text{ , } \forall d \in \{0, 1 \} \text{ , } \forall i \in \{1, ..., N \}.
\end{align*}
\end{assump}
Under Assumption~\ref{assump:SUTVA}, the observed outcome can be written as
\[Y_i^{obs}=Y_i(1)D_i+Y_i(0)(1-D_i).\]
Together with the no-interference and no-hidden-versions components of SUTVA, this implies that each unit has two well-defined potential outcomes and that unit $i$'s potential outcomes are unaffected by the treatment assignment of other units. Moreover, Assumption \ref{assump:SUTVA} removes different treatment manifestations from consideration, which is also known as treatment variation irrelevance. While Assumption \ref{assump:SUTVA} is a common but relatively strong assumption, modifications such as allowing for multiple versions of treatment \citep{vanderweele_ignorability_2008, vanderweele_causal_2013, laffers_identification_2020} or allowing for interference between individuals within certain groups of the population \citep{hudgens_toward_2008} have been proposed but are not the focus of this paper.
We further assume an unconfounded treatment assignment mechanism independent of potential outcomes in Assumption \ref{assump:unconfoundedness}. Given the covariates, the presence of unmeasured confounders is ruled out and leads to conditional independence between potential outcomes and treatment status.
\begin{assump}{Unconfoundedness.}
\label{assump:unconfoundedness}
\begin{align*}
D_i \perp\!\!\!\perp \left(Y_i(1), Y_i(0)\right) \mid X_i, \\
\text{or equivalently,} \quad
Pr(D_i \mid Y_i(1), Y_i(0), X_i) = Pr(D_i \mid X_i).
\end{align*}
\end{assump}
In perfectly randomized experiments, Assumption \ref{assump:unconfoundedness} is fulfilled by design as the assignment mechanisms are under the researcher's control. In the analysis of observational studies, this assumption and its sensitivity can, at best, be verified indirectly. Often, it remains untestable and has to be justified in the context of the research question at hand \citep{imbens_nonparametric_2004}.
Thirdly, we assume common support so that treatment effects are identifiable throughout the support of the covariates. Let the conditional treatment probability or propensity score be defined as \begin{equation}
\label{eq:propensity}
p(X_i=x) \coloneqq Pr(D_i=1|X_i=x).
\end{equation}
The following overlap condition in Assumption \ref{assump:overlap} requires at every covariate value in the support of \(X_i\), units have a nonzero probability of receiving either treatment state.
\begin{assump}{Overlap.}
\label{assump:overlap}
\begin{align*}
\exists\, \epsilon > 0 \text{ such that } \epsilon < p(X_i=x) < 1-\epsilon, \quad \forall x \text{ in the support of } X_i, \text{ with probability } 1.
\end{align*}
\end{assump}
Strong overlap improves the stability of treatment-effect estimates and is also
important for the construction of valid confidence intervals. Moreover, convergence rates of semi-parametric estimators depend heavily on the degree of common support. Two common solutions to ensure sufficient overlap are trimming based on the propensity score and the use of more robust inference procedures that estimate data adaptively based on the amount of covariate overlap \citep{damour_overlap_2021, busso_new_2014, rothe_robust_2017}. Assumptions \ref{assump:unconfoundedness} and \ref{assump:overlap} imply that treatment assignment is ignorable for the identification of treatment effects. Using Assumptions \ref{assump:SUTVA} to \ref{assump:overlap}, one can identify the conditional average treatment effect (CATE) as
\begin{align}
\begin{split}
\label{eq:CATE}
\tau(X_i = x) &\coloneqq \mathbb{E}[Y_i(1) - Y_i(0) \mid X_i = x] \\
&=
\underbrace{\mathbb{E}[Y_i \mid D_i = 1, X_i = x]}_{\coloneqq\mu_1(x)}
-
\underbrace{\mathbb{E}[Y_i \mid D_i = 0, X_i = x]}_{\coloneqq\mu_0(x)}.
\end{split}
\end{align}
The CATE in Equation~\eqref{eq:CATE} is the target estimand of this paper, and the splitting rules of Section~\ref{sec:MedianSplitting} are designed to create partitions of the covariate space that allow a robust estimation of the sample analog of \eqref{eq:CATE}. Note that the sample analogs of the conditional mean functions $\mu_1(x)$ and $\mu_0(x)$ can be estimated separately from observational data and that the ITE in Equation~\eqref{eq:ITE} should not be confused with the CATE in \eqref{eq:CATE}, since in general $\tau(X_i) \neq \tau_i$ \citep{post2024flexible, PostvandenHeuvel+2025}.
Averaging the CATE over the covariate distribution yields the Average Treatment Effect (ATE),
\begin{equation} \label{eq:ATE}
\tau = \mathbb{E}[Y_i(1) - Y_i(0)] = \mathbb{E}\left[\tau(X_i)\right],
\end{equation}
a scalar summary of the heterogeneity in \eqref{eq:CATE} that we also report in Sections~\ref{sec:simulations} and~\ref{sec:empappl}.
Because both the CATE and the ATE are functionals of the same conditional mean functions $\mu_1(x)$ and $\mu_0(x)$, a rich representation of the CATE in Equation~\eqref{eq:CATE} requires a flexible nonparametric estimator of these conditional means. Causal trees and forests provide such an estimator while supporting valid inference through honest sample splitting, as summarized in Section~\ref{sec:causal_trees}.
\section{Causal Trees and Causal Forests}
\label{sec:causal_trees}
By partitioning the feature space into rectangles, the general tree-based method builds upon the idea of fitting trivial separate models to each rectangle, aiming at estimating the expectation of the outcome variable conditional on the regressors \citep{Hastie2009}. These methods can cope very well with large-scale datasets while the structure of a single tree retains its interpretability. Throughout this paper we focus on regression trees based on the CART procedure proposed by \cite{cart84} and follow the notation of \cite{athey_recursive_2016}.\footnote{Related approaches we do not consider are, for instance, the Iterative Dichotomiser 3 (ID3) or multivariate adaptive regression splines (MARS) \citep{quinlan_induction_1986, friedman_multivariate_1991}.} A CART regression tree $\Pi$, also referred to as a partition, divides the feature space into separate segments in a recursive manner until reaching a set of leaves $\Pi = \{l_1, ..., l_L\}$ with a predefined minimum leaf size \citep{wager_estimation_2018}. Let $\mathcal{S}$ denote a sample with sample size $N$ and let $l(x; \Pi)$ be a specific leaf $l \in \Pi$ with $x \in l$. Then the leaf-level conditional mean function $\mu_d(x;\Pi)$ is the population average outcome within the leaf containing $x$,
\begin{align}
\label{eq:conditional_mean_function}
\mu_d(x; \Pi) = \mathbb{E}[Y_i|D_i=d, X_i \in l(x;\Pi) ].
\end{align}
Here $\mu_d(x; \Pi)$ is the leaf-averaged, partition-dependent counterpart of the pointwise population conditional mean $\mu_d(x)$ defined in \eqref{eq:CATE}, averaging $\mu_d(\cdot)$ over the leaf $l(x; \Pi)$ containing $x$ rather than conditioning on the single point $X_i = x$.
For a partition $\Pi$ chosen independently of a sample $\mathcal{S}$, the leaf sample mean is an unbiased estimate of this conditional mean function and is given by
\begin{equation}
\label{estim_cond_mean_func}
\hat\mu_d(x; \mathcal{S}, \Pi) \coloneqq \frac{1}{\# (i \in \mathcal{S}: D_i= d, X_i \in l(x; \Pi))} \sum_{i\in \mathcal{S}:D_i=d, X_i \in l(x; \Pi)} Y_i.
\end{equation}
Given a tree $\Pi$ and a test point $x$, the treatment effect within a specific leaf $l \in \Pi$ can thus be denoted as the difference of average outcomes between the treatment and control group in leaf $l$
\begin{equation}
\label{leaf_treatment_effect}
\hat\tau(x; \mathcal{S}, \Pi) \coloneqq \hat\mu_1(x; \mathcal{S}, \Pi) - \hat\mu_0(x; \mathcal{S}, \Pi).
\end{equation}
For simplicity and without loss of clarity, we abstract from explicitly including $\Pi$ or $\mathcal{S}$ in our further notation if the partition itself is not relevant for the main argument (i.e. $\hat\tau(x; \mathcal{S})$ instead of $\hat\tau(x; \mathcal{S}, \Pi)$).
We further use two interchangeable forms for any leaf-constant quantity: a pointwise form $\hat{\tau}(x; \mathcal{S}, \Pi)$ when predicting for individuals or summing over $i \in \mathcal{S}$, and a leaf-direct form $\hat{\tau}(l; \mathcal{S}, \Pi)$ when summing over $l \in \Pi$. Both refer to the same object via $l = l(x; \Pi)$, and the sample argument $\mathcal{S}$ is likewise suppressed when clear from context.
In the usual CART algorithm, recursive partitioning and leaf estimation use the same training data set. This adaptive procedure uses the same observations to choose the partition and to estimate the within-leaf quantities. An adaptive reuse of the data is prone to overfitting, because splits are chosen precisely to improve in-sample fit.
Tree complexity is typically controlled through cost-complexity pruning. After growing a large tree, a penalty on the number of terminal leaves is used to select a smaller subtree. The benefit of the conventional, adaptive CART procedure is its efficient data usage, as all training observations are considered for partitioning and estimation.
However, the CART algorithm tends to produce more extreme leaf means as spurious extreme values are grouped with high probability \citep{cart84, Hastie2009}.
\cite{athey_recursive_2016} address this problem by introducing honesty.
A tree is called honest if the data used to estimate within-leaf responses is disjoint from, and thereby independent of, the data used to determine the splits of the partition $\Pi$. Honesty is a property of the
estimation protocol, not of any specific algorithm: the honest causal tree of \cite{athey_recursive_2016}, the Double-Sample Tree of \cite{wager_estimation_2018}, and the Propensity Tree of \cite{wager_estimation_2018} all achieve it by slightly different constructions.
Honest leaf estimation mitigates overfitting bias that would otherwise arise from using the same data both to select $\Pi$ and to estimate $\hat\tau(l)$, and is what delivers per-leaf conditional unbiasedness for $\tau(l;\Pi)$.
Following \citet{athey_recursive_2016}, we distinguish between three disjoint samples by considering a training sample $\mathcal{S}^{tr}$ to build a tree, an estimation sample $\mathcal{S}^{est}$ to estimate effects for a given tree partition $\Pi$, and a test sample $\mathcal{S}^{te}$ to evaluate out-of-sample performance.
We denote subsample sizes by $N^{\star}_{d,l}$, where the superscript $\star \in \{\mathrm{tr}, \mathrm{est}, \mathrm{te}\}$ indicates the sample $\mathcal{S}^{\star}$, the subscript $d \in \{0,1\}$ indicates treatment status, and $l \in \Pi$ indicates the leaf. Subscripts are dropped when the corresponding restriction is not imposed, so that $N^{\mathrm{tr}}$ is the full training-sample size, $N^{\mathrm{tr}}_1$ the number of treated training observations, $N^{\mathrm{tr}}_l$ the number of training observations in leaf $l$, and $N^{\mathrm{tr}}_{1,l}$ the number of treated training observations in leaf $l$. The sample superscript is suppressed when clear from context.
Throughout the paper, we refer to a standard honest causal tree
\citep{athey_recursive_2016} that uses the following tree construction and estimation
procedure in four steps.
\begin{enumerate}[label=\textup{(}\roman*\textup{)}, ref=\roman*]
\item \label{itm:first} split the data into a training sample
$\mathcal{S}^{tr}$, an estimation sample $\mathcal{S}^{est}$, and a test
sample $\mathcal{S}^{te}$;
\item \label{itm:second} grow a tree $\Pi$ on $\mathcal{S}^{tr}$ by
maximizing the honest splitting criterion
$\hat Q^{\,\text{MSE, H}}_\tau$ in \eqref{eq:MSE_CT_H}, which trades the
in-sample heterogeneity reward
$\sum_l (N^{tr}_l/N^{tr})\,\hat\tau(l)^2$ against a leaf-variance
penalty weighted by $(1/N^{tr}_l + 1/N^{est}_l)$, and determine tree
complexity by $k$-fold cross-validation within $\mathcal{S}^{tr}$ via
cost-complexity pruning. We use $k$ here for the number of cross-validation folds, a standard convention that is unrelated to the observation indices elsewhere in the paper;
\item \label{itm:third} with $\Pi$ fixed, estimate leaf-level CATEs as the
difference-in-means on the held-out estimation sample,
$\hat{\tau}(l;\mathcal{S}^{est}) = \bar{Y}^{\,est}_{1,l} - \bar{Y}^{\,est}_{0,l}$. Because $\mathcal{S}^{est}$ and $\mathcal{S}^{tr}$ are independent, this estimator is conditionally unbiased for $\tau(l;\Pi)$, and
\item \label{itm:fourth} predict for a new $x$ by routing $x$ through
$\Pi$ to its leaf $l(x;\Pi)$ and reading
$\hat\tau(x;\Pi) = \hat\tau\bigl(l;\mathcal{S}^{est}\bigr)$.
\end{enumerate}
Step~(\ref{itm:second}) is the core component in the causal tree construction process in \citet{athey_recursive_2016} and is
where the four splitting rules in this paper differ. We introduce the
mean-based criterion of \citet{athey_recursive_2016} here in detail and use it
as baseline for constructing the median-based and heuristic rules in
Sections~\ref{sec:MSD} and \ref{sec:MAD_LMS}: each rule is presented through
an infeasible MSE on $(\mathcal{S}^{te},\mathcal{S}^{est})$ that defines
its population target, followed by an in-sample criterion on
$\mathcal{S}^{tr}$ that is the sample analog used for split selection.
Since the individual treatment effect $\tau_i = Y_i(1) - Y_i(0)$ is never
observed, the natural prediction risk
$\frac{1}{N^{te}}\sum_{i\in\mathcal{S}^{te}}(\tau_i - \hat\tau(X_i))^2$ is
infeasible. \citet{athey_recursive_2016} work with the modified MSE that
removes the unobservable component while preserving the ranking of candidate
trees, since $\tau_i^2$ depends on neither $\Pi$ nor $\mathcal{S}^{est}$:
\begin{equation}\label{eq:MSE_treatment_effect}
\text{MSE}_\tau(\mathcal{S}^{te}, \mathcal{S}^{est}, \Pi)
= \frac{1}{N^{te}} \sum_{i \in \mathcal{S}^{te}}
\Bigl[\bigl(\tau_i - \hat{\tau}(X_i; \mathcal{S}^{est}, \Pi)\bigr)^2 - \tau_i^2\Bigr].
\end{equation}
Equation~\eqref{eq:MSE_treatment_effect} is still not feasible in step~(\ref{itm:second}):
the leaf estimates $\hat\tau(X_i;\mathcal{S}^{est},\Pi)$ depend on $\mathcal{S}^{est}$, which has not been used during split selection, and the individual $\tau_i$ are unobserved on $\mathcal{S}^{te}$. To make the criterion feasible, \citet{athey_recursive_2016} take expectations of \eqref{eq:MSE_treatment_effect} over \((\mathcal S^{te},\mathcal S^{est})\) to obtain the expected mean squared error \(\mathrm{EMSE}_\tau(\Pi)\). They then construct a training-sample estimator of \(-\mathrm{EMSE}_\tau(\Pi)\) by replacing the unknown leaf CATE \(\tau(l;\Pi)\) with the training-sample difference-in-means and by subtracting the corresponding variance correction. This yields the honest splitting criterion
\begin{equation}\label{eq:MSE_CT_H}
\hat{Q}_\tau^{\,\text{MSE, H}}(\mathcal{S}^{tr}, \Pi)
= \sum_{l \in \Pi} \frac{N_l^{tr}}{N^{tr}}
\left[\hat{\tau}(l;\mathcal{S}^{tr})^2
- \left(\frac{1}{N_l^{tr}} + \frac{1}{N_l^{est}}\right)
\!\left(\frac{\hat{S}^2_{1}(l;\mathcal{S}^{tr})}{p}
+ \frac{\hat{S}^2_{0}(l;\mathcal{S}^{tr})}{1-p}\right)\right],
\end{equation}
which is maximized over $\Pi$. Here, \(\hat\tau(l;\mathcal S^{tr})=\bar Y^{tr}_{1,l}-\bar Y^{tr}_{0,l}\) is the
training-sample CATE estimate, and
\(\hat S_d^2(l;\mathcal S^{tr})\) is the within-leaf outcome variance for treatment status \(d\). The factor \(N_l^{tr}/N^{tr}\) estimates the leaf probability, while
\(N_l^{est}\) denotes the estimation-sample leaf size associated with leaf \(l\) under \(\Pi\).
The treatment share $p$ equals the known propensity score in the randomized designs considered here and should not be confused with the leaf probability.
The first term in \eqref{eq:MSE_CT_H}, $\hat\tau(l;\mathcal{S}^{tr})^2$
weighted by $N_l^{tr}/N^{tr}$, is the in-sample heterogeneity reward: splits
that produce leaves with larger differences in CATEs increase this term. The second term is the variance penalty. It corrects the upward bias of the plug-in $\hat\tau(l;\mathcal{S}^{tr})^2$ and discounting heterogeneity that will be hard to detect on
$\mathcal{S}^{est}$. The main contribution of this paper
is the modification of \eqref{eq:MSE_treatment_effect} and \eqref{eq:MSE_CT_H} by replacing its mean-based anchor with median-based alternatives that are less sensitive to heavy-tailed or skewed outcome distributions.
Although every quantity in \eqref{eq:MSE_CT_H} is computed on $\mathcal{S}^{tr}$, the penalty depends on the estimation-sample leaf size \(N_l^{est}\), which determines how strongly estimation variance is penalized. The observations in $\mathcal{S}^{est}$ are used only later for honest leafwise estimation in step~(\ref{itm:third}). The relation between \eqref{eq:MSE_CT_H} and \eqref{eq:MSE_treatment_effect} is derived in detail in Appendix~\ref{appendix:Honest_Splitting_CT}. We further note that the superscript ``H'' on $\hat Q^{\,\text{MSE, H}}_\tau$ marks the splitting criterion as honest due to the variance penalty derived under the EMSE objective. An adaptive criterion
\(\hat Q^{\,\mathrm{MSE,A}}_\tau\) omits this term. This is conceptually different from honest leaf estimation, which uses an independent $\mathcal{S}^{est}$ in step~(\ref{itm:third}) to compute $\hat\tau(l;\mathcal{S}^{est})$. Honest leaf estimation delivers per-leaf conditional unbiasedness for $\tau(l;\Pi)$, given the selected partition.
So far, we have focused on split selection and leafwise estimation for a single honest causal tree. The criterion in \eqref{eq:MSE_CT_H} is defined at the tree level. In practice, however, single trees can be unstable because small changes
in the data may lead to different partitions and high-variance leaf estimates. For this reason, the estimator we
evaluate in Sections~\ref{sec:simulations} and~\ref{sec:empappl} is the
corresponding causal forest of \citet{wager_estimation_2018}: an ensemble of
$B$ honest causal trees grown on random subsamples of the data and aggregated
by averaging,
\begin{equation}\label{eq:causal_forest}
\hat{\tau}(X_i) \;=\; \frac{1}{B}\sum_{b=1}^{B} \hat{\tau}_b(X_i),
\end{equation}
where $\hat\tau_b(X_i)$ is the prediction of the $b$-th tree, built and
re-estimated according to steps~(\ref{itm:first})--(\ref{itm:fourth}) on its
assigned subsample. The splitting rules considered in this paper (the mean-based
$\hat{Q}^{\,\text{MSE, H}}_\tau$ in \eqref{eq:MSE_CT_H}, the median-based MSD
in Section~\ref{sec:MSD}, and the heuristic MAD and LMS in
Section~\ref{sec:MAD_LMS}) differ only in step~(\ref{itm:second}) of each tree. Honest leaf estimation, subsampling, and forest aggregation are kept fixed across methods.
Aggregating the CATE estimates in \eqref{eq:causal_forest} over the full sample gives a simple plug-in estimator of the ATE in \eqref{eq:ATE},
\begin{equation} \label{eq:diff_estimator_ATE}
\hat{\tau}_{\text{plug-in}}
= \frac{1}{N} \sum_{i \in \mathcal{S}} \hat{\tau}(X_i),
\end{equation}
which is consistent but not necessarily efficient. For the ATE comparisons in Sections~\ref{sec:simulations} and~\ref{sec:empappl} we instead use the augmented inverse probability weighting (AIPW) estimator \citep{athey_policy_2021, robins_estimation_1994, chernozhukov_doubledebiased_2018}, detailed in Appendix~\ref{appendix:AIPW_ATE_CF}.
The pointwise consistency and asymptotic normality of the causal forest CATE estimates in \eqref{eq:causal_forest} are established in Theorem~4.1 of \citet{wager_estimation_2018},
which carries the generic forest central limit theorem of their Theorem~3.1 into the potential outcomes setting under unconfoundedness and overlap. In particular the trees must be honest, $\alpha$-regular, random-split,
and symmetric in the sense of \citet{wager_estimation_2018}, and the subsample size must scale appropriately.
We note that there is recent literature that scrutinizes these conditions of \citet{wager_estimation_2018}. \citet{cattaneo_honest_2025} show that adaptive causal trees generate unbalanced leaves with non-vanishing probability and argue that fixed-fraction $\alpha$-regularity may be incompatible with standard recursive partitioning, so that honesty alone may not guarantee uniform convergence. By contrast,
\citet{bladt_consistency_2026} show that honest trees remain consistent when leaves localize in feature space while retaining enough observations for stable within-leaf estimation. We treat
honesty and $\alpha$-regularity as maintained assumptions, enforced at the level of the implementation in Appendix~\ref{appendix:implementation}, and note that
the operative practical condition is the leaf-size regime rather than the splitting rule.
Replacing the mean-based criterion with the Median Squared Deviation (MSD), Median Absolute Deviation (MAD), or Least Median of Squares (LMS) rule introduced in the next section changes only the split-selection criterion in
step~(\ref{itm:second}). The honest sample split, subsampling scheme, and difference-in-means leaf estimator in step~(\ref{itm:third}) remain unchanged.\footnote{The balance constraint is Definition~4 of \citet{wager_estimation_2018}, which keeps $\alpha$-regularity intact under criteria whose unconstrained optimum can be highly unbalanced, most notably LMS.}
Therefore, their Theorem~4.1 continues to apply. The result relies on the finite-moment conditions of \citet{wager_estimation_2018}, which the heavy-tailed but finite-variance designs of Section~\ref{sec:simulations} satisfy, so the robustness we target lies in the stability of split selection under heavy-tailed or skewed outcome distributions. Appendix~\ref{appendix:implementation} documents that the implemented forest satisfies the conditions of Theorems~3.1 and~4.1 of \citet{wager_estimation_2018} and identifies the mechanism enforcing it, so the transfer rests on properties of the protocol that we briefly discuss there.
\section{Median-based Splitting Rules}
\label{sec:MedianSplitting}
We introduce median-based splitting rules as robust alternatives to the MSE-based criteria in Equations~\eqref{eq:MSE_treatment_effect} to \eqref{eq:MSE_CT_H}. \cite{roy_robustness_2012} proposed median-based splits to robustify regression forests for outcome prediction with outliers, while \cite{athey_recursive_2016} proposed the causal tree framework in Section \ref{sec:causal_trees} to allow for estimation and inference of CATE. We combine both approaches. This idea is also related to the robust M-estimation perspective of \citet{Huber1964},
where the target estimand can remain unchanged while the loss function used for estimation is modified to improve finite-sample stability (see also \citet{Hastie2009}). We adapt this logic to the causal-tree setting, the target estimand remains the CATE in \eqref{eq:CATE}, while the split-selection objective is modified by replacing its mean-based rule with median-based quantities, to quantify treatment effect heterogeneity under irregular outcome distributions. Because changing the criterion affects only which partition is chosen, and not how effects are estimated within leaves, the leaf-level unbiasedness and the forest inference of Section~\ref{sec:causal_trees} are unaffected by the choice of the splitting rule. This is a property of the honest leaf estimator, not of the splitting criterion in particular, and it holds for all considered rules. Whether the MSD criterion itself reproduces the honest variance-bias decomposition is a separate question, discussed in Section~\ref{sec:MSD}.
Before introducing the splitting rules, let us define a median-based analogue of the average treatment effect, the median treatment effect (MTE), as
\begin{align} \label{eq:MTE}
\tau_{\operatorname{med}} &\coloneqq \operatorname{med} \tau_i = \operatorname{med} (Y_i(1) - Y_i(0)),
\end{align}
Unlike the estimation of the ATE in \eqref{eq:ATE}, the estimation of the MTE in Definition \eqref{eq:MTE} is particularly challenging, because it is generally not identified from the marginal distributions of \(Y_i(1)\) and \(Y_i(0)\) alone. In particular, it depends on the joint distribution of \((Y_i(1),Y_i(0))\), while only one potential outcome is observed for each unit and crucially, two different joint distributions can share the exact same marginal distributions yet produce very different $\tau_{\operatorname{med}}$. Thus, additional assumptions on the dependence structure of the potential outcomes are required \citep{Addanki2024LimitsOA}.
The same issue remains for the conditional median treatment effect (CMTE) adapted to our leaf-wise potential outcomes setting and defined as
\begin{align} \label{eq:CMTE}
\tau_{\operatorname{med}}(X_i) &\coloneqq \operatorname{med}\bigl(Y_i(1) - Y_i(0)\mid X_i \bigr).
\end{align}
Since $\tau_i = Y_i(1) - Y_i(0)$ is not observable, we cannot directly compute the median of these unobserved quantities. The key distinction to keep in mind is that
\begin{align}\label{eq:median_of_differences}
\tau_{\operatorname{med}}(X_i) = \operatorname{med} (Y_i(1) - Y_i(0) \mid X_i )\neq \operatorname{med}\bigl(Y_i(1) \mid X_i \bigr) - \operatorname{med}\bigl(Y_i(0) \mid X_i \bigr),
\end{align}
which is the difference of conditional medians. The difference of conditional medians is identifiable from observed data \citep{leqi_median_2022}, but as \citet{Addanki2024LimitsOA} argue, it can diverge substantially from the median of individual differences because the coupling between $Y_i(1)$ and $Y_i(0)$ conditional on $X_i$ is unobservable.
\subsection{The Hodges–Lehmann estimator as splitting anchor}
Instead of using the right-hand side of \eqref{eq:median_of_differences} directly as a splitting target, we follow the literature on the Hodges--Lehmann (HL) estimator and the Wilcoxon rank-sum statistic \citep{hodges_estimates_1963, hoyland_HL, hollander_nonparametric_2014}. The HL estimator is most naturally interpreted as a rank-based estimator of a location-shift parameter. We therefore first establish its leafwise shift properties.
For a fixed \(\Pi\), define the leafwise median treatment effect in
the leaf containing \(x\) as
\[\tau_{\operatorname{med}}(x;\Pi):=\operatorname{med}\bigl(Y_i(1)-Y_i(0)\mid X_i\in l(x;\Pi)\bigr).\]
The connection between the shift parameter, $\tau_{\operatorname{med}}(x;\Pi)$ or the CATE, as well as its role in the subsequent splitting criteria, requires additional assumptions, which are discussed in the following subsections.
Let $\mathcal{S}_{1}$ and $\mathcal{S}_{0}$ denote the subsamples for the treatment and control group of sample $\mathcal{S}$, with corresponding sample sizes $N_{1}$ and $N_{0}$. Given a tree $\Pi$ and test point $x$, we define the HL estimator adapted to the leaf-wise potential outcomes setting as
\begin{equation}\label{eq:HL-leafwise}
\hat\tau_{\mathrm{HL}}(x;\Pi)
=
\operatorname{med}\Bigl\{
Y_j-Y_m
:\,
j\in\mathcal S_1,\;
m\in\mathcal S_0,\;
X_j\in l(x;\Pi),\;
X_m\in l(x;\Pi)
\Bigr\}.
\end{equation}
To study the large-sample behavior of $\hat\tau_{\mathrm{HL}}(x;\Pi)$ in the leafwise setting, we
first introduce the following assumptions.
\begin{assump}\label{ass:HL_shift}
For a fixed $\Pi$ and a target point $x$, let $l(x;\Pi)$ be the leaf containing $x$.
Define the leafwise population conditional distributions
\[
F_{1,l}(z) := Pr\bigl(Y_i(1)\le z \mid X_i\in l(x;\Pi)\bigr),
\qquad
F_{0,l}(z) := Pr\bigl(Y_i(0)\le z \mid X_i\in l(x;\Pi)\bigr).
\]
Assume that these distributions are absolutely continuous, with corresponding densities $f_{1,l}$ and
$f_{0,l}$. In addition, suppose that:
\begin{itemize}
\item[$i)$] the leafwise potential-outcome distributions differ by a location shift, that is,
there exists a leaf-specific parameter $\Delta(x;\Pi)\in\mathbb{R}$ such that
\[
F_{1,l}(z)=F_{0,l}(z-\Delta(x;\Pi))
\qquad \text{for all } z\in\mathbb{R};
\]
\item[$ii)$] under condition~(i), the common shifted density \(f_l\),
defined by
\[
f_{0,l}(z)=f_l(z),
\qquad
f_{1,l}(z)=f_l(z-\Delta(x;\Pi)),
\]
is continuous and satisfies
\[
0<\int f_l^2(u)\,du<\infty;
\]
\item[$iii)$] as the total sample size tends to infinity, the leaf containing $x$ contains
diverging numbers of treated and control observations, with
\[
N_{0,l(x;\Pi)} + N_{1,l(x;\Pi)} \to \infty
\qquad\text{and}\qquad
\frac{N_{0,l(x;\Pi)}}{N_{0,l(x;\Pi)} + N_{1,l(x;\Pi)}} \to \lambda \in (0,1).
\]
\end{itemize}
\end{assump}
Condition (i) in Assumption \ref{ass:HL_shift} is the main restriction. This condition requires the treatment to shift the leafwise outcome distribution without changing its shape. Under this restriction, the treatment-control comparison is summarized by the scalar shift parameter \(\Delta(x;\Pi)\), making the Hodges--Lehmann estimator a natural rank-based anchor. The shift model holds automatically under a within-leaf constant effect (Lemma~\ref{lemma:Delta_equals_tau_leaf}), and we return in Section~\ref{sec:MSD} to the weaker reading of the criterion when it does not.
The following results, proved in the appendix, establish that $ \hat\tau_{\mathrm{HL}}(x;\Pi)$ is a well-behaved anchor for our splitting criteria.
\begin{theorem}\label{thm:theorem_consistency}
Under Assumption~\ref{ass:HL_shift}, the Hodges--Lehmann estimator
$\hat\tau_{\mathrm{HL}}(x;\Pi)$ associated with the Wilcoxon rank-sum statistic in the
Mann--Whitney form is asymptotically normal and centered around
$\Delta(x;\Pi)$. In particular, $\hat\tau_{\mathrm{HL}}(x;\Pi)$ is a consistent estimator of
$\Delta(x;\Pi)$.
\end{theorem}
\begin{proof}
See Appendix \ref{appendix:proofs}.
\end{proof}
In practice, recursive partitioning produces progressively smaller leaves, so finite-sample behavior also matters. A useful property is median unbiasedness, meaning that the
estimator is equally likely to overestimate or underestimate its target. The following lemma gives a general
bound.
\begin{lemma} \label{lem:median_unbiased}
Let $\mathcal W_l=t(Y^0_l,Y^1_l)$ be the leafwise Wilcoxon rank-sum statistic. If
\[
Pr\bigl(t(Y^0_l,Y^1_l)=\xi_l \mid \Delta(x;\Pi)=0\bigr)=\delta,
\]
where $\xi_l$ is a symmetry point, then
\[
\frac12-\frac{\delta}{2}
\le
Pr\!\left(\hat\tau_{\mathrm{HL}}(x;\Pi)\le \Delta(x;\Pi)\right)
\le
\frac12+\frac{\delta}{2}.
\]
\end{lemma}
\begin{proof}
See Appendix \ref{appendix:proofs}.
\end{proof}
Lemma~\ref{lem:median_unbiased} shows that exact median unbiasedness obtains whenever the symmetry point is unattainable, so whenever $\delta=0$. When the symmetry point is attainable, the same lemma quantifies the possible deviation from exact median unbiasedness by the size of $\delta$. Theorem~\ref{thm:theorem_median_unbiased} applies this observation to show that $\hat\tau_{\mathrm{HL}}(x;\Pi)$ is median unbiased for $\Delta(x;\Pi)$.
\begin{theorem}\label{thm:theorem_median_unbiased}
For a fixed $\Pi$ and a target point $x$, let $l(x;\Pi)$ be the leaf containing $x$.
Suppose that, under the null shift $\Delta(x;\Pi)=0$, the distribution of the leafwise Wilcoxon rank-sum statistic $\mathcal W_l$ is symmetric about
\[
\xi_l=\frac{N_{1,l}N_{0,l}}{2}.
\]
Then the Hodges--Lehmann estimator $\hat\tau_{\mathrm{HL}}(x;\Pi)$ is median unbiased for
$\Delta(x;\Pi)$ when $N_{1,l}N_{0,l}$ is odd, in the sense that
\[
Pr\!\left(\hat\tau_{\mathrm{HL}}(x;\Pi)\le \Delta(x;\Pi)\right)=\frac12.
\]
If \(N_{1,l}N_{0,l}\) is even, the estimator is approximately median unbiased, in the sense that its deviation from exact median unbiasedness is bounded by \(\delta/2\), with \(\delta\) defined in Lemma~\ref{lem:median_unbiased}.
\end{theorem}
\begin{proof}
See Appendix \ref{appendix:proofs}.
\end{proof}
\subsection{Motivating Example: Splitting under Heavy Tails}
\label{sec:motivation}
Consider a leaf $l$ of a tree $\Pi$ grown on $\mathcal{S}^{tr}$ with $N_{1,l} = N_{0,l} = 10$
observations per treatment group, and let the data follow a constant additive treatment effect,
\begin{align}\label{eq:sharp-null}
\begin{split}
Y_i(0) &= 3 + \varepsilon_i, \\
Y_i(1) &= Y_i(0) + \tau, \qquad \tau = 2.
\end{split}
\end{align}
Under \eqref{eq:sharp-null} the conditional average and conditional median
treatment effects in $l$ both equal $\tau$, and the leafwise shift model of Assumption \ref{ass:HL_shift}(i) holds with $\Delta(x; \Pi) = \tau$ by Lemma \ref{lemma:Delta_equals_tau_leaf}. Thus, the example isolates a simple setting in which the HL anchor is well defined.
What is at stake here is not how leaf effects are estimated, but which partition is selected.
Honest leaf estimation in step (\ref{itm:third}) reports the difference in means on the held-out $\mathcal{S}^{est}$, and the median-based rules of this section leave that step untouched. The quantity a heavy-tailed draw can corrupt is the training-sample anchor $\hat\tau(l ; \mathcal{S}^{tr})$ entering the splitting criterion \eqref{eq:MSE_CT_H}.
We therefore compare two candidate anchors for that criterion, the training-sample difference-in-means estimator and the HL estimator,
\begin{align}\label{eq:two-estimators}
\begin{split}
\hat\tau(l; \Pi) &= \bar Y_{1,l} - \bar Y_{0,l}, \\
\hat\tau_{\mathrm{HL}}(l; \Pi)
&= \operatorname{med} \{ Y_{j,l} - Y_{m,l}: 1 \le j \le N^{tr}_{1,l},\,
1 \le m \le N^{tr}_{0,l} \},
\end{split}
\end{align}
where $\hat\tau_{\mathrm{HL}}(l; \Pi)$ is the median of all $N^{tr}_{1,l} \times N^{tr}_{0,l} = 100$ pairwise treated--control differences within $l$ and $\bar Y_{1,l}, \bar Y_{0,l}$ the corresponding training-sample outcome averages.
The difference-in-means estimator targets $\tau$ because the treatment and control means differ by $\tau$ in population.
The HL estimator targets the same quantity because each pairwise treated--control difference equals $Y_{j,l}-Y_{m,l}=\tau+(\varepsilon_j-\varepsilon_m)$, and $\varepsilon_j-\varepsilon_m$ is symmetric about zero when $\varepsilon_j$ and $\varepsilon_m$ are independent draws from the same baseline distribution.
Hence, the population median of the pairwise differences is $\tau$, even if the baseline distribution itself is skewed. Both anchors in \eqref{eq:two-estimators} are therefore consistent for $\tau$ but they might differ in finite samples.
To see this, let the baseline be heavy-tailed with $\varepsilon_i$ in \eqref{eq:sharp-null} as
\begin{align}
\varepsilon_i \stackrel{\mathrm{iid}}{\sim} t_\nu
\end{align}
and $\nu = 3$, in particular $\operatorname{Var}(\varepsilon) = 3$.
A single observation that is $c$ units above the typical treated outcome increases the treatment-group mean by $c/N^{tr}_{1,l}$. For $N^{tr}_{1,l}=10$, this increase equals $c/10$ and choosing $c=10$ gives an increase of $1$, so $\hat\tau(l;\mathcal{S}^{tr})$ moves from $\tau=2$ to approximately $3$.
One such draw suffices, because the sample mean has breakdown zero.
The same draw enters only $N^{tr}_{0,l}=10$ of the $100$ pairwise differences, so the median defining $\hat\tau_{\mathrm{HL}}(l;\Pi)$ is still taken over $90$ uncontaminated values and is left essentially unchanged \citep{hollander_nonparametric_2014}. We note that honesty does not fully close this gap. Step (\ref{itm:third}) removes the training-sample draw from the reported effect by re-estimating on $\mathcal{S}^{est}$, but it cannot remove it from $\Pi$, because the split has already been made and the partition is not revised. The mean-based criterion in \eqref{eq:MSE_CT_H} rewards apparent leaf-heterogeneity through $\hat\tau(l;\mathcal{S}^{tr})^2$ while the variance penalty does not offset this increase. Appendix \ref{app:motivation} provides further details for this behavior based on our setup in \eqref{eq:sharp-null}. Because \eqref{eq:MSE_CT_H} is maximized over a large number of candidate splits, and the maximum of a collection of heavy-tailed statistics is governed by its tails, the split that isolates the extreme draw is preferentially selected, yielding trees that track tail noise rather than genuine heterogeneity. Anchoring the split criterion at $\hat\tau_{\mathrm{HL}}(l;\Pi)$ suppresses this behavior while still targeting $\tau$ under \eqref{eq:sharp-null}.
\subsection{Median Squared Deviation as a Robust Splitting Criterion}
\label{sec:MSD}
We now construct the Median Squared Deviation (MSD) criterion, the main median-based alternative to the honest MSE criterion in \eqref{eq:MSE_CT_H}. The MSD rule modifies only the split-selection objective in step~(\ref{itm:second}) by replacing the mean-based leaf anchor in the sample analog of \eqref{eq:MSE_treatment_effect} with the leafwise Hodges--Lehmann estimator from \eqref{eq:HL-leafwise}. Honest leaf estimation remains difference-in-means estimation on \(\mathcal S^{est}\), so the final target remains the leafwise CATE \(\tau(l;\Pi)\). The idea is that, when the conditional outcome distributions are
heavy-tailed, $\hat\tau(l;\mathcal S^{tr},\Pi)$ is volatile and the within-leaf variances that enter the honest penalty are inflated, both of which distort split selection.
The HL estimator is consistent for the population location shift $\Delta(x;\Pi)$, and under the location-shift conditions in Appendix~\ref{appendix:HL_estimator} we have $\Delta(x;\Pi) = \tau(l;\Pi)$.
When those conditions fail, the Hodges--Lehmann anchor estimates the location shift $\Delta(x;\Pi)$, which need not equal $\tau(l;\Pi)$, so the criterion no longer targets the CATE objective exactly and may rank splits suboptimally. This affects only partition selection.
Because leaf effects are estimated by difference in means on $\mathcal{S}^{est}$ under every splitting rule that we consider in this paper, the resulting CATE estimates remain unbiased for $\tau(l;\Pi)$ on whatever partition is chosen, so a misaligned anchor costs estimation efficiency through a coarser or less informative partition, not validity of the estimate. Scenarios~S3 and~S4 in Section~\ref{sec:simulations} illustrate this violation, and the MSD forest retains the lowest bias and root-mean-square error there, which bounds the practical size of this efficiency cost on the designs we study.
Starting from the modified, infeasible MSE in Equation~\eqref{eq:MSE_treatment_effect}, we
construct a median squared deviation splitting rule by replacing the leafwise treatment effect
estimator $\hat\tau(x;\mathcal S^{est},\Pi)$ with the Hodges--Lehmann estimator
$\hat\tau_{\mathrm{HL}}(x;\mathcal S^{est},\Pi)$. The previous subsection established that
$\hat\tau_{\mathrm{HL}}(x;\Pi)$ estimates the leafwise location shift $\Delta(x;\Pi)$. This gives rise to the infeasible median squared deviation (MSD) splitting criterion based on the MSE in Equation~\eqref{eq:MSE_treatment_effect}:
\begin{equation}\label{eq:mse_tau_med}
\mathrm{MSE}_{\tau,\mathrm{med}}(\mathcal{S}^{te},\mathcal{S}^{est},\Pi)
\coloneqq \frac{1}{N^{te}} \sum_{i \in \mathcal{S}^{te}}
\Bigl[\bigl(\tau_i - \hat{\tau}_{\mathrm{HL}}(X_i; \mathcal{S}^{est}, \Pi)\bigr)^2 - \tau_i^2\Bigr],
\end{equation}
where $\hat{\tau}_{\mathrm{HL}}(X_i;\mathcal{S}^{est},\Pi)$ is the leaf-constant HL estimate assigned to unit~$i$. The subtraction of $\tau_i^2$ again preserves the relative ranking of candidate trees since $\tau_i^2$ does not depend on $\Pi$ or $\mathcal{S}^{est}$, as in \eqref{eq:MSE_treatment_effect}. Mirroring the honest derivation in Appendix~\ref{appendix:Honest_Splitting_CT} and using
Theorems~\ref{thm:theorem_consistency} and~\ref{thm:theorem_median_unbiased}, we obtain the honest objective
\begin{align}
\label{eq:honest_median}
\begin{split}
\hat{Q}^{\,\text{MSD, H}}_\tau(\mathcal{S}^{tr}, \Pi)
&= \sum_{l \in \Pi} \frac{N_l^{tr}}{N^{tr}}
\left[\hat{\tau}_{\mathrm{HL}}(l;\mathcal{S}^{tr})^2
- \left(\frac{1}{N_l^{tr}} + \frac{1}{N_l^{est}}\right) \widehat{V}_l \right]
\end{split}
\end{align}
where $\hat\tau_{\mathrm{HL}}(l;\mathcal{S}^{tr})$ is the within-leaf HL estimate on the
training sample and $\widehat V_l$ estimates its asymptotic variance. The $(1/N_l^{tr}+1/N_l^{est})$ penalty weight is identical to that of \eqref{eq:MSE_CT_H}. Further details on the construction are provided in Appendix~\ref{appendix:Honest_Splitting_MSD}.
Two features distinguish \eqref{eq:honest_median} from the mean-based honest criterion $\hat Q^{\,\text{MSE,H}}_\tau$ in \eqref{eq:MSE_CT_H}.
First, the heterogeneity reward is anchored on the training-sample estimate $\hat\tau_{\mathrm{HL}}(l;\mathcal{S}^{tr})^2$ rather than the difference in means.
Second, the variance penalty replaces the treatment and control group sample variances $\hat S^2_d / N^{est}_l$ with $\widehat V_l$, an estimated variance appropriate for a median-based estimator.
Because $\widehat V_l$ requires kernel density evaluations and is costly within a tree-growing loop, our implementation maximizes an auxiliary splitting criterion to preserve the qualitative effect of penalizing leaves with weak effective signal-to-noise while avoiding density estimation at every candidate split.
\begin{align}\label{eq:Q_MSD_code_main}
\begin{split}
\widehat Q^{\,\text{MSD}}_\tau\left(\mathcal{S}^{tr}, \Pi\right)
&=\sum_{l\in\Pi}\frac{N_l^{tr}}{N^{tr}}
\!\left[\hat\tau_{\mathrm{HL}}(l;\mathcal{S}^{tr})^2\,-\Bigl(
\hat\tau_{\mathrm{HL}}\left(l;\mathcal{S}^{tr}\right)-\hat\tau\left(l;\mathcal{S}^{tr}\right)\Bigr)^2\right]
\end{split}
\end{align}
The honest MSD criterion in \eqref{eq:honest_median} reproduces the variance-bias decomposition of the mean-based criterion in \eqref{eq:MSE_CT_H}, with the Hodges--Lehmann variance in place of the outcome-variance penalty. The implemented criterion in \eqref{eq:Q_MSD_code_main} is the leading-order surrogate of \eqref{eq:honest_median}. This surrogate drops the explicit variance penalty in favor of the squared difference between Hodges--Lehmann estimates and standard mean-based leafwise estimates while reproducing the decomposition only up to a residual of order $O(1/N_l^{est})$, derived in Appendix~\ref{appendix:Honest_Splitting_MSD}. It is this surrogate that the simulations of Section~\ref{sec:simulations} evaluate.
\subsection{Two Baseline Robust Splitting Criteria}
\label{sec:MAD_LMS}
We complement the previous discussion with two additional splitting rules that are robust by construction but differ from MSD in their theoretical role. First, we consider the median absolute deviation (MAD) rule that uses the absolute difference between the mean-based leaf estimate and the Hodges--Lehmann leaf estimate, and second a least median squares (LMS) rule built on fit-based
partitioning \citep{ZeileisModelBasedRP, athey_recursive_2016, roy_robustness_2012}. We treat both as heuristics and report them in
Section~\ref{sec:simulations} as comparison rules for MSD, since they do not yield the same honest EMSE decomposition. Appendix~\ref{appendix:MAD_LMS} provides further details about MAD and LMS.
Following \citet{roy_robustness_2012}, the absolute-deviation analog of the
infeasible CATE-targeting MSE in \eqref{eq:MSE_treatment_effect} is
\begin{equation}\label{eq:infeasible_MAD}
\mathrm{MAD}_\tau(\mathcal{S}^{tr})
\;:=\;\frac{1}{N^{tr}}\sum_{i\in\mathcal{S}^{tr}}
\bigl|\tau_i - \tau_{\mathrm{med}}(X_i)\bigr|,
\end{equation}
where $\tau_{\mathrm{med}}(X_i)$ is the CMTE from \eqref{eq:CMTE}.
\eqref{eq:infeasible_MAD} is doubly infeasible: $\tau_i$ is unobserved (as in
the MSE), and $\tau_{\mathrm{med}}(X_i)$ is not point-identified without
restrictions on the joint distribution of potential outcomes
\citep{Addanki2024LimitsOA}. We obtain a feasible rule by replacing $\tau_i$ by
the mean-based plug-in $\hat\tau(X_i;\mathcal{S}^{tr},\Pi)$ and
$\tau_{\mathrm{med}}(X_i)$ by the leafwise HL estimator
$\hat\tau_{\mathrm{HL}}(X_i;\mathcal{S}^{tr},\Pi)$:
\begin{equation}\label{eq:MAD_split}
\widehat Q^{\,\mathrm{MAD}}_\tau(\mathcal{S}^{tr},\Pi)
\;:=\;\frac{1}{N^{tr}}\sum_{i\in\mathcal{S}^{tr}}
\bigl|\hat\tau(X_i;\mathcal{S}^{tr},\Pi)
-\hat\tau_{\mathrm{HL}}(X_i;\mathcal{S}^{tr},\Pi)\bigr|.
\end{equation}
Splits are chosen by minimizing \eqref{eq:MAD_split}, so the criterion favors partitions for which the mean-based and Hodges--Lehmann leaf estimates are close.
Large values of \eqref{eq:MAD_split} indicate leaves in which the mean-based treatment-effect estimate differs substantially from the rank-based location estimate, which can occur when tail observations influence the leaf mean.
The interpretation of \eqref{eq:MAD_split} differs from that of the infeasible criterion in \eqref{eq:infeasible_MAD}. The infeasible \eqref{eq:infeasible_MAD} is defined in terms of individual treatment effects $\tau_i$ and the
CMTE $\tau_{\mathrm{med}}(X_i)$, whereas the feasible \eqref{eq:MAD_split}
compares two leaf-level estimators of $\tau(l;\Pi)$ (the mean-based
$\hat\tau(X_i;\mathcal{S}^{tr},\Pi)$ and the location-shift-based $\hat\tau_{\mathrm{HL}}(X_i;\mathcal{S}^{tr},\Pi)$) on
$\mathcal{S}^{tr}$. The two objects coincide only under the location-shift
conditions of Appendix~\ref{appendix:HL_estimator}, where
$\Delta(x;\Pi) = \tau_{\mathrm{med}}(x) = \tau(l;\Pi)$. Outside that regime,
\eqref{eq:MAD_split} is a robustness diagnostic and not a sample analog of
\eqref{eq:infeasible_MAD}.
Next, we motivate LMS as a robust analogue of fit-based splitting. The adaptive fit-based criterion of \citet{ZeileisModelBasedRP}, mirroring
\citet{athey_recursive_2016}, selects splits by within-leaf
goodness-of-fit of an outcome model $\hat\mu(D_i,X_i;\mathcal{S}^{tr},\Pi)$
consisting of an intercept and a treatment indicator:
\begin{equation}\label{eq:MSE_fit_A}
\mathrm{MSE}^{\mathrm{fit\text{-}A}}_\mu(\mathcal{S}^{tr},\mathcal{S}^{tr},\Pi)
\;:=\;\frac{1}{N^{tr}}\sum_{i\in\mathcal{S}^{tr}}
\left[\left(Y_i-\hat\mu(D_i,X_i;\mathcal{S}^{tr},\Pi)\right)^2-Y_i^2\right].
\end{equation}
Replacing the mean squared residual within each leaf by the median squared residual \citep{rousseeuw_least_1984} yields the LMS criterion,
\begin{equation}\label{eq:LMS_split}
\widehat Q_\mu^{\mathrm{LMS}}
(\mathcal S^{tr},\Pi)
:=
\sum_{l\in\Pi}
\frac{N_l^{tr}}{N^{tr}}
\underset{i\in\mathcal S_l^{tr}}{\operatorname{med}}
\left[
\left(
Y_i-
\hat\mu(D_i,X_i;\mathcal S^{tr},\Pi)
\right)^2
\right],
\end{equation}
where the median squared residual is computed separately within each leaf. Candidate splits are chosen by minimizing \eqref{eq:LMS_split}, or equivalently by maximizing the reduction in the criterion relative to the parent node.
Like \eqref{eq:MSE_fit_A}, LMS targets fit of $Y$ given $(D,X)$ within leaves, not treatment effect heterogeneity and it inherits the limitations of fit-based partitioning discussed in
\citet[\S 3.1.3]{athey_recursive_2016}.
Similar to the transformed-outcome trees described in
\citet{athey_recursive_2016}, MAD and LMS are best viewed as benchmark splitting rules with a more practical motivation. They do not mimic the honest EMSE derivation of the causal tree criterion, but they provide simple robust alternatives that can be implemented within the same causal forest algorithm. MAD focuses split selection on the absolute deviation between mean-based and rank-based leaf estimates, thereby highlighting partitions whose estimated effects are sensitive to tail observations. LMS robustifies fit-based partitioning by replacing the
mean of squared residuals with their median.
\section{Simulation Study}
\label{sec:simulations}
The simulation study builds on the non-linear treatment effect setting of
\citet{wager_estimation_2018} and extends the data-generating process to analyze the conditions of
Assumption~\ref{ass:HL_shift}. We generate the potential outcomes $Y_i(0), Y_i(1)$ as
\begin{align}\label{eq:PO_dgp}
\begin{split}
Y_i(0) &= \varepsilon_i,\\
Y_i(1) &= Y_i(0) + \tau(X_i) + U_i,
\end{split}
\end{align}
where $\varepsilon_i$ is the outcome noise, $\tau(X_i)$ is the CATE from Equation \eqref{eq:CATE}, and $U_i$ is a mean-zero individual treatment effect component independent of $(X_i,\varepsilon_i)$. The treatment indicator $D_i$ is randomized as $D_i\sim\operatorname{Bernoulli}(0.5)$ with $D_i\perp X_i$ and we consider a constant propensity score with $Pr(D_i=1|X_i=x)=0.5$.
The covariates are drawn independently as $X_i\sim\mathrm{Uniform}[0,1]^{K}$ as in \citet{wager_estimation_2018} and we use $K=10$ and a sample size of $N=1000$ as defaults throughout this section. In contrast to \citet{wager_estimation_2018}, we vary three components of the data-generating process across scenarios:
\begin{enumerate}
\item[(a)] The noise distribution $\varepsilon_i$, either standard
Gaussian or Student's $t_3$,
\begin{equation*}
\varepsilon_i\sim\mathcal{N}(0,1)
\qquad\text{or}\qquad
\varepsilon_i\sim t_3.
\end{equation*}
\item[(b)] The non-linear treatment effect function $\tau(x)$. In the
non-sparse setting, $\tau(x)=\zeta(x_1)\,\zeta(x_2)$ with, in general,
\begin{equation*}
\zeta(v) = 1 + \frac{1}{1+\exp\!\left(-20\bigl(v-\tfrac{1}{3}\bigr)\right)},
\end{equation*}
such that $\tau(x)$ depends multiplicatively on the first two covariates, $x_1$ and $x_2$, only. In the sparse setting, we add
$10\cdot\mathds{1}(x_1^2+x_2^2>1.2^2)$ to $\tau(x)=\zeta(x_1)\,\zeta(x_2)$. The sparse region $\{x:x_1^2+x_2^2>1.2^2\}$ is the upper-right corner of the unit square in the first two coordinates, so the extreme responders are a minority subpopulation rather than a bulk effect.
\item[(c)] The ITE component $U_i$, set either to zero (so that
$Y_i(1)-Y_i(0)$ is deterministic given $X_i$) or drawn from a
mean-centered log-normal,
\begin{equation*}
U_i = 2\left(\exp(Z_i) - \exp(\tfrac{1}{2})\right),\qquad
Z_i\sim\mathcal{N}(0,1),
\end{equation*}
which is right-skewed with $\mathbb{E}[U_i]=0$ but
$\operatorname{med}(U_i)=2\!\left(1-\exp\!\left(\tfrac{1}{2}\right)\right)<0$. The non-zero specification places $U_i$ in the role of the unmeasured effect modifier studied by \citet{post2024flexible} and \citet{PostvandenHeuvel+2025}, in an additive form consistent with their setup. \citet{PostvandenHeuvel+2025} emphasize that flexible CATE estimators silently miss the resulting gap between ITE and CATE whenever $U_i$ is non-degenerate. The independence $U_i\perp(\varepsilon_i,X_i)$ imposed here corresponds to their identification setup, and
the log-normal shape of $U_i$ mirrors an illustrative ITE distribution in \citet{PostvandenHeuvel+2025}.
\end{enumerate}
We consider four scenarios, summarized in Table~\ref{tab:scenarios}.
Scenario \textbf{S1} serves as the benchmark design. It uses Gaussian
outcome noise, a smooth treatment effect with no sparse region, and no residual individual effect component, which places it closest to the idealized leafwise location shift of Assumption~\ref{ass:HL_shift}. The only departure from an exact location shift arises through the smooth variation of $\tau(X_i)$ inside each leaf, and this departure vanishes as the partition refines, so S1 is the setting in which the MSE criterion is expected to perform well.
Scenario \textbf{S2} differs from S1 only in replacing the Gaussian noise with Student's $t_3$ noise. Because $\operatorname{Var}(t_3)=3$, this both thickens the tails and triples the noise variance relative to the unit-variance Gaussian, so S2 stresses robustness on both counts while leaving the treatment effect structure unchanged. Condition~(i) of Assumption~\ref{ass:HL_shift} therefore continues to hold in the same asymptotic sense as in S1, since changing the noise distribution leaves the within-leaf shift structure unchanged.
Scenario \textbf{S3} retains Gaussian noise and $U_i = 0$ but introduces the sparse region of extreme responders in $\tau(x)$. Leaves that cross the responder boundary then contain a bimodal mixture of baseline and extreme effects, which produces sharp heterogeneity within the leaf and a clear departure from the location shift model. Because the boundary is fixed and the offending leaves shrink in measure as the forest refines, this violation is a property of the partition rather than of the data generating process.
Scenario \textbf{S4} keeps the smooth treatment effect of S1 but adds a skewed component $U_i$ with mean zero. Since $U_i$ is independent of $X_i$ with $\mathbb{E}[U_i] = 0$, the pointwise CATE remains $\tau(x)$, yet the skew inflates and reshapes the treated potential outcome distribution relative to the control distribution, so the mean and median individual treatment effects no longer coincide. In contrast to S3, this violation is structural and persists at every leaf as the sample size grows.\footnote{Besides the violation of condition (i) in Assumption \ref{ass:HL_shift}, scenario~S4 also violates the symmetry condition~(ii) of
Proposition~\ref{prop:Delta_equals_tau_equals_tauMed}, since the skewed $U_i$
makes the within-leaf distribution of $Y_i(1)-Y_i(0)$ asymmetric. Table \ref{tab:scenarios}
records condition~(i) in Assumption \ref{ass:HL_shift} only, for comparability across scenarios.}
\begin{table}[htbp]
\centering
\caption{Simulation scenarios.}
\label{tab:scenarios}
\begin{tabular}{lcccc}
\toprule
Scenario & $\varepsilon_i$ & $\tau(X_i)$ & $U_i$
& Violation of condition (i) in Assumption~\ref{ass:HL_shift} \\
\midrule
\textbf{S1} & $\mathcal{N}(0,1)$ & non-sparse & $0$
& none \\
\textbf{S2} & $t_3$ & non-sparse & $0$
& none \\
\textbf{S3} & $\mathcal{N}(0,1)$ & sparse & $0$
& within-leaf effect heterogeneity at boundary \\
\textbf{S4} & $\mathcal{N}(0,1)$ & non-sparse & skewed
& non-degenerate ITE inflates and skews $F_{1,l}$ \\
\bottomrule
\end{tabular}
\begin{notes}
The error term $\varepsilon_i$, the treatment effect function $\tau(X_i)$, and the ITE component $U_i$ determine the potential outcomes in \eqref{eq:PO_dgp}. The leafwise location-shift condition~(i) of Assumption~\ref{ass:HL_shift} holds only in S1 and S2. The regularity
conditions~(ii) and~(iii) hold by construction in every scenario.
\end{notes}
\end{table}
We illustrate the CATE and ATE results in Figures~\ref{fig:CATE_precision_s1s4} and~\ref{fig:ATE_precision}. Appendix~\ref{app:add_sim_res} reports sensitivity to sample size (Tables~\ref{tab:cate_samplesize} and~\ref{tab:ate_samplesize}) and to covariate dimension (Tables~\ref{tab:cate_covariatesize} and~\ref{tab:ate_covariatesize}). Further details of the design are presented in Appendix~\ref{appendix:MC_Sim} while Table \ref{tab:runtime} compares the algorithm runtime of causal forests with different splitting criteria. Confidence intervals around the causal forest CATE estimates use an analog of the bootstrap of little bags (BLB) \citep{sexton_standard_2009, athey_generalized_2019} with scaling. Appendix~\ref{appendix:BLB} presents the implementation and compares BLB with and without scaling in Table~\ref{tab:blb_scaling_ci}.
The MSD results throughout this section use the implemented surrogate criterion \eqref{eq:Q_MSD_code_main}, which approximates the honest criterion \eqref{eq:honest_median} to leading order. Therefore, the simulations evaluate the practical implementation of MSD rather than the exact formal honest criterion. Figure~\ref{fig:CATE_precision_s1s4} shows that this surrogate MSD splitting criterion attains the lowest absolute bias in all four scenarios. In S1 and S2, which are closest to the location-shift setting in Assumption~\ref{ass:HL_shift} and where the MSE splitting criterion is asymptotically optimal for $\tau(x)$, MSD nevertheless lowers root-mean-square error (RMSE) by roughly 13 to 15 percent relative to the MSE splitting criterion in Table~\ref{tab:cate_samplesize}. The MAD and LMS rules perform worse than both MSE and MSD. MAD penalizes differences between the mean-based and Hodges--Lehmann leaf estimates, but does not reward treatment-effect heterogeneity. LMS chooses splits based on outcome fit rather than treatment-effect heterogeneity. As a result, both rules are less directly aligned with variation in $\tau(x)$ than MSD.
Empirical coverage in Figure~\ref{fig:CATE_precision_s1s4} shows the sharpest contrast between MSD and the alternatives. MSD stays closest to the nominal $0.95$ level in the designs where MSE undercovers (S1, S3), while the MAD and LMS rules undercover already in the rather well-behaved scenario S1 around $0.70$. The undercoverage of MAD and LMS confirms that their corresponding BLB variance estimators are poorly calibrated, consistent with the absence of a clean identification result for those criteria. In S1, MSD overcovers at around $0.99$, and in S4 all four methods overcover, reflecting the additional variance contributed by the skewed $U_i$, which inflates absolute interval widths uniformly across methods. The price of MSD's conservative coverage is visible in the bottom row, where MSD intervals are systematically about $10$ to $15$ percent wider than MSE intervals. Across all four scenarios, MSD therefore trades a modest width premium for improved coverage, while simultaneously delivering the best point estimate precision.
\begin{figure}[htbp]
\centering
\caption{Precision and confidence interval coverage of CATE estimates from causal forests with different splitting criteria.}
\includegraphics[width=\textwidth]{boxplot_CATE_main.pdf}
\begin{notes}
Per-replication absolute bias (top row) of CATE estimates as well as empirical coverage (middle row) and width (bottom row) of 95\% confidence intervals for the CATE estimates from causal forests with different splitting criteria across $MC = 100$ Monte Carlo replications ($N = 1000$, $K = 10$). Each panel corresponds to one of the four simulation scenarios defined in Table~\ref{tab:scenarios}. Splitting criteria within the causal forest are: MSE in \eqref{eq:MSE_CT_H}, MSD in \eqref{eq:Q_MSD_code_main} and highlighted in orange, MAD in \eqref{eq:MAD_split}, and LMS in \eqref{eq:LMS_split}. Boxes span the interquartile range (IQR), whiskers extend to $1.5 \times \text{IQR}$ and dots mark outlying replications. Appendix \ref{app:add_sim_res} reports full results in tabular format.
\end{notes}
\label{fig:CATE_precision_s1s4}
\end{figure}
Part of MSD's advantage is not directly attributable to robustness. In S1, with Gaussian noise and no outliers to resist, one would expect the Hodges--Lehmann anchor to be marginally less efficient than the mean, since robust location estimators trade efficiency at the Gaussian model for resistance to contamination \citep{huber_robust_2009}. For the Hodges--Lehmann estimator this cost is small, with an asymptotic relative efficiency under normality that is close to that of the mean \citep{hodges_estimates_1963}. The superior performance of MSD in S1 therefore suggests that its conservative split selection, which avoids aggressive splitting in response to small-sample fluctuations, provides a regularization benefit that is distinct from tail robustness. The analysis does not disentangle the effects of regularization and robustness with a separate baseline. The observed advantage in S1 should therefore be interpreted as consistent with regularization rather than as direct evidence for the Hodges--Lehmann anchor. The contribution of robustness is instead indicated by changes in the margin: it is roughly constant across S1 and S2 and widens in S3, precisely the scenario in which asymmetric within-leaf contamination pulls mean-based splits toward spurious heterogeneity. A regularization effect would not concentrate in the contaminated design, so this increment is attributable to the robust anchor rather than to conservative splitting. In S4, the skewed but mean-zero $U_i$ leaves the pointwise CATE equal to $\tau(x)$, so the difference-in-means leaf estimates remain unbiased for $\tau(x)$ under every splitting rule. The skew inflates and reshapes $F_{1,l}$, and the HL splitting anchor may be biased for $\tau(l;\Pi)$ under the failed location shift. However, this affects only which splits are chosen, not the difference-in-means leaf estimates, so the MSD point estimates carry no median-shift bias and MSD is still superior, because its conservative splitting reduces estimation variance.
\begin{figure}[htbp]
\centering
\caption{Absolute bias of ATE estimates from causal forests and competing estimators.}
\label{fig:ATE_precision}
\includegraphics[width=\textwidth]{boxplot_ATE_main.pdf}
\begin{notes}
Per-replication absolute bias of ATE estimates from causal forests across $MC = 100$ Monte Carlo replications ($N = 1000$, $K = 10$). For CF (MSE) and CF (MSD), we use cross-fitted AIPW \citep{robins_estimation_1994, chernozhukov_doubledebiased_2018} with causal forest CATE estimates as explained in Appendix \ref{appendix:AIPW_ATE_CF}. Further reported ATE estimators are: the difference-in-means estimator (Diff.Mean), OLS with interaction (Lin.OLS) proposed in \cite{lin_agnostic_2013}, the rank-based regression-adjustment estimator (R.Adj) of \cite{ghosh_robustness_2026}, and the efficient-influence-function (EIF) and weighted-average-of-quantiles (WAQ) estimators used in \cite{Athey2023}. We provide further information about these competing estimators in Appendix \ref{appendix:ATE_estimators}. Scenarios S1--S4 are as defined in Table~\ref{tab:scenarios} and the splitting criteria within the causal forests are MSE in \eqref{eq:MSE_CT_H} and MSD (coloured in orange) in \eqref{eq:Q_MSD_code_main}. Boxes span the interquartile range (IQR) and the whiskers extend to $1.5 \times \text{IQR}$ while dots mark outlying replications.
\end{notes}
\end{figure}
For the ATE, the splitting rule within the causal forest is immaterial for precision and coverage results. CF (MSD) and CF (MSE) give near-identical ATE estimates and confidence interval widths across all four scenarios (Tables~\ref{tab:ate_samplesize} and~\ref{tab:ate_covariatesize} in Appendix \ref{app:add_sim_res}), since the splitting choice washes out once CATE estimates are aggregated by AIPW and is further discussed in Appendix \ref{appendix:AIPW_ATE_CF}. This holds true even for the noisier CATE estimates of CF (MAD) and (LMS), which are not reported here for brevity due to their similar ATE estimates, comparable to CF (MSD) and CF (MSE).
The informative contrast that we present here is with the competing estimators from the simulation study in \citet{ghosh_robustness_2026} and described in Appendix \ref{appendix:ATE_estimators}. The difference-in-means estimator and OLS with interactions \citep{lin_agnostic_2013} stay
centered on the mean ATE in every scenario, holding nominal coverage in S1 to
S3 and over-covering in S4 along with the causal forests, as the heavy $U_i$ inflates all variances. The rank- and quantile-based estimators (Rosenbaum's adjustment (R.Adj) \citep{ghosh_robustness_2026}, efficient-influence-function (EIF) and weighted-average-of-quantiles (WAQ) estimators \citep{Athey2023}) instead lose coverage in S3 and S4. These undercoverage issues are an estimand mismatch, not an efficiency loss. Each of the rank- and quantile-based estimators targets a median, rank, or quantile contrast that equals the mean ATE only under symmetry, and the sparse responders of S3 and the skewed individual effects of S4 break that equality, so the intervals concentrate around the wrong target. Both causal forests stay centered on the mean ATE because MSD confines its robustness to split selection and estimates leaves by the difference in means. Robust partitioning therefore buys outlier resistance without the estimand drift that pure median- or quantile-based estimators incur under skew.
\section{Empirical Applications}
\label{sec:empappl}
We provide two empirical applications that instantiate the skewness contrast of the simulation study. The Progresa vote share \citep{de_la_o_conditional_2013} in the first empirical study in Subsection \ref{sec:progresa} is a strongly right-skewed outcome variable. This case aligns with scenario~S4 of Table \ref{tab:scenarios} in our simulation study and predicts the choice of the splitting rule to reshape the CATE distribution, while the AIPW-aggregated ATE is unaffected. However, rank- and quantile-based estimators drift from the mean ATE due to estimand mismatch. The second application in Subsection \ref{sec:ACTG_175} re-visits a study of antiretroviral treatments in HIV-positive adults \citep{hammer_trial_1996, leqi_median_2022}. The outcome variable, CD4
count, is closer to symmetric and matches the scenarios in cases S1 and S2 of Table \ref{tab:scenarios}, in which all methods are expected to agree regarding ATE estimation results. Appendix \ref{append:emp_appl} provides supplementary material for both empirical applications.
\subsection{Progresa}
\label{sec:progresa}
We examine data from the randomized rollout of Mexico’s conditional cash transfer program, Progresa, which was designed as a social policy intervention and provides an opportunity to study its political consequences, as analyzed in \cite{de_la_o_conditional_2013}. In the original experiment, eligible villages were randomly assigned to begin receiving program benefits either well before (early treatment) or shortly before (delayed control) the 2000 presidential election, thereby generating plausibly exogenous variation in program exposure at the local level. Following the empirical strategy in \cite{de_la_o_conditional_2013}, we aggregate the data to the precinct level (417 observations) and focus on electoral support for the incumbent party in the 2000 election, measured as its vote share among eligible voters as the outcome variable. Figure \ref{fig:outcome_distribution_progresa} in Appendix \ref{append:emp_appl} illustrates the right-skewed nature of the target variable.
We note that the Progresa outcome is a vote share bounded in $[0,100]$, so the leafwise location-shift condition~(i) of Assumption~\ref{ass:HL_shift} cannot hold exactly: a pure shift would move probability mass beyond the support boundary.
By the leaf-estimation argument of Section~\ref{sec:MSD}, this does not bias the MSD point estimates, which are differences in means on the estimation sample, and degrades only the interpretation of the splitting anchor to a robust proxy for $\tau(l;\Pi)$.
The analysis adjusts for pre-treatment socioeconomic and political characteristics, including measures of poverty, historical population size, prior voter turnout, and past party vote shares, along with village fixed effects, to improve precision and account for baseline differences. This design leverages the randomized timing of program implementation to identify the treatment effect of exposure to Progresa on subsequent electoral behavior.
\cite{ghosh_robustness_2026} reanalyze the study in \cite{de_la_o_conditional_2013} to apply its proposed regression-adjusted estimator based on \cite{Rosenbaum1993}. In our analysis, reported in Table~\ref{tab:progresa_estimates}, the difference-in-means estimate of 3.62 is borderline, with a 95\% interval of $[-0.05, 7.30]$ that only just includes zero, and OLS with interaction \citep{lin_agnostic_2013} is the single estimator returning a significant positive effect at 4.21 with interval $[0.24, 8.19]$.
Every estimator that adjusts for the skewed outcome attenuates the effect to between $1.3$ and $2.2$ and finds it statistically insignificant, including Rosenbaum's regression adjustment, the EIF and WAQ estimators, and both causal forests. Rosenbaum's adjustment gives the narrowest interval of these. The two forests agree closely, at $1.47$ for MSE and $1.48$ for MSD, which matches the simulation finding in Figure \ref{fig:ATE_precision} that the splitting rule has little effect on the aggregated ATE estimate.
That a statistically significant positive treatment effect appears only under an estimator that is exposed to the skewed tail and vanishes under every robust adjustment matches the broader reassessment of this experiment. \citet{imai_nonpartisan_2020} reanalyze the Progresa rollout of \citet{de_la_o_conditional_2013} under a bias-corrected setup and find no electoral effect on incumbent support, so our robust estimates should be read as further evidence that median-based splitting does not manufacture an outlier-driven treatment effect.
\begin{table}[htbp]
\centering
\caption{ATE estimates of the effect of early Progresa on PRI support rates with corresponding standard errors, 95\% confidence intervals, and interval widths.}
\label{tab:progresa_estimates}
\begin{tabular}{lcccc}
\hline\hline
Estimator & Estimate & Std.\ Error & 95\% CI & CI Width \\
\hline
Diff.Mean & 3.622 & 1.875 & $[-0.052,\; 7.297]$ & 7.348 \\
Lin.OLS \citep{lin_agnostic_2013} & 4.214 & 2.027 & $[\phantom{-}0.240,\; 8.187]$ & 7.947 \\
R.Adj \citep{ghosh_robustness_2026} & 2.185 & 1.338 & $[-0.439,\; 4.808]$ & 5.246 \\
EIF \citep{Athey2023} & 1.953 & 1.648 & $[-1.276,\; 5.182]$ & 6.459 \\
WAQ \citep{Athey2023} & 1.306 & 1.632 & $[-1.892,\; 4.504]$ & 6.396 \\
$\text{CF (MSE)}$ \citep{wager_estimation_2018} & 1.474 & 1.485 & $[-1.436,\; 4.385]$ & 5.822 \\
$\text{CF (MSD)}$ & 1.482 & 1.487 & $[-1.432,\; 4.396]$ & 5.828 \\
\hline\hline
\end{tabular}
\begin{notes}
ATE estimates for the Progresa application in Section~\ref{sec:progresa}. Estimator definitions are given in Appendix~\ref{appendix:ATE_estimators}.
\end{notes}
\end{table}
Figure~\ref{fig:histogram_cate_progresa} shows where the two causal forests differ in terms of CATE estimates. Both place the bulk of precinct-level CATE estimate distribution near the small positive ATE estimate, but the shapes differ. MSE concentrates its estimates on two sharp modes, while MSD spreads them more evenly and carries more mass into both tails. The discreteness of the MSE estimates reflects mean-based splitting settling on a small number of leaf values, whereas MSD's smoother spread follows from its more conservative partitioning. Neither distribution gives evidence of strong structured heterogeneity, and the agreement in location matches the near-identical ATEs. Only the shape and the tails separate the two, which is where the skewed outcome makes the splitting rule matter.
We further note that the pointwise confidence intervals behind CATE estimates in Figure \ref{fig:histogram_cate_progresa}, reported in Figure~\ref{fig:plot_cate_progresa} of Appendix \ref{append:emp_appl}, are wide and almost all contain the ATE estimate, so the causal forest CATE estimates for the Progresa data provide no evidence of substantial heterogeneity.
\begin{figure}[htbp]
\centering
\caption{Distribution of precinct-level CATE estimates for the Progresa data.}
\includegraphics[width=\textwidth]{histogram_CATE_progresa.pdf}
\label{fig:histogram_cate_progresa}
\begin{notes}
Causal forests with MSD (orange) and MSE (grey) splitting. The dashed, vertical line in black marks the similar ATE estimate of both causal forests.
\end{notes}
\end{figure}
\subsection{ACTG 175}
\label{sec:ACTG_175}
The second empirical study applies our methods to data from a randomized trial studying antiretroviral treatments in HIV-positive adults \citep{hammer_trial_1996, leqi_median_2022, juraska_speff2trial_2022}. Participants were assigned either to zidovudine monotherapy ($D_i=0$) or to alternative antiretroviral therapies ($D_i=1$), and their CD4 (also known as T helper cells) count was recorded approximately 96 weeks post-randomization. As in the Progresa application, condition~(i) of Assumption~\ref{ass:HL_shift}
is not directly verifiable, though the CD4 count is unbounded above and less obviously in conflict with a location-shift model than the bounded vote share.
We adjust for a rich set of baseline characteristics, including immune markers (baseline CD4 and CD8 counts), demographics (age, weight, race, gender), clinical indicators (Karnofsky score, symptomatic status, hemophilia), behavioral variables (homosexual activity, drug use history), and prior treatment exposure (zidovudine and antiretroviral use). Table \ref{tab:desc_stats_actg175} reports descriptive statistics of the considered variables. The randomized design yields a known, constant propensity score of $p(x) = 0.75$ across all subjects. Dropping individuals with missing 96-week measurements leaves an analysis sample of $N=1{,}342$. Dropout in ACTG 175 is substantial and has been linked to baseline characteristics, so this complete-case analysis rests on the outcome being missing at random given the covariates we adjust for, which we take as a maintained assumption rather than a tested one.
The ACTG 175 results for the ATE estimates in Table \ref{tab:ate_antiretroviral} are more uniform compared to the results for the Progresa data in Subsection \ref{sec:progresa}. Every ATE estimator finds a large and significant increase in CD4 count under the alternative antiretroviral therapies, with point estimates ranging from about $46$ for WAQ to $67$ for the causal forests. The two causal forests again agree almost exactly, at $66.60$ for MSE and $66.62$ for MSD. The outcome here is less skewed than the Progresa vote share, and the robust and mean-based estimators move together rather than splitting apart.
Figure~\ref{fig:histogram_cate_actg175} confirms the rather uniform findings for causal forest ATE estimates also at the CATE level. The MSD and MSE distributions are both unimodal, centered near the ATE, and overlap almost entirely, with MSD only slightly more dispersed toward lower values. Individual-level CATEs spread from roughly $50$ to $82$, but the two splitting rules are nearly indistinguishable. Where Progresa shows the rules diverging in shape, ACTG 175 shows them coinciding. Similar to Subsection \ref{sec:progresa}, we note that the pointwise confidence intervals of the CATE estimates, reported in Figure~\ref{fig:plot_cate_actg175} of Appendix \ref{append:emp_appl}, are wide and almost all contain the ATE estimates. At the resolution of the conservative scaled-BLB intervals, we cannot reject a constant treatment effect.
We interpret this as a statement about the precision of our procedure on the outcome variable measuring CD4 at 96 weeks, not as a claim that treatment effects are homogeneous. This result is consistent with \citet{leqi_median_2022}, who also study the ACTG 175 data and likewise find only limited evidence of effect heterogeneity on this outcome.
\begin{table}[htbp]
\centering
\caption{ATE estimates of the effect of antiretroviral treatments in HIV-positive adults with corresponding standard errors, 95\% confidence intervals, and interval widths.}
\label{tab:ate_antiretroviral}
\begin{tabular}{lcccc}
\hline\hline
Estimator & Estimate & Std.\ Error & 95\% CI & CI Width \\
\hline
Diff.Mean & 53.830 & 10.901 & $[32.465,\; 75.195]$ & 42.730 \\
Lin.OLS \citep{lin_agnostic_2013} & 63.485 & 8.798 & $[46.240,\; 80.730]$ & 34.489 \\
R.Adj \citep{ghosh_robustness_2026} & 62.188 & 8.709 & $[45.118,\; 79.258]$ & 34.140 \\
EIF \citep{Athey2023} & 49.146 & 10.749 & $[28.078,\; 70.214]$ & 42.135 \\
WAQ \citep{Athey2023} & 46.052 & 10.996 & $[24.500,\; 67.605]$ & 43.105 \\
$\text{CF (MSE)}$ \citep{wager_estimation_2018} & 66.603 & 10.265 & $[46.484,\; 86.722]$ & 40.238 \\
$\text{CF (MSD)}$ & 66.621 & 10.271 & $[46.489,\; 86.753]$ & 40.264 \\
\hline\hline
\end{tabular}
\begin{notes}
ATE estimates for the ACTG 175 application in Section~\ref{sec:ACTG_175}. Estimator definitions are given in Appendix~\ref{appendix:ATE_estimators}.
\end{notes}
\end{table}
\begin{figure}[htbp]
\centering
\caption{Distribution of CATE estimates for the ACTG 175 data.}
\includegraphics[width=\textwidth]{histogram_CATE_actg175.pdf}
\label{fig:histogram_cate_actg175}
\begin{notes}
Causal forests with MSD (orange) and MSE (grey) splitting. The dashed, vertical line in black marks the similar ATE estimate of both causal forests.
\end{notes}
\end{figure}
\section{Conclusion}
\label{sec:conclusion}
Heavy-tailed and skewed outcomes are common in economic and biomedical
applications and standard CATE estimators might be imprecise under these circumstances \citep{li_robust_2021, ghosh_robustness_2026, Athey2023}.
Instead of relying on the mean-squared-error splitting rule \citep{athey_recursive_2016, wager_estimation_2018}, we have introduced a median-based splitting rule within the causal forest estimator that targets the CATE while resisting heavy-tailed and contaminated outcomes.
Under the leafwise location-shift model, the Hodges--Lehmann estimator \citep{hodges_estimates_1963} is consistent for the leaf effect and median unbiased, and the resulting Median Squared Deviation (MSD) criterion in Section \ref{sec:MSD} preserves the variance-bias structure of the mean-based rule. An implementable cross-product form of the criterion avoids density estimation during tree growth. Two further rules, MAD and LMS, accompany it as robust baselines without an analogous honest MSE-based derivation.
The simulation study in Section \ref{sec:simulations} compares causal forests under the four splitting criteria across the four scenarios of Table \ref{tab:scenarios}.
For the CATE, MSD attains the lowest RMSE and absolute bias across all four designs. Its coverage is closest to nominal in the designs where the mean-based rule undercovers, most clearly under the sparse extreme responders of scenario S3, while it overcovers conservatively in the more well-behaved designs.
For the ATE the splitting choice washes out under AIPW, so MSD and MSE coincide and both remain centered on the true ATE. The rank- and quantile-based competitors instead drift to a different estimand once skew breaks the symmetry under which their targets agree with the mean. Robust partitioning buys resistance to outliers in split selection without the estimand drift that pure median- or quantile-based estimators incur under skew. The two applications in Section \ref{sec:empappl}, chosen for their skewed outcomes, reproduce this pattern. MSD and MSE agree on the ATE estimate, and the wide CATE confidence intervals provide no evidence of heterogeneity distinguishable from a constant effect.
Several limitations qualify our presented results. The implemented MSD criterion is the cross-product surrogate. It estimates the leafwise second moment by the squared difference in means and is exact only to leading order, trading a small bias for the avoidance of exact density estimation. Identifying the Hodges--Lehmann anchor with the leaf CATE rests on the within-leaf location-shift condition, which fails under the bimodal and skewed leaf distributions of S3 and S4 in Table \ref{tab:scenarios}.
Where it fails, the criterion still selects stable splits, but its target degrades to a robust proxy for the leaf effect. Robustness also carries a cost in conservatism, as MSD intervals are systematically wider than the mean-based intervals and overcover in the well-behaved designs.
We further note that our scaled BLB confidence intervals around the causal forest CATE estimates (see Appendix \ref{appendix:BLB}) are rather conservative. An exact calibration against alternative variance estimators remains for future work.
Finally, the Hodges--Lehmann anchor is rebuilt from all within-leaf pairwise differences at each candidate split. Table \ref{tab:runtime} in Appendix \ref{app:add_sim_res} shows that this computational effort makes MSD (and MAD) several times slower than the mean-based MSE, an overhead that a selection algorithm avoiding the pairwise set might remove without altering any result \citep{monahan_algorithm_1984}. Evaluating the density-penalty form and broadening the analysis to observational designs with estimated propensity scores are the natural next steps for further extensions.
\acks{We would like to thank Christoph Hanck for valuable feedback on earlier drafts of this manuscript and we are grateful to seminar participants at the RuhrMetrics Research Seminar and the EuroCIM 2023 for helpful comments and discussions. The authors acknowledge partial financial support from TRR 391 Spatio-temporal Statistics for the Transition of Energy and Transport (520388526) by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) and from the Rhine-Ruhr Center for Scientific Data Literacy (DKZ.2R) by the German Federal Ministry of Education and Research (BMBF). During the preparation of this manuscript, we used Anthropic's and OpenAI's large language models (up to and including Claude Opus 4.8 and GPT-5.5) to assist with proofreading and editing. We reviewed and verified all resulting content and are solely responsible for any remaining errors.
}
\bibliography{paper_bibtex_file}