Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
57,410 characters · 8 sections · 69 citation commands
Flexible machine learning estimation of conditional average treatment effects: a blessing and a curse
\abstract{ Causal inference from observational data requires untestable identification assumptions. If these assumptions apply, machine learning (ML) methods can be used to study complex forms of causal effect heterogeneity. Recently, several ML methods were developed to estimate the conditional average treatment effect (CATE). If the features at hand cannot explain all heterogeneity, the individual treatment effects (ITEs) can seriously deviate from the CATE. In this work, we demonstrate how the distributions of the ITE and the CATE can differ when a causal random forest (CRF) is applied. We extend the CRF to estimate the difference in conditional variance between treated and controls. If the ITE distribution equals the CATE distribution, this estimated difference in variance should be small. If they differ, an additional causal assumption is necessary to quantify the heterogeneity not captured by the CATE distribution. The conditional variance of the ITE can be identified when the individual effect is independent of the outcome under no treatment given the measured features. Then, in the cases where the ITE and CATE distributions differ, the extended CRF can appropriately estimate the variance of the ITE distribution while the CRF fails to do so. }
The increasing availability of (big) observational data has tremendously boosted the field of machine learning (ML) Mooney2018. ML provides us with flexible, non-parametric methods to study the observed outcome $Y$, given features $\boldsymbol{X}$, that may involve a treatment (or exposure) $A$, by statistical inference on the (conditional) distributions of $Y \mid \boldsymbol{X}, A{{=}}a$. Therefore, ML methods are excellent at predicting future observations that arise from the same (factual) distribution Dickerman2020. However, it is essential to realize that these models cannot be automatically used to answer `what if' questions for the treatment $A$, i.e. for counterfactual prediction, as associations found in the data are not necessarily causal Hernan2019b, Dickerman2020, Prosperi2020, VanGeloven2020, Mooney2022, Cui2022, Dickerman2022. Statistical inference of associations is thus only one step in causal inference and, as such, in counterfactual prediction Balzer2021.
The critical step for causal inference is linking the distribution of outcomes in a universe where everyone was treated with $a$, i.e. potential outcomes Hernan2019, $Y^{a} \mid \boldsymbol{X}$ to the distribution of the observed data. When working with observational data, we have to make assumptions for the identification of causal quantities of interest that cannot be verified with the data, so ML is insufficient. Instead, we have to rely on the knowledge of experts Hernan2019. Suppose these identification assumptions can be made so that the distributions of potential and observed outcomes can be linked. In that case, causal estimands (targets) of interest can be connected to estimands of the data-generating distribution and estimated with statistical inference. Next to the validity of the identification assumptions, accurate statistical inference is thus necessary for causal inference. When the assumptions are applicable, but the statistical inference is off, e.g., when using misspecified models for the observed outcomes, the causal inference will also be invalid. The flexibility offered by ML methods can thus improve statistical inference Mooney2022, Blakely2019. More precisely, ML methods can be exploited to learn nuisance parameters of the data generating distribution, such as conditional means and propensity scores, which in turn can be used to estimate the causal estimand, as is, for example, done in targeted maximum likelihood estimation (TMLE) Laan2011, Schuler2017.
The increasing availability of diverse data makes studying effect heterogeneity among individuals more feasible. The field of precision medicine aims to understand this heterogeneity to improve individual treatment decisions Kosorok2019. The average treatment effect (ATE), $\mathbb{E}[Y^{1}-Y^{0}]$, might seriously differ from the individual treatment effect (ITE), i.e. the actual change in outcome caused by the exposure for a particular individual $i$ ($Y_{i}^{1}-Y_{i}^{0}$) Kravitz2004. However, it is well known that an ITE is not identifiable because of the fundamental problem of causal inference Holland1986, i.e. it is impossible to observe the different potential outcomes for one individual jointly. On the other hand, marginalized effects like the ATE and the conditional average treatment effect (CATE) become identifiable in the absence of unmeasured confounding. In randomized experiments, this unconfoundedness assumption holds by design. Treatment effect heterogeneity studies thus focus on the estimation of the individualized CATEs, $\mathbb{E}[Y^{1}-Y^{0} \mid \boldsymbol{X}]$, given measured features $\boldsymbol{X}$, but aggregated over remaining unmeasured features, as a proxy for the ITEs Robertson2020. The functional form of effect modification by different levels of the measured features might be very complex, so ML methods are promising tools for estimating CATEs Bica2021.
In recent years several meta-learning strategies for CATE estimation have been proposed. These strategies decompose the CATE estimation into regression problems that can be solved with any suitable ML method (see Caron2022 for a detailed review). T-learners fit separate models for treated and controls and estimate CATEs as the plug-in difference of the conditional mean estimates (see, e.g. Athey2016 and Powers2018). The performance of T-learners will depend on the levels of sparsity and smoothness of conditional means for treated and controls, as well as the choice of the base learner. T-learners generally fail for subgroups where the treated and control samples differ in size, as illustrated by Kunzel2019. X-learners have been proposed to deal with the difference in sample size by first using a T-learner to predict individual treatment effects that are subsequently used to derive the CATEs for treated and controls separately from which a weighted average is derived Kunzel2019. S-learners include treatment assignment as another covariate next to other features, and the CATE is estimated as the difference of the estimated conditional means for treated and controls (see e.g. Hill2011, Foster2011, Green2012 and Imai2013). Estimation with S-learners might suffer from serious finite-sample bias because they do not involve the CATE directly but focus on the conditional means, a problem also known for ATE estimation Chernozhukov2018. To remedy this issue, Hahn2020 introduced the Bayesian causal forest model that extends the work of Hill2011 by including the CATE as an explicit parameter in the model with its own prior. The R-learner directly identifies the CATE by regressing transformed outcomes on transformed treatment assignment using estimates of nuisance parameters in a first step (as we will elaborate on in Section (ref)) Nie2020. The R-learner is also called `double machine learning' and may give unbiased estimates of the average causal effect for finite samples. At the same time, a one-step approach (S-learner) would still be biased Chernozhukov2018. Similarly, the DR-learner deals with augmented inverse probability weighted transformation Robins1994 of observations after constructing estimates of the propensity score and conditional means in a first step Kennedy2020optimal, Fan2022. The cost of making weaker modelling assumptions using the flexible ML methods is slower convergence rates for the estimators, known as the curse of dimensionality Naimi2022. Therefore, much of the ongoing research is focused on comparing the different methods for CATE estimation to derive whether and when they are optimal (see e.g. Wendling2018, Knaus2020, Kennedy2020optimal, Curth2021 and kennedy2023minimax).
The aim of this work is different, as we want to emphasize the difference between the CATE and the ITE. The CATE is much more personalized than the ATE and, thus, an important step towards precision medicine. However, it concerns us that the CATE is sometimes perceived as equivalent to the ITE (see, e.g. Lu2018). Whether the CATE can appropriately approximate the ITE depends on the remaining variability of causal effects given the considered modifiers, e.g. a CATE $\geq 0$ given $\boldsymbol{X}{=}\boldsymbol{x}$ Talisa2021 does not imply that all ITEs $\geq 0$ for those individuals Hand1992. In this work, we investigate whether we can use a causal random forest (CRF) Athey2019 to estimate the variance of the marginal ITE distribution. More specifically, we investigate the performance of the CRF to estimate $\text{var}(Y^1-Y^0\mid \boldsymbol{X}{=}\boldsymbol{x})$ and $\text{var}(Y^1-Y^0)$ that could be essential characteristics next to the ATE and CATE. To do so, we have simulated data from a causal system based on estimates from a real case study and fit the CRF to estimate individual CATEs and compared the distributions of the (random) conditional expectation, $\mathbb{E}[Y^1-Y^0 \mid \boldsymbol{X}]$, and the ITE, $Y^1-Y^0$.
To open up the field of ITE distribution estimation, we derive what identification assumption should hold, additionally to those necessary for marginal causal inference, to identify other characteristics of the conditional ITE distribution. We show that under a conditional independent effect deviation assumption (ITE independent of $Y^{0}$ given the CATE), the (conditional) variance of the ITE becomes identifiable. This new identifiability assumption can also not be verified with data, but contrary to the unconfoundedness assumption, the assumption might be violated in a randomized experiment. To give an idea of how an assumption on the joint distribution of potential outcomes can evolve the field of treatment effect heterogeneity, we extend the CRF algorithm to estimate the variance of the ITE given the measured features. Suppose one would additionally be willing to assume a Gaussian distribution. In that case, the algorithm outputs a conditional ITE distribution (centred at the CATE) that will be degenerate in the absence of remaining effect heterogeneity.
Section (ref) introduces our notation and presents the identification assumptions necessary for CATE estimation. Furthermore, we describe the methodology behind the CRF algorithm and our reality-based data simulation. Section (ref) presents the results of fitting the CRF to datasets simulated under different settings. In Section (ref), we introduce the new causal assumption such that the conditional variance of the ITE becomes identifiable and extend the CRF to estimate this. Furthermore, we present the results of analyzing the simulated datasets with this new algorithm. Finally, we present some concluding remarks and ideas for future research in Section (ref).
Probability distributions of factual and counterfactual outcomes are defined in the potential outcome framework Neyman1990, Rubin1974. Let $Y_{i}$ and $A_{i}$ represent the (factual) stochastic outcome and the random treatment assignment level of individual $i$. Let $Y_{i}^{a}$ equal the potential outcome under an intervention on the treatment to level $a$ ($Y_{i}^{a}$ is counterfactual when $A_{i} \neq a$). We thus rely on a deterministic potential outcome framework, where each level of treatment corresponds to only one outcome for each individual (but its value typically differs between individuals) Robins1989, VanderWeele2012
We will consider only two treatment levels $\{0,1\}$ with $0$ indicating no treatment. Thus, the individual causal effect of an arbitrary individual $i$ is defined as $Y_{i}^{1}-Y_{i}^{0}$ Hernan2019. When we discuss the random variable describing the heterogeneity of the potential outcomes in the population, we do not subscript the variables. We must make some identification assumptions to relate the distribution of potential outcomes to the distribution of observed outcomes. First of all, it is necessary to have access to a set of measured features $\boldsymbol{X}$ so that the treatment assignment is conditionally independent of the potential outcomes.
This independence is called conditional exchangeability (or unconfoundedness) and implies the absence of unmeasured confounding that cannot be verified with observational data Hernan2019. Then there are no features, other than $\boldsymbol{X}$, that $Y^0$ or $Y^1$ depend on and that differ in distribution between individuals with $A~{=}~1$ and $A~{=}~0$. Since we are interested in causal effect heterogeneity, the set of features $\boldsymbol{X}$ will also contain modifiers $\boldsymbol{X}_{\text{m}}$ (i.e. $\exists \boldsymbol{x}_{1}, \boldsymbol{x}_{2}{:}~\mathbb{E}[Y^1-Y^0 \mid \boldsymbol{X}_{\text{m}}=\boldsymbol{x}_{1}] \neq \mathbb{E}[Y^1-Y^0 \mid \boldsymbol{X}_{\text{m}}=\boldsymbol{x}_{2}]$ VanderWeele2009) next to the confounders that are necessary to obtain the independence. A feature can be only a modifier, only a confounder, or both, all on the additive scale. For a feature $L$ that is only a confounder but not a modifier $\forall l{:}~\mathbb{E}[Y^1-Y^0 \mid L=l, \boldsymbol{X}_{\text{m}}=\boldsymbol{x}] = \mathbb{E}[Y^1-Y^0 \mid \boldsymbol{X}_{\text{m}}=\boldsymbol{x}]$, where $\boldsymbol{X}_{\text{m}}$ represents the feature set $\boldsymbol{X}$ without $L$.
Furthermore, we need to assume that the observed outcome of an individual equals the potential outcome for the assigned treatment, referred to as causal consistency Cole2009.
Causal consistency is also referred to as the stable unit treatment value assumption (SUTVA)Imbens2015. Causal consistency implies that potential outcomes are independent of the treatment levels of other individuals (no interference) and that there are no different versions of the exposure levels. Causal consistency can also not be verified with data.
Finally, the probability of receiving treatment should be bounded away from 0 and 1 for all levels of $\boldsymbol{X}$, referred to as positivity Hernan2019.
Positivity is also known as overlap Imbens2015.
As in Athey2019, by causal consistency, we use the parameterization
where $b_{i}$ is the ITE of individual $i$, so that $Y_{i}^{1}~{=}~Y^{0}_{i}+b_{i}$. The conditional mean of $b_{i}$ given features $\boldsymbol{X}_{i}$ equals the CATE $\tau(\boldsymbol{X}_{i})$, where $\tau(\boldsymbol{x}) = \mathbb{E}[Y^1-Y^0 \mid \boldsymbol{X}=\boldsymbol{x}]$. The ITE can thus be divided into $\tau(\boldsymbol{X}_{i})$ and the individual deviation from the CATE that is referred to as $U_{1i}$. For our purposes, it helps to rewrite Equation (ref) as
where $\theta_{0}(\boldsymbol{x})~{=}~\mathbb{E}[Y_{i}^{0} \mid \boldsymbol{X}_{i}{=}\boldsymbol{x}]$, $N_{Yi}$ represents the deviation of $Y_{i}^{0}$ from $\theta_{0}(\boldsymbol{X}_{i})$, $\tau(\boldsymbol{x})~{=}~\mathbb{E}[b_{i} \mid \boldsymbol{X}_{i}{=}\boldsymbol{x}]$, $\mathbb{E}[N_{Yi} \mid \boldsymbol{X}_{i}{=}\boldsymbol{x}]~{=}~0$ and $\mathbb{E}[U_{1i} \mid \boldsymbol{X}_{i}{=}\boldsymbol{x}]~{=}~0$. In this parameterization, the individual $Y^{0}$ and effect $b$ have been rewritten as the sum of their conditional expectations and zero mean deviations from these expectations. Note that other characteristics (different from the mean) of the $N_{Y} \mid \boldsymbol{X}{=}\boldsymbol{x}$ and $U_{1} \mid \boldsymbol{X}{=}\boldsymbol{x}$ distributions can depend on the value of $\boldsymbol{x}$. Furthermore, $U_{1}$ and $N_{Y}$ can be dependent. The latter can never be studied from data without making additional assumptions due to the fundamental problem of causal inference, i.e. we cannot observe the pair $(Y^{0}, Y^{1})$.
To illustrate how the random conditional expectation $\mathbb{E}[Y^{1}-Y^{0} \mid \boldsymbol{X}]$ and $Y^{1}-Y^{0}$ may differ in distribution, we simulate data based on the Framingham Heart Study (FHS) Mahmood2014. We focus on the heterogeneity in the effect of non-alcoholic fatty liver disease on a clinical precursor to heart failure, the left ventricular filling pressure Chiu2020. The association found in the original work was adjusted for age, sex, smoking, alcohol use, diabetes, systolic blood pressure (SBP), antihypertensive-med use, lipid-lowering med use, total cholesterol, high-density lipoprotein cholesterol, triglycerides and fasting glucose. However, for this illustration, we assume that only sex (male $=0$ and female $=1$) and SBP are confounders. We will simulate the following cause-effect relations
where $X_{\text{sex},i}\sim \text{Ber}(p)$, $X_{\text{SBP},i}\sim \mathcal{N}(0, 1)$, $U_{1i}\sim \mathcal{N}(0,\sigma_{1}^{2})$, $N_{Yi} \sim \mathcal{N}(0, \sigma_{0}^{2})$, $N_{Ai} \sim \text{Uni}[0,1]$, and $U_{1i} \protect\mathpalette{\protect\independenT}{\perp} N_{Yi}$. Moreover, there is no unmeasured confounding, i.e. $N_{Ai} \protect\mathpalette{\protect\independenT}{\perp} N_{Yi}, U_{1i}$ so that $A_{i} \protect\mathpalette{\protect\independenT}{\perp} (Y^{1}_{i}, Y^{0}_{i}) \mid X_{\text{sex},i}, X_{\text{SBP},i}$. By causal consistency, the observed outcome $Y_{i}~{=}~Y_{i}^{A}$ equals $Y_{i}^{1}$ when $A_{i}~{=}~1$ and $Y_{i}^{0}$ when $A_{i}~{=}~0$. The parameter values are obtained by fitting a linear mixed model for the relation of fatty liver disease and the left ventricular filling pressure adjusted for standardized SBP and sex,
to the subset of the FHS participants $(n=2356)$ as used by Chiu2020. The values obtained with PROC LOGISTIC and PROC MIXED in SAS equal $\alpha_{0}=-1.7$, $\alpha_{\text{sex}}=-0.1$, $\alpha_{\text{SBP}}=0.4$ (so that $\mathbb{P}(A~{=}~1 \mid X_{\text{sex}}=1)~{=}~0.15$ and $\mathbb{P}(A~{=}~1 \mid X_{\text{sex}}=0)~{=}~0.16$), $\beta_{0}=5.9$, $\beta_{\text{sex}}=0.8$, $\beta_{\text{SBP}}=0.5$, $\tau_{0}=0.45$, $\tau_{\text{sex}}=0.1$, $\tau_{\text{SBP}}=0.15$, $\sigma_{0}^{2}=1.6^2$ and $\sigma_{1}^2=1.4^2$. The distribution of $Y^{1}-Y^{0}$ is shown in Figure (ref), where $\mathbb{E}[Y^{1}-Y^{0}]~{=}~0.5$, $\sqrt{\text{var}(Y^{1}-Y^{0})}~{=}~1.41$ and $\mathbb{P}(Y^{1}-Y^{0}>0)~{=}~0.64$. Furthermore, the distribution of the conditional expectation $\mathbb{E}[Y^{1}-Y^{0} \mid X_{\text{SBP}}, X_{\text{sex}}]$ is shown in Figure (ref) with a standard deviation equal to $0.16$ and $\mathbb{P}\left( \mathbb{E}[ Y^{1} - Y^{0} \mid X_{\text{SBP}}, X_{\text{sex}} ] > 0 \right) ~{=}~ 1.00$. The conditional expectation distribution seriously differs from that of the ITE due to the unmeasured (remaining) effect heterogeneity $(U_{1})$. For completeness, the distributions of $Y^1$ and $Y^0$ are presented in Figure (ref). Moreover, we simulate $X_{0}$, which is a measured variable associated with the level of the individual modifier $U_{1}$, $(U_{1}, X_{0})^{T} \sim \mathcal{N}\left(\boldsymbol{0},
\right)$. For $\rho>0$, $X_{0}$ is another measured modifier. Varying $\rho$ can thus be used to investigate cases where more of the latent individual effect modification can be explained while preserving the distribution of $Y^{1}-Y^{0}$. All programming codes used for this work can be found online at \url{https://github.com/RAJP93/CATE}.
Since the actual causal effects are not observed, defining an appropriate loss function is not straightforward, so using ML methods to study causal effect heterogeneity is challenging Athey2016. Causal trees Athey2016 and causal forests Wager2018 have been introduced to draw inferences in case of complex causal effect heterogeneity. These papers mainly focus on data from randomized experiments and suggest adding traditional propensity score methods Athey2016 or using a different algorithm less sensitive to the complexity of the treatment effect function Wager2018 for observational studies. The CRF procedure as implemented in the causal_forest function from the R-package grf is an example of a generalized random forest (GRF) Athey2019. This CRF can be used on data from randomized experiments and observational studies without unmeasured confounding. This ML method is a random forest-based variant of the R-learner Nie2020. As mentioned in the introduction, R-learners focus on the relation between normalized outcome and treatment assignment as originally used by Robinson1988
where $m(\boldsymbol{x})=\mathbb{E}[Y\mid \boldsymbol{X}{=}\boldsymbol{x}]$, $e(\boldsymbol{x})=\mathbb{E}[A\mid \boldsymbol{X}{=}\boldsymbol{x}]$ and $\forall \boldsymbol{x}, a{:}~\mathbb{E}\left[\widetilde{{}N}_{Yi}~{\mid}~A_{i}{{=}}a, \boldsymbol{X} _{i}{{=}}\boldsymbol{x}\right]=0$. Associational quantities are presented with a tilde, and causal quantities are without. Using this convention, $\widetilde{{}\tau}(\boldsymbol{x})$ represents the conditional association measure of $A$ and $Y$ given $\boldsymbol{X}{=}\boldsymbol{x}$. This association will equal the CATE for $\boldsymbol{X}{=}\boldsymbol{x}$, $\tau(\boldsymbol{x})$, under certain identification assumptions, as we will show next. Using parameterization (ref) one can derive, for $\boldsymbol{X}_{i}=\boldsymbol{x}$,
where the distribution of $A_{i}, U_{1i}$ and $N_{Yi}$ can depend on the value of $\boldsymbol{x}$. If $\mathbb{E}[U_{1}+N_{Y}~{\mid}~A{{=}}1, \boldsymbol{X}{=}\boldsymbol{x}] \neq \mathbb{E}[N_{Y}~{\mid}~A{=}0, \boldsymbol{X}{=}\boldsymbol{x}]$, then $\widetilde{{}\tau}(\boldsymbol{x}) \neq \tau(\boldsymbol{x})$. This can occur when there is remaining (unmeasured) confounding after adjusting for $\boldsymbol{X}$. However, in absence of unmeasured confounding, {i.e.} $A \protect\mathpalette{\protect\independenT}{\perp} N_{Y}, U_{1}$, for $\boldsymbol{X}_{i}=\boldsymbol{x}$,
where $\forall \boldsymbol{x}{:}~\mathbb{E}[U_{1i}+N_{Yi}~{\mid}~A_{i}{=}1, \boldsymbol{X} _{i}{=}\boldsymbol{x}]{=}\mathbb{E}[N_{Yi}~{\mid}~A_{i}{=}0 , \boldsymbol{X} _{i}{=}\boldsymbol{x}]{=}0$. Then, $\widetilde{{}\tau}(\boldsymbol{x})$ equals $\tau(\boldsymbol{x})$, the CATE for $\boldsymbol{X}{=}\boldsymbol{x}$.
In the absence of unmeasured confounding, the R-learner thus allows us to estimate (or predict) the CATEs, $\tau(\boldsymbol{x})$, from observational data. The GRF implementation of the CRF starts by predicting $m(\boldsymbol{x}_{i})$ and $e(\boldsymbol{x}_{i})$ for each $i$ by fitting two separate regression forests consisting of honest trees (each tree is fitted on a random subsample of half the sample size) Wager2018. The out-of-bag predictions $\hat{{}m}^{-i}(\boldsymbol{x}_{i})$ and $\hat{{}e}^{-i}(\boldsymbol{x}_{i})$ are only based on those trees that did not use individual $i$ for training, and used to create the centered outcomes $\widetilde{{}Y}_{i}~{=}~Y_{i}-\hat{{}m}^{-i}(\boldsymbol{x}_{i})$ and $\widetilde{{}A}_{i}~{=}~A_{i}-\hat{{}e}^{-i}(\boldsymbol{x}_{i})$ for individual $i$. Subsequently, for a new set of features $\boldsymbol{x}$, similarity weights $\alpha_{j}(\boldsymbol{x})$ are produced for each observation in the sample, and $\widetilde{{}\tau}(\boldsymbol{x})$ (that can be a complex function of $\boldsymbol{x}$) is estimated as
see Athey2019b for more details. The $\alpha_{j}(\boldsymbol{x})$ are obtained by first growing a set of $B$ (user-specified, default is $2000$) trees for $\widetilde{Y}$. For each tree, a random subsample $\mathcal{I}$ of the available data is taken (fraction is user-specified, default is $0.5$). The subsample is randomly divided into (by default) equally sized $\mathcal{J}_{1}$ and $\mathcal{J}_{2}$. The honest decision tree, in the sense of Wager2018, is only fitted on $\mathcal{J}_{1}$ and optimizes the heterogeneity in the effect of $\widetilde{{}A}$ on $\widetilde{{}Y}$ between the different nodes using gradient-based approximations of treatment-effect estimates in candidate children notes, see Athey2019. The similarity weights, $\alpha_{bj}$, are first estimated per tree $b$ and are non-zero (and equal) for those elements of $\mathcal{J}_{2}$ that fall in the same leaf as $\boldsymbol{x}$, and are averaged over all trees to obtain $\alpha_{j}$. For individuals from the original dataset, (out-of-bag) predictions are made by averaging the similarity weights only over trees that did not use this particular observation during training. The ATE is estimated using the augmented inverse probability weighting (AIPW) estimator Robins1995 and equals
where $\hat{{}\mu}_{1i}(\boldsymbol{x} _{i})~{=}~ \hat{{}m}^{-i}(\boldsymbol{x} _{i})+(1-\hat{{}e}^{-i}(\boldsymbol{x}_{i})) \hat{{}\widetilde{{}\tau}}(\boldsymbol{x} _{i})$ and $\hat{{}\mu}_{0i}(\boldsymbol{x} _{i})~{=}~ \hat{{}m}^{-i}(\boldsymbol{x} _{i})-\hat{{}e}^{-i}(\boldsymbol{x} _{i}) \hat{{}\widetilde{{}\tau}}(\boldsymbol{x} _{i})$. The AIPW estimator is double robust, i.e. consistent if $ \hat{{}\widetilde{{}\tau}}(\boldsymbol{x})$ or $\hat{{}e}(\boldsymbol{x})$ is a consistent estimator of $\tau(\boldsymbol{x})$ or $e(\boldsymbol{x})$ respectively.
In this work, we fit a CRF to the simulated data as described in Section (ref) to estimate the $(X_{\text{sex}}, X_{\text{SBP}}, X_{0})$-CATE for each individual. We vary the sample size, $n \in \{200, 2000, 20000 \}$, and the correlation between the unmeasured modifier $U_{1}$ and the measured $X_{0}$, $\rho \in \{0, 0.25, 0.5, 0.75, 1 \}$, while fixing $\delta=2$. We use the default settings of the causal_forest function, except for the $n=200$ settings where we set min.node.size{=}1. For each simulation, we compute the empirical standard deviation (SD) and positive effect probability (PEP), $\mathbb{P}(Y^1-Y^0~{>}~0)$, of the estimated CATE distribution to estimate the SD and PEP of the ITE distribution, respectively. As presented in Section (ref), their actual values equal $1.41$, and $0.64$, respectively. The ATE (equal to $0.5$) is estimated with the AIPW estimator using the average_treatment_effect function of the grf package. Furthermore, based on $1000$ bootstrap samples, we estimate $95\%$ confidence intervals (CIs) for all three characteristics. Based on $1000$ simulations, we estimate the bias, mean squared error (MSE) and coverage for the different settings. Finally, we estimate the ITE distribution per simulation with a Gaussian kernel density estimator over the estimated CATEs using the $\texttt{density}$ function in R with the default settings.
The bias, MSE and coverage for the ATE, SD and PEP of the ITE distribution based on the CATE distribution, estimated with the CRF, are presented in Table (ref) for the different settings.
In the absence of features ($X_{0}$) that are associated with the unmeasured modifier $U_{1}$, i.e. when $\rho~{=}~0$, the variability in the ITE is seriously underestimated when using the CATE distribution as a proxy as shown in the first row of Figure (ref). Therefore, the SD and PEP of the distribution of the conditional expectation ($\mathbb{E}[Y^1-Y^0 \mid \boldsymbol{X}]$) are biased estimators of the characteristics of the ITE distribution. The bias is the lowest for a small sample size due to a finite-sample effect for both the SD and PEP. For $n=200$, the coverage of the PEP is not much off. For larger sample sizes, the CATE distribution can be estimated more precisely. Then, the bias increases, and the coverage decreases.
We observe the same trend for $\rho~{=}~0.25$. However, for $\rho\geq 0.50$, the bias is more extensive for small sample sizes. In the latter cases, the finite-sample effect of the CATE distribution estimator no longer compensates for the difference between the ITE and CATE distribution. Nevertheless, the uncertainty in the estimate for small sample sizes still results in higher coverage of the PEP. For $\rho~{=}~0.75$, the CATE distribution becomes a reasonable proxy for the ITE distribution, as seen from the fourth row in Figure (ref). Finally, for $\rho~{=}~1$, $U_{1}$ equals $X_{0}$, and there is thus no unmeasured effect modification. In this case, the flexible ML estimation of the CATE distribution can be used to estimate the ITE distribution and understand the variability in the treatment effect. Indeed, the bias of the SD and the PEP become small, and the coverage approaches the nominal probability. The bias of the SD is still not neglectable, so the coverage of the SD deviates from the nominal probability.
We presented examples in which the distribution of $\mathbb{E}[Y^{1}-Y^{0} \mid \boldsymbol{X}]$ differs from that of $Y^{1}-Y^{0}$. To overcome this issue, we should consider the remaining effect heterogeneity. In this section, we show that the variance of $Y^{1}-Y^{0} \mid \boldsymbol{X}{=}\boldsymbol{x}$ is only identifiable when we are willing to make another causal assumption. Under this assumption, we can extend the CRF algorithm by also estimating the conditional variance of the effect for each individual. For the parameterization in Equation (ref), the variance of $Y^{1}-Y^{0} \mid \boldsymbol{X}{=}\boldsymbol{x}$ equals $\sigma_{1}^{2}(\boldsymbol{x})=\mathbb{E}[(U_{1})^{2} \mid \boldsymbol{X}{=}\boldsymbol{x}]$. As derived in Appendix (ref), in the absence of unmeasured confounding, the Robinson decomposition of the squared observations enables us to estimate $$\Delta(\boldsymbol{x}) = \tau(\boldsymbol{x})^{2} + \sigma_{1}^{2}(\boldsymbol{x}) + 2\tau(\boldsymbol{x})\theta_{0}(\boldsymbol{x}) + 2\mathbb{E}[N_{Y} U_{1} \mid \boldsymbol{X}{=}\boldsymbol{x}] $$ from observational data. Subsequently, via estimation of $\tau(\boldsymbol{x})$ and $\theta_{0}(\boldsymbol{x})$, $$\widetilde{\sigma_{1}}^{2}(\boldsymbol{x})= \sigma_{1}^{2}(\boldsymbol{x}) + 2\mathbb{E}[N_{Y} U_{1} \mid \boldsymbol{X}{=}\boldsymbol{x}]$$ can be estimated by substraction. Since $\mathbb{E}[U_{1}\mid\boldsymbol{X}=\boldsymbol{x}]$ and $\mathbb{E}[N_{Y}\mid\boldsymbol{X}=\boldsymbol{x}]$ equal 0, $\widetilde{\sigma_{1}}^{2}(\boldsymbol{x})$ represents the sum of the conditional variance of the ITE and twice the conditional covariance of $Y^{0}$ and the ITE. However, as a result of the fundamental problem of causal inference, $\mathbb{E}[N_{Y} U_{1} \mid \boldsymbol{X}{=}\boldsymbol{x}]$, the conditional expectation of the product of the deviation of $Y^{0}$ from $\theta(\boldsymbol{x})$ and the deviation of $Y^{1}-Y^{0}$ from $\tau(\boldsymbol{x})$, is not identifiable. So, we cannot estimate $\sigma_{1}^{2}(\boldsymbol{x})$ without an additional (cross-world) assumption.
If we can assume that $U_{1} \protect\mathpalette{\protect\independenT}{\perp} N_{Y} \mid \boldsymbol{X}{=}\boldsymbol{x}$, then $\mathbb{E}[N_{Y} U_{1} \mid \boldsymbol{X}{=}\boldsymbol{x}]~{=}~0$ and $\widetilde{\sigma_{1}}^{2}(\boldsymbol{x})=\sigma_{1}^{2}(\boldsymbol{x})$, so that the variance of the causal effect given $\boldsymbol{X}$ becomes identifiable. The assumption implies conditional independence of the outcome under no treatment and the effect, i.e. the deviation of $Y^0 \mid \boldsymbol{X}{=}\boldsymbol{x}$ from $\mathbb{E}[Y^0 \mid \boldsymbol{X}{=}\boldsymbol{x}]$ is independent of the deviation of $Y^1-Y^0 \mid \boldsymbol{X}{=}\boldsymbol{x}$ from the CATE. Therefore, we refer to this assumption as conditional independent effect deviation.
Assumption (ref) implies that all features that affect both $Y^0$ and $Y^1-Y^0$ should be contained in $\boldsymbol{X}$. As an example, one could think of the effectiveness of medical drugs that depends on the amount of enzyme present for an individual, while the presence of the enzyme itself does not inform on the outcome of interest in the absence of the drug. The antiplatelet medicine Clopidogrel reduces the risk of stroke and myocardial infarction in individuals with acute coronary syndrome, but its effect depends on its conversion to an active metabolite which is accomplished by the cytochrome P450 2C19 (CYP2C19) enzyme Craig2022. For individuals with a CYP2C19 gene mutation, the drug is known to have a reduced antiplatelet effect; the CYP219 gene thus results in effect heterogeneity. However, there is no reason to believe that the phenotype affects platelet aggregation in the absence of the drug. In cases where $Y^0$ is still expected to inform on the value of $Y^1-Y^0$ given the levels of $\boldsymbol{X}$, the identification assumption does not apply. Similar to Assumption (ref), this causal assumption cannot be verified with data as it concerns unmeasured features that affect both $Y^{0}$ and $Y^{1}-Y^{0}$ and should be judged by experts in the field of application. However, in contrast to Assumption (ref), no reason guarantees that Assumption (ref) holds in a randomized experiment.
Also, in the case where sufficient features are measured so that the ITE equals the CATE, Assumption (ref) holds as $\forall i{:}~ U_{1i}~{=}~0$ and thus independent of $Y_{i}^{0}$.
If Assumption (ref) (in addition to (ref), (ref) and (ref)) holds, then $\text{var}(Y^1-Y^0 \mid \boldsymbol{X}{=}\boldsymbol{x})$ can be estimated with an extended CRF as presented in Algorithm (ref).
\\
Algorithm (ref) provides us with an estimate for both the CATE and the conditional variance of the ITE given the measured features. The ATE estimate remains the same as for the original CRF. The SD of the effect in the total population equals $$\sqrt{\mathbb{E}\left[\text{var}(Y^{1}-Y^{0} \mid \boldsymbol{X}) + \mathbb{E}\left[Y^{1}-Y^{0} \mid \boldsymbol{X} \right]^{2}\right] - \mathbb{E}[Y^{1}-Y^{0}]^{2}}$$ and is therefore estimated as $$\sqrt{\max\left\{0,n^{-1} \left(\sum_{i{=}1}^{n} \hat{{}\widetilde{{}\sigma_{1}}}^{2}((\boldsymbol{x}_{i})) + \hat{{}\widetilde{{}\tau}}(\boldsymbol{x}_{i})^2 \right) - \widehat{\text{ATE}}^{2}\right\}}.$$
Only when the conditional ITE distribution can be well approximated with a Gaussian distribution the distribution of $Y^{1}-Y^{0} \mid \boldsymbol{X}{=}\boldsymbol{x}$ is identified by the CATE and the conditional variance. Then, by the law of total probability,
For illustration, we will assume the Gaussianity of the conditional ITE distribution in our example to use the extended CRF to estimate the ITE distribution from the simulated datasets. The PEP is now estimated as $n^{-1}\sum_{i{=}1}^{n} \mathbb{P}(Z_{i}>0)$, where $Z_{i} \sim \mathcal{N}\left(\hat{{}\widetilde{{}\tau}}(\boldsymbol{x}_{i}), \max\left\{0,{\hat{{}\widetilde{\sigma_{1}}}}^{2}(\boldsymbol{x})\right\}\right)$.
The Gaussianity assumption plays a different role than the identification assumptions (ref) to (ref). The focus of this work is on the conditional variance (and the CATE) that is only identifiable under assumptions (ref) to (ref). We resort to the Gaussianity assumption to also estimate the conditional effect distribution. In Section (ref), we will discuss that under violation of the Gaussianity assumption, the SD of the effect can still be appropriately estimated with Algorithm (ref), but the PEP and ITE distribution estimates will be off.
The bias, MSE and coverage for the ATE $(=0.5)$, SD $(=1.41)$ and PEP $(=0.64)$ of the ITE distribution, respectively, using the extended CRF while assuming Gaussian distributed $Y^{1}-Y^{0} \mid \boldsymbol{X}{=}\boldsymbol{x}$ are presented in Table (ref) for the different settings of the simulation study described in Section (ref).
In the case of remaining heterogeneity ($\rho<1$), the bias of the SD estimator using the extended CRF is much lower than using the CRF. For larger sample sizes ($n=2000$ and $n=20000$), the small bias is of opposite sign to the one using the CRF. The bias of the PEP is also seriously decreased. For all settings, the MSE of the extended estimator is smaller for both the SD and PEP. Also, the coverage of the SD and PEP did considerably improve. However, for $n=20000$ and $\rho~{=}~0.50$ or $\rho~{=}~0.75$, the coverage probability of the SD did deviate from the nominal level due to the small bias and the narrow CIs.
The extended CRF still performs well for the $\rho~{=}~1$ case, where all variability in causal effect could be explained with the measured features. Only when $n=20000$ the bias of the SD estimator using the extended CRF is slightly higher than for the original estimator due to an overspecified model. The difference is so small that the MSE is of the same magnitude. In this case, the coverage again deviates from the nominal level and is now slightly lower than the coverage using the traditional CRF.
The pointwise mean (and $95\%$ CI), from $1000$ simulations, of the estimated probability density function of the ITE, is presented in Figure (ref) for the different settings.
ML methods are of great value in understanding effect heterogeneity using CATEs. In this work, we have shown that there might be individual effect modifications that cannot be explained by the features collected. As a result, the individualized CATE can still seriously differ from the ITE. Then, the ITE distribution cannot be identified by the distribution of the (random) conditional expectation alone, and applied researchers must be aware of this possible discrepancy. For example, remaining effect heterogeneity beyond heterogeneity in the CATEs can result in a lack of generalizability since the distribution of unmeasured effect modifiers in other populations might seriously differ from that in the sample Seamans2021.
Studying the remaining effect heterogeneity is challenging as the fundamental problem of causal inference prevents us from learning the joint distribution of potential outcomes. Nevertheless, the conditional second moments of the treated and the controls should be similar under remaining effect homogeneity. As an example, we have extended the CRF algorithm Athey2019 also to estimate the difference in conditional variance between treated and controls. If variances are different, under assumptions (ref), (ref), and (ref), the ITE distribution cannot be explained by the CATEs alone. The increased variance among the treated is due to the ITE's conditional variance and the covariance of the ITE and $Y^{0}$. Therefore, to estimate the (conditional) variance of the ITE, we need to assume how the ITE and $Y^{0}$ are correlated. In the examples presented in this work, $Y^{1}-Y^{0} \protect\mathpalette{\protect\independenT}{\perp} Y^{0} \mid \boldsymbol{X}{=}\boldsymbol{x}$ so that the causal assumption of conditional independent effect deviation applies. Under this assumption, the conditional variance of the ITE can be estimated next to the expected effect for each individual. As a result, in contrast to the CRF, the extended CRF can be used to estimate the ITE distribution's variance unbiasedly. It should be clear that for settings where Assumption (ref) is violated, the estimate of the (conditional) ITE variance based on the extended CRF will be biased as we illustrate with several scenarios as presented in Appendix (ref). The (conditional) variance of the ITE is also identifiable when the conditional independent effect deviation is violated, but the joint distribution of $Y^1-Y^0$ and $Y^0$ is known. To identify the (conditional) variance of the ITE, the dependence structure of $Y^1-Y^0$ and $Y^0$ should be known since this cannot be learned from the data. The independence Assumption (ref) presented in this work is an example of this, but other assumptions on the dependence would also suffice.
In the absence of remaining effect heterogeneity, the estimated individual variance will be small, indicating that the CATE can be used as an appropriate proxy for the ITE. When assuming that the conditional ITE distributions can be approximated with Gaussian distributions, as done in our example, the ITE distribution can also be estimated. Note that when the conditional ITE distributions are not Gaussian, the distribution is not captured by the CATE and conditional variance alone. Then, other distributional properties like the PEP estimate will be off, as we have demonstrated for a scenario with non-Gaussian conditional treatment effects in Appendix (ref). However, the (conditional) ITE variance can still be appropriately estimated in this scenario. In this paper, we have focussed on the identification of the (conditional) variance of the causal effect. However, under Assumption (ref) (and assumptions (ref), (ref) and (ref)), also higher moments of the (conditional) ITE distribution are identifiable. One could derive the Robinson decomposition for $Y^3$ similarly to the derivation in Appendix (ref).
The effect modification relations in the example presented in the main paper were linear. However, as illustrated with the scenario presented in Appendix (ref), the linearity is not needed to estimate the (conditional) ITE variance using the extended CRF. We have presented the extended CRF just as an example, and other ML algorithms can be extended in a similar way to estimate the (conditional) variance of causal effects under Assumptions (ref), (ref), (ref) and (ref). This is only possible when the ML method appropriately estimates the CATEs and is not prone to overfitting Naimi2022, Balzer2022. Furthermore, it is important to realize that the objective function considered is chosen to maximize heterogeneity in CATEs, and an algorithm might thus not account for confounders that are no effect modifiers (on the additive scale). Therefore, one should always discuss whether the distribution of potential outcomes is correctly linked to the observed distribution. In Appendix (ref), we present scenarios with confounders that are no effect modifiers (on the additive scale) to illustrate that for the CRF, the orthogonalization step facilitates this link. Without this step, even the ATE estimate can be biased. For studies where confounding is absent (e.g. a randomized experiment), the orthogonalization would thus not be necessary, as illustrated in the scenario presented in Appendix (ref)
Although the extended CRF algorithm can be used in practice, our main aim was to emphasize that the CATE and ITE distributions can differ and present what assumption is necessary to identify the (conditional) variance of the causal effect. Since assumptions on the conditional dependence of $Y^0$ and the ITE cannot be tested with factual data, it will be challenging for epidemiologists and other experts in the field of an application of interest to review such assumptions. Judgement should be made based on reasoning about the causal pathways involved. It will be impossible for some applications to reason this way, but in others, it might be possible, like in the Clopidogrel example we mentioned in this paper. Reasoning about more examples will be an important topic of future interdisciplinary research. With this paper, we hope to open up the field of (conditional) ITE distribution estimation under assumptions like conditional independent effect deviation. The latter is necessary to fully understand to what degree individualized CATEs are informative at the actual individual level.