<figure><img src="image_1.png" /></figure>

HAL
open science

# How to select predictive models for decision making or causal inference?

Matthieu Doutreligne, Gaël Varoquaux

► To cite this version:

Matthieu Doutreligne, Gaël Varoquaux. How to select predictive models for decision making or causal inference?. 2023. hal-03946902v4

HAL Id: hal-03946902

https://hal.science/hal-03946902v4

Preprint submitted on 8 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Open licence - etalab
# How to select predictive models for decision making or causal inference?

Matthieu Doutreligne, Gaël Varoquaux

► To cite this version:

Matthieu Doutreligne, Gaël Varoquaux. How to select predictive models for decision making or causal inference?. GigaScience, 2025, 14, 10.1093/gigascience/giaf016. hal-03946902v5

HAL Id: hal-03946902

https://hal.science/hal-03946902v5

Submitted on 3 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# How to select predictive models for decision-making or causal inference?

Matthieu Doutreligne¹,²,* and Gaël Varoquaux¹

¹Inria, Soda, Saclay, France

²Mission Data, Haute Autorité de Santé, Saint-Denis, France

*Corresponding author: matt.dout@gmail.com

December 29, 2023

## Abstract

**Objective:** We investigate which procedure selects the predictive model most trustworthy to reason on the effect of an intervention and support decision making.

**Methods:** We study such a variety of model selection procedures in practical settings: finite samples settings and without theoretical assumption of *well-specified* models. Beyond standard cross-validation or internal validation procedures, we also study elaborate causal risks. These build proxies of the causal error using “nuisance” re-weighting to compute it on the observed data. We evaluate whether empirically estimated nuisances, which are necessarily noisy, add noise to model selection. We compare different metrics for causal model selection in an extensive empirical study based on a simulation and three healthcare datasets based on real covariates.

**Results:** Among all metrics, the mean squared error, classically used to evaluate predictive modes, is worse. Re-weighting it with propensity score does not bring much improvements. The *R*-risk, which uses as nuisances a model of mean outcome and propensity scores, leads to the best performances. Nuisance corrections are best estimated with flexible estimators such as a super learner.

**Conclusions:** When predictive models are used to reason on the effect of an intervention, they must be evaluated with different procedures than standard predictive settings; using the *R*-risk from causal inference.

**Keywords** Model Selection, Predictive model, Treatment Effect, G-formula, Machine Learning,
How to select predictive models for decision-making or
causal inference?

Matthieu Doutreligne¹,* and Gaël Varoquaux²

¹Soda, Inria Saclay, France and ²Mission Data, Haute Autorité de Santé, France

*matt.dout@gmail.com; matthieu.doutreligne@inria.fr

Abstract

**Background:** We investigate which procedure selects the predictive model most trustworthy to reason on the effect of an intervention and support decision making.

**Methods:** We study a large variety of model selection procedures in practical settings: finite samples settings and without theoretical assumption of *well-specified* models. Beyond standard cross-validation or internal validation procedures, we also study elaborate causal risks. These build proxies of the causal error using “nuisance” re-weighting to compute it on the observed data. We evaluate whether empirically estimated nuisances, which are necessarily noisy, add noise to model selection. We compare different metrics for causal model selection in an extensive empirical study based on a simulation and three healthcare datasets based on real covariates.

**Results:** Among all metrics, the mean squared error, classically used to evaluate predictive modes, is worse. Re-weighting it with propensity score does not bring much improvements in most cases. On average, the *R*-risk, which uses as nuisances a model of mean outcome and propensity scores, leads to the best performances. Nuisance corrections are best estimated with flexible estimators such as a super learner.

**Conclusions:** When predictive models are used to reason on the effect of an intervention, they must be evaluated with different procedures than standard predictive settings; using the *R*-risk from causal inference.

**Key words:** Model Selection, Predictive model, Treatment Effect, G-computation, Machine Learning

Introduction

Extending prediction to prescription needs causality

Prediction models have long been used in biomedical settings, as with risk score or prognostic models [1, 2]. While these have historically been simple models on simple data, this is changing with progress in machine learning and richer medical data [3, 4]. Health predictions can now integrate medical images [5, 6, 7, 8, 9], patient records [10, 11, 12] or clinical notes [13, 14, 15]. Complex data is difficult to control and model, but these models are validated by verifying the accuracy of the prediction on left-out data [16, 17, 18]. Crucial to the clinical adoption of a model predicting a health outcome is that it “can support decisions about patient care” [19]. Precision medicine is about guiding decisions: *eg* will an individual benefit from an intervention such as surgery [20]? An estimate of

the effect of the treatment can be obtained by contrasting model predictions with and without the treatment, but statistical validity requires causal inference [21, 22, 23].

Indeed, concluding on the effect of a treatment is a difficult causal-inference task, as it can be easily compromised by confounding: spurious associations between treatment allocation and baseline health, e.g. only prescribing a drug to mild cases [24, 25]. Predictive modeling is linked to causal inference theory by the concept of *outcome models* (or *g-computation*, *g-estimation*, *g-formula* [26], *Q-model* [21], conditional mean regression [27]). Medical statistics and epidemiology have mostly used other causal-inference methods, modeling treatment assignment with propensity scores [28, 29, 30, 31]. Outcome modeling brings the benefit of going beyond average effects, estimating individualized or conditional average treatment effects (CATE), central to precision medicine. For this purpose, such methods are also invaluable for
# Introduction

## Extending prediction to prescription needs causality

Prediction models have long been used in medicine, as with risk score or prognostic models<sup>47,82</sup>. While these have historically been simple models on simple data, this is changing with progress in machine learning and richer medical data<sup>9,61</sup>. Health predictions can now integrate medical images<sup>36,93,91,73,50</sup>, patient records<sup>46,17,77</sup> or clinical notes<sup>29,88,79</sup>. Complex data is difficult to control and model, but these models are validated by verifying the accuracy of the prediction on left-out data<sup>3,58,86</sup>. Crucial to the clinical adoption of a model predicting a health outcome is that it “can support decisions about patient care”<sup>90</sup>. Precision medicine is about guiding decisions: *eg* will an individual benefit from an intervention such as surgery<sup>20</sup>? An estimate of the effect of the treatment can be obtained by contrasting model predictions with and without the treatment, but statistical validity requires causal inference<sup>78,80,10</sup>.

Indeed, concluding on the effect of a treatment is a difficult causal-inference task, as it can be easily compromised by confounding: spurious associations between treatment allocation and baseline health, *e.g.* only prescribing a drug to mild cases<sup>26,85</sup>. Predictive modeling bridges to causal inference theory under the name of *outcome models* (or G-computation, G-formula<sup>63</sup>, Q-model<sup>78</sup>, conditional mean regression<sup>89</sup>). Medical statistics and epidemiology have mostly used other causal-inference methods, modeling treatment assignment with propensity scores<sup>66,7,12,23</sup>. Outcome modeling brings the benefit of going beyond average effects, estimating individualized or conditional average treatment effects (CATE), central to precision medicine. For this purpose, such methods are also invaluable on randomized trials<sup>83,41,28</sup>.

Outcome-modeling methods, even when specifically designed for causal inference, are numerous: Bayesian Additive Regression Trees<sup>27</sup>, Targeted Maximum Likelihood Estimation<sup>40,71</sup>, causal boosting<sup>59</sup>, causal multivariate adaptive regression splines<sup>59</sup>, random forests<sup>87,5</sup>, Meta-learners<sup>37</sup>, R-learners<sup>52</sup>, Doubly robust estimation<sup>14</sup>... The wide variety of methods raises the problem of selecting between different estimators based on the data at hand. Indeed, estimates of treatment effects can vary markedly across different predictive models<sup>19,18</sup> (illustration in Appendix A).

Given complex health data, which predictive model is to be most trusted to yield valid causal estimates needed to motivate individual treatment decisions? As no single machine-learning method performs best across all data sets, there is a pressing need for clear guidelines to select outcome models for causal inference.

**Objectives and structure of the paper** We study *model selection procedures* in practical settings, without theoretical assumptions often made in statistical literature such as *infinite data* or *well-specified models* (Appendix B). Asymptotic causal-inference theory recommends complex risks, but a practical question is whether model-selection procedures, that rely on data split, can estimate these risks reliably enough. Indeed, these risks come with more quantities to estimate, which may bring additional variance, leading to worse model selection.

We first illustrate the problem of causal model selection and briefly review prior art. Then we anchor causal model selection in the *potential outcome* framework and details the causal risks and model-selection procedure. We then rewrite the so-called *R*-risk as a reweighted version of mean squared difference between the true and estimated individualized treatment effect. Finally, we conduct a thorough empirical study comparing the different metrics on diverse datasets, using a family of simulations and real health data, going beyond prior work limited to specific simulation settings<sup>70,1</sup> (Appendix B).
randomized trials [32, 33, 34].

Outcome-modeling methods, even when specifically designed for causal inference, are numerous: Bayesian Additive Regression Trees [35], Targeted Maximum Likelihood Estimation [36, 37], causal boosting [38], causal multivariate adaptive regression splines [38], random forests [39, 40], Meta-learners [41], R-learners [42], Doubly robust estimation [43]... The wide variety of methods raises the problem of selecting between different estimators based on the data at hand. Indeed, estimates of treatment effects can vary markedly across different predictive models [44, 45, 46, 47] (illustration in Appendix A.1).

Given complex health data, which predictive model is to be most trusted to yield valid causal estimates needed to motivate individual treatment decisions? As no single machine-learning method performs best across all datasets, there is a pressing need for clear guidelines to select outcome models for causal inference.

*Objectives and structure of the paper.* The intersection between machine learning and causal inference is growing rapidly [48, 49]. We focus on *model selection procedures* in practical settings, without theoretical assumptions often made in statistical literature such as *infinite data* or *well-specified models* (Appendix A.2). Asymptotic causal-inference theory recommends complex risks, but a practical question is whether *model-selection procedures*, that rely on data split, can estimate these risks reliably enough. Indeed, these risks come with more quantities to estimate, which may bring additional variance, leading to worse model selection.

We first illustrate the problem of causal model selection. Then we anchor causal model selection in the *potential outcome* framework and detail the causal risks and model-selection procedure. We then rewrite the so-called R-risk as a reweighted version of mean squared difference between the true and estimated individualized treatment effect. Finally, we conduct a thorough empirical study comparing the different metrics on diverse datasets, using a family of simulations and real health data, going beyond prior work limited to specific simulation settings [50, 51] (Appendix A.2).

## Illustration: the best predictor may not estimate best causal effects

Using a predictor to reason on causal effects relies on contrasting the prediction of the outcome for a given individual with and without the treatment. Given various predictors of the outcome, which one should we use? Standard predictive modeling or machine-learning practice selects the predictor that minimizes the expected error on the outcome [17, 18]. However, this predictor may not be the best model to reason about causal effects of an intervention as Figure 1 illustrates. Consider the probability Y of an undesirable outcome (e.g.death), a binary treatment A ∈ {0, 1}, and a covariate X ∈ ℝ summarizing the patient health status (e.g. the Charlson index [52]). We simulate a treatment beneficial (decreases mortality) for patients with high Charlson scores (bad health status) but with little effect for patients in good condition (low Charlson scores).

Figure 1a shows a random forest predictor with a counter-intuitive behavior: it predicts well on average the outcome (as measured by a regression R² score) but perform poorly to estimate causal quantities: the average treatment effect τ (as visible via the error |τ − ̂τ|) or the conditional average treatment effect (the error E[(τ(x) − ̂τ(x))²], called CATE). On the contrary, Figure 1b shows a linear model with smaller R² score but better causal inference.

The problem is that causal estimation requires controlling an error on both treated and non-treated outcome for the same individual: the observed outcome, and the non-observed *counterfactual* one. The linear model is misspecified – the outcome functions are not linear–, leading to poor R²; but it interpolates better to regions where there are few untreated individuals – high Charlson score– and thus gives better causal estimates. Conversely, the random

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Illustration: a) a random-forest predictor with high performance for standard prediction (high R²) but that yields poor causal estimates (large error between true effect τ and estimated ̂τ), b) a linear predictor with smaller prediction performance leading to better causal estimation.</figcaption></figure>

Selecting the predictor with the smallest error to the individual treatment effect E[(τ(x) − ̂τ(x))²] – the τ-risk, eq. 10 – would lead to the best causal estimates; however computing this error is not feasible: it requires access to unknown quantities: τ(x).

While the random forest fits the data better than the linear model, it gives worse causal inference because its error is inhomogeneous between treated and untreated. The R² score does not capture this inhomogeneity.

forest puts weaker assumptions on the data, thus has higher R² score but is biased by the treated population in the poor-overlap region, leading to bad causal estimates.

This toy example illustrates that the classic minimum mean squared error (MSE) criterion is not suited to choosing a model among candidate estimators for causal inference.

## Methods

### Neyman-Rubin Potential Outcomes framework

We first expose the classic construction of the outcome modeling (or G-computation) estimators of causal effect [53, 21, 24].

*Settings.* The Neyman-Rubin Potential Outcomes framework [54, 55] enables statistical reasoning on causal treatment effects: Given an outcome Y ∈ ℝ (e.g. mortality risk or hospitalization length), function of a binary treatment A ∈ A = {0, 1} (e.g. a medical procedure), and baseline covariates X ∈ X ⊂ ℝd, we observe the factual distribution, O = (Y(A), X, A) ~ D = P(y, x, a). However, we want to model the existence of potential observations (unobserved ie. counterfactual) that correspond to a different treatment. Thus we want quantities on the counterfactual distribution O* = (Y(1), Y(0), X, A) ~ D* = P(y(1), y(0), x, a).

Popular quantities of interest (estimands) are: at the population level, the Average Treatment Effect

$$ \text{ATE} = \underset{\tau = \mathbb{E}[Y(1), Y(0) \sim D^*] [Y(1) - Y(0)]}{\text{def}} $$
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 1: Illustration:** a) a random-forest predictor with high performance for standard prediction (high $R^2$) but that yields poor causal estimates (large error between true effect $\tau$ and estimated $\hat{\tau}$), b) a linear predictor with smaller prediction performance leading to better causal estimation.</figcaption></figure>

Selecting the predictor with the smallest error to the individual treatment effect $E[(\tau(x) - \hat{\tau}(x))^2]$ –the $\tau$-risk, eq. 9 – would lead to the best causal estimates; however computing this error is not feasible: it requires access to unknown quantities: $\tau(x)$.

While the random forest fits the data better than the linear model, it gives worse causal inference because its error is inhomogeneous between treated and untreated. The $R^2$ score does not capture this inhomogeneity.

## Illustration: the best predictor may not estimate best causal effects

Using a predictor to reason on causal effects relies on contrasting the prediction of the outcome for a given individual with and without the treatment. Given various predictors of the outcome, which one should we use? Standard predictive modeling or machine-learning practice selects the predictor that minimizes the expected error on the outcome⁵⁸,⁸⁶. However, this predictor may not be the best model to reason about causal effects of an intervention as Figure 1 illustrates. Consider the probability $Y$ of an undesirable outcome (e.g.death), a binary treatment $A \in \{0, 1\}$, and a covariate $X \in \mathbb{R}$ summarizing the patient health status (e.g.the Charlson index¹³). We simulate a treatment beneficial (decreases mortality) for patients with high Charlson scores (bad health status) but with little effect for patients in good condition (low Charlson scores).

Figure 1a shows a random forest predictor with a counter-intuitive behavior: it predicts well on average the outcome (as measured by a regression $R^2$ score) but performs poorly to estimate causal quantities: the average treatment effect $\tau$ (as visible via the error $|\tau - \hat{\tau}|$) or the conditional average treatment effect (the error $E[(\tau(x) - \hat{\tau}(x))^2]$, called CATE). On the contrary, Figure 1b shows a linear model with smaller $R^2$ score but better causal inference.

The problem is that causal estimation requires controlling an error on both treated and non-treated
outcome for the same individual: the observed outcome, and the non-observed *counterfactual* one. The linear model is misspecified—the outcome functions are not linear—, leading to poor $R^2$; but it interpolates better to regions where there are few untreated individuals—high Charlson score— and thus gives better causal estimates. Conversely, the random forest puts weaker assumptions on the data, thus has higher $R^2$ score but is biased by the treated population in the poor-overlap region, leading to bad causal estimates.

This toy example illustrates that the classic minimum Mean Square Error criterion is not suited to choosing a model among candidate estimators for causal inference.

# Methods

## Neyman-Rubin Potential Outcomes framework

**Settings** The Neyman-Rubin Potential Outcomes framework<sup>49,31</sup> enables statistical reasoning on causal treatment effects: Given an outcome $Y \in \mathbb{R}$ (e.g. mortality risk or hospitalization length), function of a binary treatment $A \in \mathcal{A} = \{0, 1\}$ (e.g. a medical procedure), and baseline covariates $X \in \mathcal{X} \subset \mathbb{R}^d$, we observe the factual distribution, $O = (Y(A), X, A) \sim \mathcal{D} = \mathbb{P}(y, x, a)$. However, we want to model the existence of potential observations (unobserved ie. counterfactual) that correspond to a different treatment. Thus we want quantities on the counterfactual distribution $O^* = (Y(1), Y(0), X, A) \sim \mathcal{D}^* = \mathbb{P}(y(1), y(0), x, a)$.

Popular quantities of interest (estimands) are: at the population level, the Average Treatment Effect

$$ \text{ATE} \qquad \tau \stackrel{\text{def}}{=} \mathbb{E}_{Y(1),Y(0)\sim\mathcal{D}^*}[Y(1) - Y(0)]; $$

at the individual level, to model heterogeneity, the Conditional Average Treatment Effect

$$ \text{CATE} \qquad \tau(x) \stackrel{\text{def}}{=} \mathbb{E}_{Y(1),Y(0)\sim\mathcal{D}^*}[Y(1) - Y(0)|X = x]. $$

**Causal assumptions** A given data needs to meet a few assumptions to enable identifying causal estimands<sup>67</sup>. The usual *strong ignorability* assumptions are (details in C): 1) *unconfoundedness* {$Y(0), Y(1)$} $\perp$ $A|X$, 2) *strong overlap* ie. every patient has a strictly positive probability to receive each treatment, 3) *consistency*, and 4) *generalization*.

**Estimating treatment effects with outcome models – g-computation**<sup>62</sup> Should we know the two expected outcomes for a given $X$, we could compute the difference between them, which gives the causal effect of the treatment. These two expected outcomes can be estimated from observed data: the consistency 3 and unconfoundedness 1 assumptions imply the following equality:

$$ \mathbb{E}_{Y(a)\sim\mathcal{D}^*}[Y(a)|X=x] = \mathbb{E}_{Y\sim\mathcal{D}}[Y|X=x, A=a] \qquad (1) $$

On the left, the expectation is taken on the counterfactual unobserved distribution. On the right, the expectation is taken on the factual observed distribution conditionally on the treatment. For the rest of the paper, the expectations will always be taken on the factual observed distribution $\mathcal{D}$. This identification leads to outcome based estimators (ie. g-computation estimators<sup>78</sup>):

$$ \begin{aligned} \tau &= \mathbb{E}_{Y\sim\mathcal{D}^*}[Y(1) - Y(0)|X=x] \\ &= \mathbb{E}_{Y\sim\mathcal{D}}[Y|A=1] - \mathbb{E}_{Y\sim\mathcal{D}}[Y|A=0] \end{aligned} \qquad (2) $$
<table><caption>**Table 1.** Review of causal risks — The R-risk* is called τ-risk<sub>R</sub> in [50].</caption><thead><tr><th>Risk</th><th>Equation</th><th>Reference</th></tr></thead><tbody><tr><td>τ-risk = MSE(τ(X), τ<sub>f</sub>(X))</td><td>E<sub>X~p(X)</sub>[(τ(X) - ţ<sub>f</sub>(X))<sup>2</sup>]</td><td>Eq. 10 [35]</td></tr><tr><td>μ-risk = MSE(Y, f(X))</td><td>E<sub>(Y,X,A)~D</sub>[(Y - f(X;A))<sup>2</sup>]</td><td>Def. 1 [50]</td></tr><tr><td>μ-risk<sub>IPW</sub></td><td>E<sub>(Y,X,A)~D</sub>[(<sup>A</sup>/<sub>e(X)</sub> + <sup>1-A</sup>/<sub>1-e(X)</sub>)(Y - f(X;A))<sup>2</sup>]</td><td>Def. 2 [58]</td></tr><tr><td>τ-risk<sub>IPW</sub></td><td>E<sub>(Y,X,A)~D</sub>[(Y (<sup>A</sup>/<sub>e(X)</sub> - <sup>1-A</sup>/<sub>1-e(X)</sub>) - ţ<sub>f</sub>(X))<sup>2</sup>]</td><td>Def. 3 [39]</td></tr><tr><td>U-risk*</td><td>E<sub>(Y,X,A)~D</sub>[(<sup>Y-m(X)</sup>/<sub>A-e(X)</sub> - ţ<sub>f</sub>(X))<sup>2</sup>]</td><td>Def. 4 [42]</td></tr><tr><td>R-risk*</td><td>E<sub>(Y,X,A)~D</sub>[(<sup>(Y - m(X))</sup>/<sub>(A - e(X))</sub> - ţ<sub>f</sub>(X))<sup>2</sup>]</td><td>Def. 5 [42]</td></tr></tbody></table>

## Model selection procedure

Causal model selection (eq 9) may involve estimating various quantities from the observed data: the outcome model $f$, its induced risk as introduce in the previous section, and possibly nuisances required by the risk. Given a dataset with $N$ samples, we split out a train and a test sets ($\mathcal{T}$, $\mathcal{S}$). We fit each candidate estimator $f \in \mathcal{F}$ on $\mathcal{T}$. We also fit the nuisance models $(\tilde{e}, \tilde{m})$ on the train set $\mathcal{T}$, setting hyperparameters by a nested cross-validation before fitting the nuisance estimators with these parameters on the full train set. Causal quantities are then computed by applying the fitted candidates estimators $f \in \mathcal{F}$ on the test set $\mathcal{S}$. Finally, we compute the model-selection metrics for each candidate model on the test set. This procedure is described in Algorithm 1 and Figure 2.

**Algorithm 1 Model selection procedure**

Given train and test sets ($\mathcal{T}$, $\mathcal{S}$) ~ $\mathcal{D}$, a candidate estimator $f$, a causal metrics $\ell$:

i. Prefit: Learn estimators for unknown nuisance quantities $(\tilde{e}, \tilde{m})$ on the training set $\mathcal{T}$

ii. Fit: learn $\hat{f}(\cdot, a)$ on $\mathcal{T}$

iii. Model selection: $\forall x \in \mathcal{S}$ predict $\hat{f}(x, 1), \hat{f}(x, 0)$ and evaluate the estimator storing the metric value: $\ell(f, \mathcal{S})$ - possibly function of $\tilde{e}$ and $\tilde{m}$

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Estimation procedure for causal model selection.</figcaption></figure>

## R-risk as reweighted oracle metric

The R-risk can be rewritten as a rebalanced τ-risk.

This rewriting involves reweighted residuals: for each potential outcome, $a \in \{0;1\}$, the variance conditionally on $x$ is [60]:

$$ \sigma_y^2(x; a) \stackrel{\text{def}}{=} \int_y (y - \mu_a(x))^2 p(y \mid x = x; A = a) dy $$

Integrating over the population, we get the Bayes squared error: $\sigma_B^2(a) = \int_X \sigma_y^2(x; a)p(x)dx$ and its propensity weighted version: $\tilde{\sigma}_B^2(a) = \int_X \sigma_y^2(x; a)p(x; a)dx$. In case of a purely deterministic

link between the covariates, the treatment, and the outcome, these residual terms are null.

**Proposition 1 (R-risk as reweighted τ-risk)** *Given an outcome model $f$, its R-risk appears as weighted version of its τ-risk (Proof in A.5):*

$$ R\text{-risk}^*(f) = \int_X e(x)(1 - e(x)) (\tau(x) - \tau_f(x))^2 p(x) dx + \tilde{\sigma}_B^2(1) + \tilde{\sigma}_B^2(0) \quad (11) $$

The R-risk targets the oracle at the cost of an overlap reweighting and the addition of the reweighted Bayes residuals, which are independent of $f$. In good overlap regions the weights $e(x)(1 - e(x))$ are close to $\frac{1}{|\mathcal{F}|}$, hence the R-risk is close to the desired gold-standard τ-risk. For randomized control trials, this weight is constant making the R-risk particularly suited for exploring heterogeneity (Appendix A.5)

## Empirical Study

We evaluate the following causal metrics, oracle and feasible versions, presented in Table 1:

$\tilde{\mu}$-risk<sub>IPW</sub>, $\tilde{R}$-risk*, $\tilde{U}$-risk*, $\tilde{\tau}$-risk<sub>IPW</sub>, $\tilde{\mu}$-risk, $\tilde{\mu}$-risk<sub>IPW</sub>, $\tilde{R}$-risk, $\tilde{U}$-risk, $\tilde{\tau}$-risk<sub>IPW</sub>. We benchmark the metrics in a variety of settings: many different simulated data generation processes and three semi-simulated datasets<fnref n="1" />.

The simulations, designed to evaluate the effect of the overlap parameter, also explore more diverse and noisy covariate distributions. They cover a diversity of causal settings such as different ratio of causal effect to background responses, and functional links between covariates, outcome and treatment.

## Caussim: Extensive simulation settings

**Data Generation.** We use simulated data, on which the ground-truth causal effect is known. Going beyond prior empirical studies of causal model selection [50, 51], we use many generative processes, which is needed to reach general conclusions (Appendix A.7).

We generate the response functions using random bases extension, a common method in biostatistics, e.g. functional regression with splines [61, 62]. By allowing the function to vary at specific knots, we control the complexity of the non-linear outcome models. We use random approximation of Radial Basis Function (RBF) kernels [63] to generate the outcome and treatment functions. RBF use the same process as polynomial splines but replace polynomial by Gaussian kernels. Unlike polynomial, Gaussian kernels have decreasing influences in the input space. This avoids unrealistic divergences of the functions at the ends of the feature space. We generate 1000 datasets based on these functions, with random overlap parameters. Example shown in Figure 13 and details in A.7.

<footnote n="1">1 Scripts for the simulations and the selection procedure are available at https://github.com/soda-inria/caussim.</footnote>
<table><caption>**Table 1:** Review of causal risks — The R-risk\* is called τ-risk<sub>R</sub> in Schuler et al. [70].</caption><thead><tr><th>Risk</th><th>Equation</th><th>Reference</th></tr></thead><tbody><tr><td>τ-risk = MSE(τ(X), τ<sub>f</sub>(X))</td><td>E<sub>X~p(X)</sub>[(τ(X) - τ̂<sub>f</sub>(X))<sup>2</sup>]</td><td>Eq. 9<sup>27</sup></td></tr><tr><td>μ-risk = MSE(Y, f(X))</td><td>E<sub>(Y,X,A)~D</sub>[(Y - f(X; A))<sup>2</sup>]</td><td>Def. 1<sup>70</sup></td></tr><tr><td>μ-risk<sup>*</sup><sub>IPW</sub></td><td>E<sub>(Y,X,A)~D</sub>[(<sup>A</sup>/<sub>e(X)</sub> + <sup>1-A</sup>/<sub>1-e(X)</sub>)(Y - f(X; A))<sup>2</sup>]</td><td>Def. 2<sup>38</sup></td></tr><tr><td>τ-risk<sup>*</sup><sub>IPW</sub></td><td>E<sub>(Y,X,A)~D</sub>[(Y (<sup>A</sup>/<sub>e(X)</sub> - <sup>1-A</sup>/<sub>1-e(X)</sub>) - τ̂<sub>f</sub>(X))<sup>2</sup>]</td><td>Def. 3<sup>87</sup></td></tr><tr><td>U-risk<sup>*</sup></td><td>E<sub>(Y,X,A)~D</sub>[(<sup>Y-m(X)</sup>/<sub>A-e(X)</sub> - τ̂<sub>f</sub>(X))<sup>2</sup>]</td><td>Def. 4<sup>52</sup></td></tr><tr><td>R-risk<sup>*</sup></td><td>E<sub>(Y,X,A)~D</sub>[(<sup>(Y - m(X))</sup>/<sub>(A - e(X))</sub>τ̂<sub>f</sub>(X))<sup>2</sup>]</td><td>Def. 5<sup>52</sup></td></tr></tbody></table>

This equation builds on two quantities: the conditional expectancy of the outcome given the covariates and either treatment or no no treatment, called response function:

Response function

$$ \mu_a(x) \stackrel{\text{def}}{=} \mathbb{E}_{Y \sim \mathcal{D}}[Y | X = x, A = a] $$

Given a sample of data and the oracle response functions $\mu_0$, $\mu_1$, the finite sum version of Equation 2 leads to an estimator of the ATE written:

$$ \hat{\tau} = \frac{1}{n} \left( \sum_{i=1}^{n} \mu_1(x_i) - \mu_0(x_i) \right) \qquad (3) $$

This estimator is an oracle finite sum estimator by opposition to the population expression of $\tau$, $\mathbb{E}[\mu_1(x_i) - \mu_0(x_i)]$, which involves an expectation taken on the full distribution $\mathcal{D}$, which is observable but requires infinite data. For each estimator $\ell$ taking an expectation over $\mathcal{D}$, we use the symbol $\hat{\ell}$ to note its finite sum version.

Similarly to the ATE, at the individual level, the CATE:

$$ \tau(x) = \mu_1(x) - \mu_0(x) \qquad (4) $$

**Robinson decomposition** The R-decomposition of the outcome model<sup>64</sup> plays an important role: introducing two quantities, the conditional mean outcome and the probability to be treated (known as propensity score<sup>66</sup>):

Conditional mean outcome

$$ m(x) \stackrel{\text{def}}{=} \mathbb{E}_{Y \sim \mathcal{D}}[Y | X = x] \qquad (5) $$

Propensity score

$$ e(x) \stackrel{\text{def}}{=} \mathbb{P}[A = 1 | X = x] \qquad (6) $$

the outcome can be written

R-decomposition

$$ y(a) = m(x) + (a - e(x))\tau(x) + \varepsilon(x; a) \\ \text{with } \mathbb{E}[\varepsilon(X; A)|X, A] = 0 \qquad (7) $$

$m$ and $e$ are often called nuisances<sup>14</sup>; they are unknown.
**Family of candidate estimators.** We test model selection across different candidate estimators that approximate imperfectly the data-generating process. To build such estimators, we first use an RBF expansion similar to that used for data generation. We choose two random knots and transform the raw data features with a Gaussian kernel. This step is referred as the featurization. Then, we fit a linear regression on these transformed features. We consider two ways of combining these steps for outcome model; we use common nomenclature [41, 64] to refer to these different meta-learners that differ on how they model, jointly or not, the treated and the non treated:

* S Learner: A single learner for both populations, taking the treatment as a supplementary covariate.

* SftLearner: A single set of basis functions is sampled at random for both populations, leading to a given feature space used to model both the treated and controls, then two separate different regressors are fitted on this shared representation.

* TLearner: Two completely different learners for each population, hence separate feature representations and regressors.

For the regression step, we fit a Ridge regression on the transformed features with 6 different choices of the regularization parameter $\lambda \in [10^{-3}, 10^{-2}, 10^{-1}, 1, 10^{1}, 10^{2}]$, coupled with a TLearner or a SftLearner. We sample 10 different random basis for learning and featurization yielding a family $\mathcal{F}$ of 120 candidate estimators.

## Semi-simulated datasets

**Datasets.** We also use three semi-simulated data adding a known synthetic causal effect to real-non synthetic-healthcare covariate. ACIC 2016 [45] is based on the Collaborative Perinatal Project [65], a RCT studying infants' developmental disorders containing 4,802 individuals and 55 features. We used 770 dataset instances: 10 random seeds for each of the 77 simulated settings for the treatment and outcomes. ACIC 2018 [66] simulated treatment and outcomes for the Linked Births and Infant Deaths Database (LBIDD) [67] with $D = 177$ covariates. We used all 432 datasets of size $N = 5000$. Twins [68] is an augmentation of real data on twin births and mortality rates [69]. There are $N = 11984$ samples, and $D = 50$ covariates for which we simulated 1,000 different treatment allocations. Appendix A.7 gives datasets details.

**Family of candidate estimators.** For these three datasets, the family of candidate estimators are gradient boosting trees for both the response surfaces and the treatment<fnref n="2" /> with S-learner, learning rate in {0.01, 0.1, 1}, and maximum number of leaf nodes in {25, 27, 30, 32, 35, 40} resulting in a family of size 18.

**Nuisance estimators.** Drawing from the TMLE literature that uses combination of flexible machine learning methods [37], we model the nuisances ě (respectivley �) with a meta-learner: a stacked estimator of ridge and boosting classifiers (respectively regressions) (hyperparameter selection in Appendix A.7).

## Measuring overlap between treated and non treated

Good overlap between treated and control population is crucial for causal inference (Assumption 3). We introduce the Normalized Total Variation (NTV), a divergence based on the propensity score summarizing the overlap between both populations (Appendix A.6).

## Results: factors driving good model selection

**The R-risk is the best metric on average.** Figure 3 shows the agreement between the ideal ranking of outcome models given the oracle $\tau$-risk and the different feasible causal metrics. We measure this agreement with relative<fnref n="3" /> Kendall tau $\kappa$ (eq. 20) [70]. Given the importance of overlap in how well metrics approximate the oracle $\tau$-risk, we separate strong and weak overlap.

Among all metrics, the classical mean squared error (ie. factual $\mu$-risk) is worse and reweighting it with propensity score ($\mu$-risk$_{IPW}$) does not bring much improvements. The R-risk, which includes a model of mean outcome and propensity scores, leads to the best performances. Interestingly, the U-risk, which uses the same nuisances, deteriorates in weak overlap, probably due to variance inflation when dividing by extreme propensity scores.

Beyond rankings, the differences in terms of absolute ability to select the best model are large: The R-risk selects a model with a $\tau$-risk only 1% higher than the best possible candidate for strong overlap on Caussim, but selecting with the $\mu$-risk or $\mu$-risk$_{IPW}$ –as per machine-learning practice– leads to 10% excess risk and using $\tau$-risk$_{IPW}$ –as in some causal-inference methods [71, 72]– leads to 100% excess risk (Figure 16). Across datasets, the R-risk consistently decreases the risk compared to the $\mu$-risk: from 0.1% to 1% on ACIC2016, 1% from to 20% on ACIC2018, and 0.05% from to 1% on Twins.

**Model selection is harder for low population overlap.** Model selection for causal inference becomes more and more difficult with increasingly different treated and control populations (Figure 4). The absolute Kendall's coefficient correlation with $\tau$-risk drops from 0.9 (excellent agreement with oracle selection) to 0.6 on both Caussim and ACIC 2018 (15).

**Nuisances can be estimated on the same data as outcome models.** Using the train set $\mathcal{T}$ both to fit the candidate estimator and the nuisance estimates is a form of double dipping which can lead errors in nuisances correlated to that of outcome models [42]. In theory,

<figure><img src="image_4.png" /><figcaption>**Figure 3.** The R-risk is the best metric: Relative Kendall's $\tau$ agreement with $\tau$-risk. Strong and Weak overlap correspond to the first and last tertiles of the overlap distribution measured with Normalized Total Variation eq. 17. A.7 presents the same results by adding semi-oracle risks in Figure 14, measured with absolute Kendall's in Figure 15 and with $\tau$-risk gains in Figure 16. Table 4 gives median and IQR of the relative Kendall.</figcaption></figure>

<footnote n="2">2 Scikit-learn regressor, HistGradientBoosting Regressor, and classifier, Hist-GradientBoostingClassifier.</footnote>

<footnote n="3">3 To remove the variance across datasets (some datasets lead to easier model selection than others), we report values for one metric relative to the mean of all metrics for a given dataset instance: Relative $\kappa(\ell, \tau$-risk) = $\kappa(\ell, \tau$-risk) - mean$_\ell$ ($\kappa(\ell, \tau$-risk))</footnote>
## Model-selection risks, oracle and feasible

**Causal model selection** We formalize model selection for causal estimation. Thanks to the g-formula identification (Equation 1), a given outcome model $f: \mathcal{X} \times \mathcal{A} \to \mathcal{Y}$ –learned from data or built from domain knowledge– induces feasible estimates of the ATE and CATE (eqs 3 and 4), $\hat{\tau}_f$ and $\hat{\tau}_f(x)$. Let $\mathcal{F} = \{f: \mathcal{X} \times \mathcal{A} \to \mathcal{Y}\}$ be a family of such estimators. Our goal is to select the best candidate in this family for the observed dataset $O$ using a risk $\ell$:

$$f_{\ell}^{*} = \underset{f \in \mathcal{F}}{\operatorname{argmin}} \ell(f, O) \qquad (8)$$

We now detail possible risks $\ell$, risks useful for causal model selection, and how to compute them.

**The $\tau$-risk: an oracle error risk** As we would like to target the CATE, the following evaluation risk is natural (also called PEHE$^{69,27}$):

$$\tau\text{-risk}(f) \stackrel{\text{def}}{=} \mathbb{E}_{X \sim p(X)}[(\tau(X) - \hat{\tau}_f(X))^2] \qquad (9)$$

Given observed data from $p(X)$, the expectation is computed with a finite sum, as in eq. 3, to give an estimated value $\hat{\tau}\text{-risk}(f)$. However this risk is not feasible as the oracles $\tau(x)$ are not accessible with the observed data $(Y, X, A) \sim \mathcal{D}$.

**Feasible error risks** Table 1 lists feasible risks (Detailed in Appendix D), based on the prediction error of the outcome model and observable quantities. These observable, called nuisances are e-propensity score, eq 6– and m-conditional mean outcome, eq 5. We give the definitions as semi-oracles, function of the true unknown nuisances, but later instantiate them with estimated nuisances, noted $(\tilde{e}, \tilde{m})$. Semi-oracles risks are superscripted with the * symbol.

## Estimation and model selection procedure

Causal model selection (eq 8) may involve estimating various quantities from the observed data: the outcome model $f$, its induced risk as introduced in the previous section, and possibly nuisances required by the risk. Given a dataset with $N$ samples, we split out a train and a test sets $(\mathcal{T}, \mathcal{S})$. We fit each candidate estimator $f \in \mathcal{F}$ on $\mathcal{T}$. We also fit the nuisance models $(\tilde{e}, \tilde{m})$ on the train set $\mathcal{T}$, setting hyperparameters by a nested cross-validation before fitting the nuisance estimators with these parameters on the full train set. Causal quantities are then computed by applying the fitted candidates estimators $f \in \mathcal{F}$ on the test set $\mathcal{S}$. Finally, we compute the model-selection metrics for each candidate model on the test set. This procedure is described in Algorithm 1 and Figure 2.

**Algorithm 1 Model selection procedure**

Given train and test sets $(\mathcal{T}, \mathcal{S}) \sim \mathcal{D}$, a candidate estimators $f$, a causal metrics $\ell$:

1. Prefit: Learn estimators for unknown nuisance quantities $(\tilde{e}, \tilde{m})$ on the training set $\mathcal{T}$

2. Fit: learn $\hat{f}(\cdot, a)$ on $\mathcal{T}$

3. Model selection: $\forall x \in \mathcal{S}$ predict $(\hat{f}(x, 1), \hat{f}(x, 0))$ and evaluate the estimator storing the metric value: $\ell(f, \mathcal{S})$ – possibly function of $\tilde{e}$ and $\tilde{m}$
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Model selection is harder for low population overlap: Kendall's τ agreement with τ-risk. Strong, medium and Weak overlap are the tertiles of the overlap measured with NTV eq. 17. Supplementary materials presents results for all metrics in Figure 18 in absolute Kendall's and continuous overlap values in Figure 15.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Nuisances can be estimated on the same data as outcome models: Results for the R-risk are similar between the shared nuisances/candidate set and the separated nuisances set procedures. Figure 17 details results for all metrics.</figcaption></figure>

these correlations can bias model selection and, strictly speaking, push to split out a third separated dataset –a “nuisance set”– to fit the nuisance models. The drawback is that it depletes the data available for model estimation and selection. However, Figure 5 shows no substantial difference between a procedure with a separated nuisance set and the simpler shared nuisance-candidate set procedure.

Empirically, the best split is 90%/10%: using 90% of the data to estimate both the nuisances and candidates, then computing the risks on the remaining test set for model selection (experiments in Appendix A.8).

**Stacked models are good overall estimators of nuisances.** Stacked nuisances estimators (boosting and linear) lead to feasible metrics with close performances to the oracles ones: the corresponding estimators recover well-enough the true nuisances. One may wonder if simpler models for the nuisance could be useful, in particular in data-poor settings or when the true models are linear. Figure 6 compares causal model selection estimating nuisances with stacked estimators or linear model. It comprises the Twins data, where the true propensity model is linear, and a downsampled version of this data, to study a situation favorable to linear models. In these settings, stacked and linear estimations of the nuisances performs equivalently. Detailed analysis (Figure 20) confirms that using adaptive models –as built by stacking linear models and gradient-boosted trees– suffices to estimate nuisance.

**R-risk is robust to a wide range of effect ratio values.** Beyond overlap, we study for caussim simulations, the effect on model selection of different causal effect ratio to baseline. We vary the empirical mean absolute difference between the causal effect and the baseline, $\Delta_{\mu} = \frac{1}{N} \sum_{i=1}^{N} \left| \frac{\mu_1(x_i) - \mu_0(x_i)}{\mu_0(x_i) + \mu_1(x_i) - \frac{1}{N} \sum_{j=1}^{N} \mu_0(x_j) + \mu_1(x_j)} \right|$, covering a ratio

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Stacked models are good overall estimators of the nuisances: Results are shown only for the R-risk; Figure 19 details every metrics. For Twins, where the true propensity model is linear, stacked and linear estimations of the nuisances performs equivalently, even for a downsampled version (N=4,794).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** R-risk is robust to a wide range of effect ratio: Kendall's τ agreement with τ-risk. Strong, medium and Weak Causal effect ratio are the tertiles of the absolute ratio causal effect to baseline response, $\Delta_{\mu}: \text{Low} [0.04; 2.86[, \text{Medium} [2.86; 16.65[, \text{High} [16.65; 206.53[. Appendix A.9 details this simulation.</figcaption></figure>

range from 0.04 to 206 (median = 9.1). Appendix A.9 details this setup as well as an alternative measure of effect ratio. Figure 7 shows that for high values of the ratio, R-risk is outperformed by the $\mu$-risk$_{IPW}$ and the $\tau$-risk$_{IPW}$. However, on average, the R-risk is still the better risk.

## Discussion and conclusion

**Nuisance models: more gain than pain.** Predictive models are increasingly used to reason about treatment effects, for instance in precision medicine to drive individualized decision. Our results highlight that they should be selected, validated, and tuned using different procedures and error measures than those classically used to assess prediction. Rather, selecting the best outcome model according to the R-risk (eq. Definition 5) leads to more valid causal estimates on average. Estimating the R-risk requires a more complex procedure than standard cross-validation used e.g.in machine learning: it involves fitting nuisance models necessary for model evaluation. Our results show that these can be learned on the same set of data as the outcome models evaluated. The nuisance models must be well estimated (Figure 6). Our results show that using for nuisance models a flexible stacking-based family of estimator suffices for good model selection. To select propensity score models, we used the Brier score, minimized by the true individual probability. An easy mistake is to use calibration errors popular in machine learning [73, 74, 75, 76] as these select not for the individual posterior
<figure><img src="image_5.png" /><figcaption>Figure 2: Estimation procedure for causal model selection.</figcaption></figure>

## R-risk as reweighted oracle metric

The R-risk can be rewritten as a rebalanced τ-risk.

This rewriting involves reweighted residuals: for each potential outcome, $a \in \{0; 1\}$, the variance conditionally on $x$ is$^{72}$:

$$ \sigma_y^2(x; a) \stackrel{\text{def}}{=} \int_y (y - \mu_a(x))^2 p(y | x = x; A = a) dy $$

Integrating over the population, we get the Bayes squared error: $\sigma_B^2(a) = \int_X \sigma_y^2(x; a)p(x)dx$ and its propensity weighted version: $\tilde{\sigma}_B^2(a) = \int_X \sigma_y^2(x; a)p(x; a) dx$. In case of a purely deterministic link between the covariates, the treatment, and the outcome, these residual terms are null.

**Proposition 1 (R-risk as reweighted τ-risk)** Given an outcome model $f$, its R-risk appears as weighted version of its τ-risk (Proof in E.1):

$$ \begin{aligned} \text{R-risk}^*(f) = & \int_x e(x)(1-e(x))(\tau(x) - \tau_f(x))^2 p(x)dx \\ & + \tilde{\sigma}_B^2(1) + \tilde{\sigma}_B^2(0) \end{aligned} \quad (10) $$

The R-risk targets the oracle at the cost of an overlap re-weighting and the addition of the reweighted Bayes residuals, which are independent of $f$. In good overlap regions the weights $e(x)(1-e(x))$ are close to $\frac{1}{4}$, hence the R-risk is close to the desired gold-standard τ-risk. For randomized control trials, this weight is constant making the R-risk particularly suited for exploring heterogeneity (Appendix E.2)

## Empirical Study

We evaluate the following causal metrics, oracle and feasible versions, presented in Table 1:

$\widehat{\mu-risk}_{IPW}$, $\widehat{R-risk}^*$, $\widehat{U-risk}^*$, $\widehat{\tau-risk}_{IPW}^*$, $\widehat{\mu-risk}$, $\widehat{\mu-risk}_{IPW}$, $\widehat{R-risk}$, $\widehat{U-risk}$, $\widehat{\tau-risk}_{IPW}$. We benchmark the metrics in a variety of settings: many different simulated data generation processes and three semi-simulated datasets$^1$.

## Caussim: Extensive simulation settings

**Data Generation** We use simulated data, on which the ground-truth causal effect is known. Going beyond prior empirical studies of causal model selection$^{70,1}$, we use many generative processes, which

<footnote>$^1$Scripts for the simulations and the selection procedure are available at https://github.com/soda-inria/caussim.</footnote>
probability but for an aggregate error rate [77].

*More R-risk to select models driving decisions.* Increasingly complex prediction models integrating richer medical data have flourished because their predictions can be easily demonstrated and validated on left-out data. But using them to underpin a decision on whether to treat or not requires more careful validation, using a metric accounting for the putative intervention, the *R*-risk. On average, the *R*-risk brings a sizeable benefit to select the most adequate model, even when model development is based on treated and untreated population with little differences, as in RCTs. Our conclusions are that without prior knowledge, the *R*-risk is a good default. However, there is much remaining variation, and the *R*-risk will not be optimal for every situation. We have identified one such specific situation: when the causal effect is large compared to the variation of the baseline effect, the $\mu$-risk$_{IPW}$ performs slightly better.

To facilitate better model selection, we provide Python code⁴. This model-selection procedure puts no constraints on the models used to build predictive models: it opens the door to evaluating a wide range of models, from gradient boosting to convolutional neural network, or language models.

## Availability of source code and requirements

Lists the following:

- Project name: Caussim

- Project home page: https://github.com/soda-inria/caussim

- Operating system(s): Platform independent

- Programming language: Python

- License: BSD 3-Clause License

## Declaration

### Competing interests

No competing interest is declared.

### Author contributions statement

M.D. conceived and conducted the experiments, M.D. and G.V. analyzed the results. M.D. and G.V. wrote and reviewed the manuscript.

### Acknowledgments

We acknowledge fruitful discussions with Bénédicte Colnet.
is needed to reach general conclusions (Appendix G.4).

We generate the response functions using random bases extension, a common method in biostatistics, e.g. functional regression with splines<sup>30,56</sup>. By allowing the function to vary at specific knots, we control the complexity of the non-linear outcome models. We use random approximation of Radial Basis Function (RBF) kernels<sup>60</sup> to generate the outcome and treatment functions. RBF use the same process as polynomial splines but replace polynomial by Gaussian kernels. Unlike polynomial, Gaussian kernels have decreasing influences in the input space. This avoids unrealistic divergences of the functions at the ends of the feature space. We generate 1000 datasets based on these functions, with random overlap parameters. Example shown in Figure 12 and details in G.1.

**Family of candidate estimators** We test model selection across different candidate estimators that approximate imperfectly the data-generating process. To build such estimators, we first use a RBF expansion similar to that used for data generation. We choose two random knots and transform the raw data features with a Gaussian kernel. This step is referred as the featurization. Then, we fit a linear regression on these transformed features. We consider two ways of combining these steps for outcome model; we use common nomenclature<sup>37,74</sup> to refer to these different meta-learners that differ on how they model, jointly or not, the treated and the non treated:

* SLearner: A single learner for both population, taking the treatment as a supplementary covariate.

* SftLearner: A single set of basis functions is sampled at random for both populations, leading to a given feature space used to model both the treated and controls, then two separate different regressors are fitted on this shared representation.

* TLearner: Two completely different learners for each population, hence separate feature representations and regressors.

For the regression step, we fit a Ridge regression on the transformed features with 6 different choices of the regularization parameter $\lambda \in [10^{-3}, 10^{-2}, 10^{-1}, 1, 10^{1}, 10^{2}]$, coupled with a TLearner or a SftLearner. We sample 10 different random basis for learning and featurization yielding a family $\mathcal{F}$ of 120 candidate estimators.

## Semi-simulated datasets

**Datasets** We also use three semi-simulated data adding a known synthetic causal effect to real-non synthetic- healthcare covariate. ACIC 2016<sup>18</sup> is based on the Collaborative Perinatal Project<sup>53</sup>, a RCT studying infants’ developmental disorders containing 4,802 individuals and 55 features. We used 770 dataset instances: 10 random seeds for each of the 77 simulated settings for the treatment and outcomes. ACIC 2018<sup>76</sup> simulated treatment and outcomes for the Linked Births and Infant Deaths Database (LBIDD)<sup>44</sup> with $D = 177$ covariates. We used all 432 datasets of size $N = 5000$. Twins<sup>43</sup> is an augmentation of real data on twin births and mortality rates<sup>2</sup>. There are $N = 11984$ samples, and $D = 50$ covariates for which we simulated 1,000 different treatment allocations. Appendix G.2 gives datasets details.

**Family of candidate estimators** For these three datasets, the family of candidate estimators are gradient boosting trees for both the response surfaces and the treatment<sup>2</sup> with S-learner, learning rate in

<footnote n="2"><sup>2</sup>Scikit-learn regressor, HistGradientBoosting Regressor, and classifier, HistGradientBoostingClassifier.</footnote>
# A Additional files

## A.1 Variability of ATE estimation on ACIC 2016

Figure 8 shows ATE estimations for six different models used in g-computation estimators on the 76 configurations of the ACIC 2016 dataset. Outcome models are fitted on half of the data and inference is done on the other half –ie. train/test with a split ratio of 0.5. For each configuration, and each model, this train test split was repeated ten times, yielding non parametric variance estimates [78]. Figure 8 shows large variations obtained across different outcome estimators on semi-synthetic datasets [45]. Flexible models such as random forests are doing well in most settings except when treated and untreated populations differ noticeably, in which case a linear model (ridge) is to be preferred. However random forests with different hyper-parameters (max depth= 2) yield poor estimates. A simple rule of thumb such as preferring flexible models does not work in general; model selection is needed.

Outcome models are implemented with scikit-learn [79] and the following hyper-parameters:

<table><caption>**Table 2.** Hyper-parameters grid used for ACIC 2016 ATE variability</caption><thead><tr><th>Outcome Model</th><th>Hyper-parameters grid</th></tr></thead><tbody><tr><td>Random Forests</td><td>Max depth: [2, 10]</td></tr><tr><td>Ridge regression without treatment interaction</td><td>Ridge regularization: [0.1]</td></tr><tr><td>Ridge regression with treatment interaction</td><td>Ridge regularization: [0.1]</td></tr></tbody></table>

## A.2 Prior work : model selection for outcome modeling (g-computation)

A natural way to select a predictive model for causal inference would be an error measure between a causal quantity such as the CATE and models' estimate. But such error is not a “feasible” risk: it cannot be computed solely from observed data and requires oracle knowledge.

*Simulation studies of causal model selection.* Using eight simulations setups from [38], where the oracle CATE is known, [50] compare four causal risks, concluding that for CATE estimation the best model-selection risk is the so-called *R-risk* [42] –def. 5, below. Their empirical results are clear for randomized treatment allocation but less convincing for observational settings where both simple MSE –μ-risk(f) def. 1– and reweighted MSE –μ-risk<sub>IPW</sub> def. 2– appear to perform better than *R-risk* on half of the simulations. Another work [51] studied empirically both MSE and reweighted MSE risks on the semi-synthetic ACIC 2016 datasets [45], but did not include the *R-risk*. We complete these prior empirical work by studying a wider variety of data generative processes and varying the influence of overlap, an

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Different outcome models lead to different estimation errors on the Average Treatment Effects, on 77 classic simulations with known true causal effect [45]. The different models are ridge regression and random forests with different hyper-parameters (details A.1). The different configurations are plotted as a function of increasing difference between treated and untreated population –see subsection. There is no systematic best performer; data-driven model selection is important.</figcaption></figure>
{0.01, 0.1, 1}, and maximum number of leaf nodes in {25, 27, 30, 32, 35, 40} resulting in a family of size 18.

**Nuisance estimators** Drawing from the TMLE literature that uses combination of flexible machine learning methods<sup>71</sup>, we model the nuisances  (respectivley ṁ) with a meta-learner: a stacked estimator of ridge and boosting classifiers (respectively regressions) (hyperparameter selection in Appendix G.3).

Measuring overlap between treated and non treated

Good overlap between treated and control population is crucial for causal inference (Assumption 2). We introduce the Normalized Total Variation (NTV), a divergence based on the propensity score summariz-
ing the overlap between both populations (Appendix F).

**Results: factors driving good model selection**

**The R-risk is the best metric** Figure 3 shows the agreement between the ideal ranking of outcome
models given the oracle τ-risk and the different feasible causal metrics. We measure this agreement with
relative<fnref n="3" /> Kendall tau κ (eq. 19)³⁵. Given the importance of overlap in how well metrics approximate the
oracle τ-risk, we separate strong and weak overlap.

Among all metrics, the classical mean squared error (ie. factual μ-risk) is worse and reweighting it
with propensity score (μ-risk<sub>IPW</sub>) does not bring much improvements. The R-risk, which includes a

<figure><img src="image_6.png" /><figcaption>**Figure 3: The R-risk is the best metric:** Relative Kendall's τ agreement with τ-risk. Strong and Weak overlap correspond to the first and last tertiles of the overlap distribution measured with Normalized Total Variation eq. 16. G.4 presents the same results by adding semi-oracle risks in Figure 13, measured with absolute Kendall's in Figure 14 and with τ-risk gains in Figure 15. Table 4 gives median and IQR of the relative Kendall.</figcaption></figure>

<footnote n="3">³To remove the variance across datasets (some datasets lead to easier model selection than others), we report values for one metric relative to the mean of all metrics for a given dataset instance: Relative $\kappa(l, \tau$-risk) = $\kappa(l, \tau$-risk) - mean$_l(\kappa(l, \tau$-risk))$</footnote>
important parameter of the data generation process which makes a given causal metric appropriate [80]. We also study how to best adapt cross-validation procedures to causal metrics which themselves come with models to estimate.

**Theoretical studies of causal model selection.** Several theoretical works have proposed causal model selection procedures that are consistent: select the best model in a family given asymptotically large data. These work rely on introducing a CATE estimator in the testing procedure: matching [81], an IPW estimate [72], a doubly robust estimator [82], or debiasing the error with influence functions [51]. However, for theoretical guarantees to hold, the test-set correction needs to converge to the oracle: it needs to be flexible enough –well-posed– and asymptotic data. From a practical perspective, meeting such requirements implies having a good CATE estimate, thus having solved the original problem of causal model selection.

**Statistical guarantees on causal estimation procedures.** Much work in causal inference has focused on procedures that guarantee asymptotically consistent estimators, such as Targeted Machine Learning Estimation (TMLE) [36, 37] or Double Machine Learning [43]. Here also, theories require asymptotic regimes and models to be well-specified.

By contrast, without assuming that estimators are well specified, there exists an upper bound on the oracle error to the CATE (τ-risk) that involves the error on the outcome and the similarity of the distributions of treated and control patients [83]. However, they use this upper bound for model optimization, and do not give insights on model selection. In addition, for hyperparameter selection, they rely on a plugin estimate of the τ-risk built with counterfactual nearest neighbors, which has been shown ineffective [50]. An interesting direction is taken in [84] where the authors derive convergence rates for orthogonal losses such as the R-loss, or the DR-Loss without assuming well-specification of the model for target parameter.

## A.3 Causal assumptions

We assume the following four assumptions, referred as strong ignorability and necessary to assure identifiability of the causal estimands with observational data [56]:

**Assumption 1 (Consistency)** The observed outcome is the potential outcome of the assigned treatment:

$$Y = AY(1) + (1 - A)Y(0)$$

Here, we assume that the intervention A has been well defined. This assumption focuses on the design of the experiment. It clearly states the link between the observed outcome and the potential outcomes through the intervention [24].

**Assumption 2 (Unconfoundedness)**

$$\{Y(0), Y(1)\} \perp A|X$$

This condition –also called ignorability– is equivalent to the conditional independence on e(X) [28]: $\{Y(0), Y(1)\} \perp A|e(X)$.

**Assumption 3 (Overlap, also known as Positivity))**

$$\eta < e(x) < 1 - \eta \quad \forall x \in \mathcal{X} \text{ and some } \eta > 0$$

The treatment is not perfectly predictable. Or with different words, every patient has a chance to be treated and not to be treated. For a given set of covariates, we need examples of both to recover the ATE.

As noted by [80], the choice of covariates X can be viewed as a trade-off between these two central assumptions. A bigger covariates set generally reinforces the ignorability assumption. In the contrary, overlap can be weakened by large X because of the potential inclusion of instruments: variables only linked to the treatment which could lead to arbitrarily small propensity scores.

**Assumption 4 (Generalization)** The training data on which we build the estimator and the test data on which we make the estimation are drawn from the same distribution $\mathcal{D}^*$, also known as the “no covariate shift” assumption [85].

## A.4 Definitions of feasible risks

**Definition 1 (Factual μ-risk)** [60] This is the usual Mean Squared Error on the target y. It is what is typically meant by “generalization error” in supervised learning:

$$\mu\text{-risk}(f) = \mathbb{E} \left[ (Y - f(X;A))^2 \right]$$

**Definition 2 (μ-risk$^{\star}_{IPW}$)** [58] Let the inverse propensity weighting function w(x, a) = $\frac{a}{e(x)} + \frac{1-a}{1-e(x)}$, we define the semi-oracle Inverse Propensity Weighting risk,

$$\mu\text{-risk}^{\star}_{IPW}(f) = \mathbb{E} \left[ \left( \frac{A}{e(X)} + \frac{1-A}{1-e(X)} \right) (Y - f(X;A))^2 \right]$$
model of mean outcome and propensity scores, leads to the best performances. Interestingly, the U-risk, which uses the same nuisances, deteriorates in weak overlap, probably due to variance inflation when dividing by extreme propensity scores.

Beyond rankings, the differences in terms of absolute ability to select the best model are large: The R-risk selects a model with a τ-risk only 1% higher than the best possible candidate for strong overlap on Caussim, but selecting with the μ-risk or μ-risk<sub>IPW</sub> –as per machine-learning practice– leads to 10% excess risk and using τ-risk<sub>IPW</sub> –as in some causal-inference methods<sup>4,25</sup>–leads to 100% excess risk (Figure 15). Across datasets, the R-risk consistently decreases the risk compared to the μ-risk: from 0.1% to 1% on ACIC2016, 1% from to 20% on ACIC2018, and 0.05% from to 1% on Twins.

**Model selection is harder for low population overlap** Model selection for causal inference becomes more and more difficult with increasingly different treated and control populations (Figure 4). The absolute Kendall’s coefficient correlation with τ-risk drops from 0.9 (excellent agreement with oracle selection) to 0.6 on both Caussim and ACIC 2018 (G.4).

<figure><img src="image_7.png" /><figcaption>**Figure 4: Model selection is harder for low population overlap:** Kendall's τ agreement with τ-risk. Strong, medium and Weak overlap are the tertiles of the overlap measured with NTV eq. 16. G.4 presents results for all metrics in Figure 17 in absolute Kendall's and continuous overlap values in Figure 14.</figcaption></figure>

**Nuisances can be estimated on the same data as outcome models** Using the train set $\mathcal{T}$ both to fit the candidate estimator and the nuisance estimates is a form of double dipping which can lead errors in nuisances correlated to that of outcome models<sup>52</sup>. In theory, these correlations can bias model selection and, strictly speaking, push to split out a third separated data set –a “nuisance set”– to fit the nuisance models. The drawback is that it depletes the data available for model estimation and selection. However, Figure 5 shows no substantial difference between a procedure with a separated nuisance set and the simpler shared nuisance-candidate set procedure.

Empirically, the best split is 90%/10%: using 90% of the data to estimate both the nuisances and candidates, then computing the risks on the remaining test set for model selection (experiments in Appendix H).
**Definition 3 (τ-risk<sub>IPW</sub><sup>*</sup>)** [39] The CATE τ(x) can be estimated with a regression against inverse propensity weighted outcomes [71, 72, 39], the τ-risk<sub>IPW</sub>.

$$ \tau\text{-risk}_{IPW}^{*}(f) = \mathbb{E}\left[\left(Y\frac{A - e(X)}{e(X)(1 - e(X))} - \tau_{f}(X)\right)^{2}\right] $$

**Definition 4 (U-risk<sup>*</sup>)** [41, 42] Based on the Robinson decomposition –eq. 8, the U-learner uses the A – e(X) term in the denominator. The derived risk is:

$$ U\text{-risk}^{*}(f) = \mathbb{E}\left[\left(\frac{Y - m(X)}{A - e(X)} - \tau_{f}(X)\right)^{2}\right] $$

Note that extreme propensity weights in the denominator term might inflate errors in the numerator due to imperfect estimation of the mean outcome m.

**Definition 5 (R-risk<sup>*</sup>)** [42, 50] The R-risk also uses two nuisance m and e:

$$ R\text{-risk}^{*}(f) = \mathbb{E}\left[\left(Y - m(X) - (A - e(X))\tau_{f}(X)\right)^{2}\right] $$

It is also based on the Robinson decomposition –eq. 8.

## A.5 Proofs: Links between feasible and oracle risks

**Reformulation of the R-risk as reweighted τ-risk**

**Proposition 1 (R-risk as reweighted τ-risk)** *Proof 1* We consider the R-decomposition: [57],

$$ y(a) = m(x) + (a - e(x))\tau(x) + \varepsilon(x; a) \qquad (12) $$

Where $\mathbb{E}[\varepsilon(X;A)|X,A] = 0$ We can use it as plug in the R-risk formula:

$$ \begin{aligned} R\text{-risk}(f) &= \int_{Y \times X \times A} [(y - m(x)) - (a - e(x))\tau_f(x)]^2 p(y; x; a) dy dx da \\ &= \int_{Y \times X \times A} [(a - e(x))\tau(x) + \varepsilon(x; a) - (a - e(x))\tau_f(x)]^2 p(y; x; a) dy dx da \\ &= \int_{X \times A} (a - e(x))^2 (\tau(x) - \tau_f(x))^2 p(x; a) dx da \\ &+ 2 \int_{Y \times X \times A} (a - e(x)) (\tau(x) - \tau_f(x)) \int_{Y} \varepsilon(x; a) p(y | x; a) dy p(x; a) dx da \\ &+ \int_{X \times A} \int_{Y} \varepsilon^2(x; a) p(y | x; a) dy p(x; a) dx da \end{aligned} $$

The first term can be decomposed on control and treated populations to force e(x) to appear:

$$ \begin{aligned} & \int_{X} (\tau(x) - \tau_{f}(x))^{2} [e(x)^{2} p(x; 0) + (1 - e(x))^{2} p(x; 1)] dx \\ &= \int_{X} (\tau(x) - \tau_{f}(x))^{2} [e(x)^{2} (1 - e(x)) p(x) + (1 - e(x))^{2} e(x) p(x)] dx \\ &= \int_{X} (\tau(x) - \tau_{f}(x))^{2} (1 - e(x)) e(x) [1 - e(x) + e(x)] p(x) dx \\ &= \int_{X} (\tau(x) - \tau_{f}(x))^{2} (1 - e(x)) e(x) p(x) dx. \end{aligned} $$

The second term is null since, $\mathbb{E}[\varepsilon(x, a)|X,A] = 0$.

The third term corresponds to the modulated residuals : $\tilde{\sigma}_{B}^{2}(0) + \tilde{\sigma}_{B}^{2}(1)$

**Interesting special cases**

**Randomization special case.** If the treatment is randomized as in RCTs, $p(A = 1 | X = x) = p(A = 1) = p_{A}$, thus μ-risk<sub>IPW</sub> takes a simpler form:

$$ \mu\text{-risk}_{IPW} = \mathbb{E}_{(Y,X,A) \sim D}\left[\left(\frac{A}{p_A} + \frac{1-A}{1-p_A}\right)(Y - f(X;A))^2\right] $$

However, we still can have large differences between τ-risk and μ-risk<sub>IPW</sub> coming from heterogeneous errors between populations as shown experimentally in [50] and our results below.

Concerning the R-risk, replacing e(x) by its randomized value p<sub>A</sub> in Proposition 1 yields the oracle τ-risk up to multiplicative and additive constants:

$$ R\text{-risk} = p_A (1 - p_A) \tau\text{-risk} + (1 - p_A) \sigma_B^2(0) + p_A \sigma_B^2(1) $$
<figure><img src="image_8.png" /><figcaption>**Figure 5:** Nuisances can be estimated on the same data as outcome models: Results for the R-risk are similar between the *shared nuisances/candidate set* and the *separated nuisances set* procedures. Figure 16 details results for all metrics.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6:** Stacked models are good overall estimators of the nuisances: Results are shown only for the R-risk; Figure 18 details every metric. For Twins, where the true propensity model is linear, *stacked* and *linear* estimations of the nuisances performs equivalently, even for a downsampled version (N=4,794).</figcaption></figure>

**Stacked models are good overall estimators of nuisances** Stacked nuisances estimators (boosting and linear) lead to feasible metrics with close performances to the oracles ones: the corresponding estimators recover well-enough the true nuisances. One may wonder if simpler models for the nuisance could be useful, in particular in data-poor settings or when the true models are linear. Figure 6 compares causal model selection estimating nuisances with stacked estimators or linear model. It comprises the Twins data, where the true propensity model is linear, and a downsampled version of this data, to study a situation favorable to linear models. In these settings, stacked and linear estimations of the nuisances performs equivalently. Detailed analysis (Figure 19) confirms that using adaptive models –as built by stacking linear models and gradient-boosted trees– suffices to estimate nuisance.

# Discussion and conclusion

**Nuisance models: more gain than pain** Predictive models are increasingly used to reason about treatment effects, for instance in precision medicine to drive individualized decision. Our results highlight that they should be selected, validated, and tuned using different procedures and error measures than those classically used to assess prediction. Rather, selecting the best outcome model according to the *R*-risk (eq. Definition 5) leads to more valid causal estimates. Estimating the *R*-risk requires a more complex procedure than standard cross-validation used e.g.in machine learning: it involves fitting nuisance models necessary for model evaluation. Our results show that these can be learned on the same set of data as the outcome models evaluated. The nuisance models must be well estimated (Figure 6). Our results show that using for nuisance models a flexible stacking-based family of estimator suffices for
Thus, selecting estimators with R-risk* in randomized setting controls the τ-risk. This explains the strong performances of R-risk in randomized setups [50] and is a strong argument to use it to estimate heterogeneity in RCTs.

*Oracle Bayes predictor.* If we have access to the oracle Bayes predictor for the outcome ie. $f(x, a) = \mu(x, a)$, then all risks are equivalent up to the residual variance:

$$ \tau\text{-risk}(\mu) = \mathbb{E}_{X-p(X)}[(\tau(X) - \tau_{\mu}(X))^2] = 0 \qquad (13) $$

$$ \begin{aligned} \mu\text{-risk}(\mu) &= \mathbb{E}_{(Y,X,A)-p(Y;X;A)}[(Y - \mu_A(X))^2] \\ &= \int_{X,A} \epsilon(x, a)^2 p(a|x) p(x) dx da \le \sigma_B^2(0) + \sigma_B^2(1) \end{aligned} \qquad (14) $$

$$ \mu\text{-risk}_{IPW}(\mu) = \sigma_B^2(0) + \sigma_B^2(1) \quad \text{from Lemma ??} \qquad (15) $$

$$ R\text{-risk}(\mu) = \tilde{\sigma}_B^2(0) + \tilde{\sigma}_B^2(1) \le \sigma_B^2(0) + \sigma_B^2(1) \quad \text{from Proposition 1} \qquad (16) $$

Thus, differences between causal risks only matter in finite sample regimes. Universally consistent learners converge to the Bayes risk in asymptotic regimes, making all model selection risks equivalent. In practice however, choices must be made in non-asymptotic regimes.

## A.6 Measuring overlap

*Motivation of the Normalized Total Variation.* Overlap is often assessed by comparing visually population distributions as in Figure 1 or computing standardized difference on each feature [86, 29]. While these methods are useful to decide if positivity holds, they do not yield a single measure. Rather, we compute the divergence between the population covariate distributions $\mathbb{P}(X|A=0)$ and $\mathbb{P}(X|A=1)$ [80, 83]. Computing overlap when working only on samples of the observed distribution, outside of simulation, requires a sophisticated estimator of discrepancy between distributions, as two data points never have the same exact set of features. Maximum Mean Discrepancy [87] is typically used in the context of causal inference [60, 83]. However it needs a kernel, typically Gaussian, to extrapolate across neighboring observations. We prefer avoiding the need to specify such a kernel, as it must be adapted to the data which is tricky with categorical or non-Gaussian features, a common situation for medical data.

For simulated and some semi-simulated data, we have access to the probability of treatment for each data point, which sample both densities in the same data point. Thus, we can directly use distribution discrepancy measures and rely on the Normalized Total Variation (NTV) distance to measure the overlap between the treated and control propensities. This is the empirical measure of the total variation distance [88] between the distributions, $TV(\mathbb{P}(X|A=1), \mathbb{P}(X|A=0))$. As we have both distribution sampled on the same points, we can rewrite it a sole function of the propensity score, a low dimensional score more tractable than the full distribution $\mathbb{P}(X|A)$:

$$ \widehat{NTV}(e, 1-e) = \frac{1}{2N} \sum_{i=1}^{N} \left| \frac{e(x_i)}{p_A} - \frac{1-e(x_i)}{1-p_A} \right| \qquad (17) $$

Formally, we can rewrite NTV as the Total Variation distance between the two population distributions. For a population $O = (Y(A), X, A) \sim D$:

$$ \begin{aligned} \text{NTV}(O) &= \frac{1}{2N} \sum_{i=1}^{N} \left| \frac{e(x_i)}{p_A} - \frac{1-e(x_i)}{1-p_A} \right| \\ &= \frac{1}{2N} \sum_{i=1}^{N} \left| \frac{P(A=1|X=x_i)}{p_A} - \frac{P(A=0|X=x_i)}{1-p_A} \right| \end{aligned} $$

Thus NTV approximates the following quantity in expectation over the data distribution $D$:

$$ \begin{aligned} \text{NTV}(D) &= \int_{\mathcal{X}} \left| \frac{p(A=1|X=x)}{p_A} - \frac{p(A=0|X=x)}{1-p_A} \right| p(x) dx \\ &= \int_{\mathcal{X}} \left| \frac{p(A=1, X=x)}{p_A} - \frac{p(A=0, X=x)}{1-p_A} \right| dx \\ &= \int_{\mathcal{X}} \left| p(X=x|A=1) - p(X=x|A=0) \right| dx \end{aligned} $$

For countable sets, this expression corresponds to the Total Variation distance between treated and control populations covariate distributions: $TV(p_0(x), p_1(x))$.
good model selection. To select propensity score models, we used the Brier score, minimized by the true
individual probability. An easy mistake is to use calibration errors popular in machine learning<sup>57,92,51,45</sup>
as these select not for the individual posterior probability but for an aggregate error rate<sup>55</sup>.

**More R-risk to select models driving decisions** Increasingly complex prediction models integrating richer medical data have flourished because their predictions can be easily demonstrated and validated on left-out data. But using them to underpin a decision on whether to treat or not requires more careful validation, using a metric accounting for the putative intervention, the *R*-risk. The *R*-risk brings a sizeable benefit to select the most adequate model, even when model development is based on treated and untreated population with little differences, as in RCTs. To facilitate better model selection, we provide Python code<sup>4</sup>. This model-selection procedure puts no constraints on the models used to build predictive models: it opens the door to evaluating a wide range of models, from gradient boosting to convolutional neutral, or language models.

Acknowledgments

We acknowledge fruitful discussions with Bénédicte Colnet.

Financial disclosure

None reported.

Conflict of interest

The authors declare no potential conflict of interests.

References

[1] Alaa, Ahmed and Schaar, Mihaela Van Der. "Validating Causal Inference Models via Influence Functions". In: *International Conference on Machine Learning* (May 24, 2019), pp. 191–201.

[2] Almond, Douglas, Chay, Kenneth Y., and Lee, David S. "The Costs of Low Birth Weight". In: *The Quarterly Journal of Economics* 120.3 (2005).

[3] Altman, Douglas G, Vergouwe, Yvonne, Royston, Patrick, and Moons, Karel GM. "Prognosis and prognostic research: validating a prognostic model". In: *Bmj* 338 (2009).

[4] Athey, Susan and Imbens, Guido. "Recursive partitioning for heterogeneous causal effects". In: *Proceedings of the National Academy of Sciences* 113.27 (2016), pp. 7353–7360.

[5] Athey, Susan, Tibshirani, Julie, and Wager, Stefan. "Generalized random forests". In: *Annals of Statistics* 47.2 (Apr. 2019), pp. 1148–1178.

[6] Austin, Peter C. "An Introduction to Propensity Score Methods for Reducing the Effects of Confounding in Observational Studies". In: *Multivariate Behavioral Research* 3 (May 2011), pp. 399–424.

<footnote n="4">⁴https://github.com/soda-inria/causal_model_selection</footnote>
*Measuring overlap without the oracle propensity scores:* For ACIC 2018, or for non-simulated data, the true propensity scores are not known. To measure overlap, we rely on flexible estimations of the Normalized Total Variation, using gradient boosting trees to approximate the propensity score. Empirical arguments for this plug-in approach is given in Figure 9.

*Empirical arguments.* We show empirically that NTV is an appropriate measure of overlap by:

* Comparing the NTV distance with the MMD for Caussim which is gaussian distributed in Figure 11,

* Verifying that setups with penalized overlap from ACIC 2016 have a higher total variation distance than unpenalized setups in Figure 10.

* Verifying that the Inverse Propensity Weights extrema (the inverse of the v overlap constant appearing in the overlap Assumption 3) positively correlates with NTV for Caussim, ACIC 2016 and Twins in Figure 12. Even if the same value of the maximum IPW could lead to different values of NTV, we expect both measures to be correlated: the higher the extrem propensity weights, the higher the NTV.

*Estimating NTV in practice.* Finally, we verify that approximating the NTV distance with a learned plug-in estimates of e(x) is reasonable. We used either a logistic regression or a gradient boosting classifier to learn the propensity models for the three datasets where we have access to the ground truth propensity scores: Caussim, Twins and ACIC 2016. We respectively sampled 1000, 1000 and 770 instances of these datasets with different seeds and overlap settings. We first run a hyperparameter search with cross-validation on the train set, then select the best estimator. We refit on the train set this estimator with or without calibration by cross validation and finally estimate the normalized TV with the obtained model. This training procedure reflects the one described in Algorithm 1 where nuisance models are fitted only on the train set.

The hyper parameters are: learning rate ∈ [1e-3, 1e-2, 1e-1, 1], minimum samples leaf ∈ [2, 10, 50, 100, 200] for boosting and L2 regularization ∈ [1e-3, 1e-2, 1e-1, 1] for logistic regression.

Results in Figure 9 comparing bias to the true normalized Total Variation of each dataset instances versus growing true NTV indicate that calibration of the propensity model is crucial to recover a good approximation of the NTV.

## A.7 Experiments

### Details on the data generation process

We use Gaussian-distributed covariates and random basis expansion based on Radial Basis Function kernels. A random basis of RBF kernel enables modeling non-linear and complex relationships between covariates in a similar way to the well known spline expansion. The estimators of the response function are learned with a linear model on another random basis (which can be seen as a stochastic approximation of the full data kernel [63]). We carefully control the amount of overlap between treated and control populations, a crucial assumption for causal inference. Algorithm 2 describes the generation process for one simulation. Figure 13 illustrates 2D examples of the simulation.

* The raw features for both populations are drawn from a mixture of Gaussians: $P(X) = p_A P(X|A=1) + (1-p_A)P(X|A=0)$ where $P(X|A=a)$ is a rotated Gaussian:

$$P(x|A=a) = W \cdot \mathcal{N} \left( \begin{bmatrix} (1-2a)\theta \\ 0 \end{bmatrix}; \begin{bmatrix} \sigma_0 & 0 \\ 0 & \sigma_1 \end{bmatrix} \right) \quad (18)$$

with θ a parameter controlling overlap (bigger yields poorer overlap), W a random rotation matrix and σ₀² = 2; σ₁² = 5. This generation process allows to analytically compute the oracle propensity scores e(x), to simply control for overlap with the parameter θ, the distance between the two Gaussian main axes and to visualize response surfaces.

* A basis expansion of the raw features increases the problem dimension. Using Radial Basis Function (RBF) Nystroem transformation⁵, we expand the raw features into a transformed space. The basis expansion samples randomly a small number of representers in the raw data. Then, it computes an approximation of the full N-dimensional kernel with these basis components, yielding the transformed features z(x). The number of basis functions -ie. knots-, controls the complexity of the ground-truth response surfaces and treatment. We first use this process to draw the non-treated response surface μ₀ and the causal effect τ. We then draw the observations from a mixture two Gaussians, for the treated and non treated. We vary the separation between the two Gaussians to control the overlap between treated and non-treated populations, an important parameter for causal inference (related to η in Assumption 3). Finally, we generate observed outcomes adding Gaussian noise.

More formally, we generate the basis following the original data distribution, $[b_1..b_D] \sim P(x)$, with D=2 in our simulations. Then, we compute an approximation of the full kernel of the data generation process with these representers: $z(x) = [RBF_γ(x, b_d)]_{d=1..D} \cdot Z^T \in \mathbb{R}^D$ with RBF_γ being the Gaussian kernel K(x, y) = $exp(-γ||x - y||^2)$ and Z the normalization constant of the kernel basis, computed as the root inverse of the basis kernel $Z = [K(b_i, b_j)]_{i,j=1..D}^{-1/2}$

* Functions μ₀, τ are distinct linear functions of the transformed features:

$$\mu_0(x) = [z(x); 1] \cdot \beta_{\mu}^T$$

$$\tau(x) = [z(x); 1] \cdot \beta_{\tau}^T$$

Where β_μ and β_τ are sampled from two normal distributions with mean 0 and unit variances N(0, I_D+1).
[7] Austin, Peter C. and Stuart, Elizabeth A. “Moving towards best practice when using inverse probability of treatment weighting (IPTW) using the propensity score to estimate causal treatment effects in observational studies”. In: *Statistics in Medicine* 34.28 (Dec. 10, 2015), pp. 3661–3679.

[8] Battocchi, Keith, Dillon, Eleanor, Hei, Maggie, Lewis, Greg, Oka, Paul, Oprescu, Miruna, and Syrgkanis, Vasilis. *EconML: A Python Package for ML-Based Heterogeneous Treatment Effects Estimation*. https://github.com/microsoft/EconML.

[9] Beam, Andrew L and Kohane, Isaac S. “Big data and machine learning in health care”. In: *Jama* 319.13 (2018), pp. 1317–1318.

[10] Blakely, Tony, Lynch, John, Simons, Koen, Bentley, Rebecca, and Rose, Sherri. “Reflection on modern methods: when worlds collide–prediction, machine learning and causal inference”. In: *International journal of epidemiology* 49.6 (2020), pp. 2058–2064.

[11] Bouthillier, Xavier, Delaunay, Pierre, Bronzi, Mirko, Trofimov, Assya, Nichyporuk, Brennan, Szeto, Justin, Mohammadi Sepahvand, Nazanin, Raff, Edward, Madan, Kanika, Voleti, Vikram, Ebrahimi Kahou, Samira, Michalski, Vincent, Arbel, Tal, Pal, Chris, Varoquaux, Gael, and Vincent, Pascal. “Accounting for Variance in Machine Learning Benchmarks”. In: *Proceedings of Machine Learning and Systems* 3 (Mar. 15, 2021), pp. 747–769.

[12] Casucci, Sabrina, Lin, Li, Hewner, Sharon, and Nikolaev, Alexander. “Estimating the causal effects of chronic disease combinations on 30-day hospital readmissions based on observational Medicaid data”. In: *Journal of the American Medical Informatics Association* 25.6 (2018), pp. 670–678.

[13] Charlson, Mary E., Pompei, Peter, Ales, Kathy L., and MacKenzie, C.Ronald. “A new method of classifying prognostic comorbidity in longitudinal studies: Development and validation”. In: *Journal of Chronic Diseases* 40.5 (Jan. 1987), pp. 373–383.

[14] Chernozhukov, Victor, Chetverikov, Denis, Demirer, Mert, Duflo, Esther, Hansen, Christian, Newey, Whitney, and Robins, James. “Double/Debiased Machine Learning for Treatment and Structural Parameters”. In: *The Econometrics Journal* (2018), p. 71.

[15] Curth, Alicia, Svensson, David, and Weatherall, James. “Really Doing Great at Estimating CATE? A Critical Look at ML Benchmarking Practices in Treatment Effect Estimation”. In: *Neurips Process 2021* (2021), p. 14.

[16] D'Amour, Alexander, Ding, Peng, Feller, Avi, Lei, Lihua, and Sekhon, Jasjeet. “Overlap in observational studies with high-dimensional covariates”. In: *Journal of Econometrics* 221.2 (2021), pp. 644–654.

[17] Desai, Rishi J, Wang, Shirley V, Vaduganathan, Muthiah, Evers, Thomas, and Schneeweiss, Sebastian. “Comparison of machine learning methods with traditional models for use of administrative claims with electronic medical records to predict heart failure outcomes”. In: *JAMA network open* 3.1 (2020), e1918962–e1918962.

[18] Dorie, Vincent, Hill, Jennifer, Shalit, Uri, Scott, Marc, and Cervone, Dan. “Automated versus do-it-yourself methods for causal inference: Lessons learned from a data analysis competition”. In: *Statistical Science* 34.1 (2019), pp. 43–68.

[19] Fang, Gang, Annis, Izabela E, Elston-Lafata, Jennifer, and Cykert, Samuel. “Applying machine learning to predict real-world individual treatment effects: insights from a virtual patient cohort”. In: *Journal of the American Medical Informatics Association* 26.10 (2019), pp. 977–988.
(a) Uncalibrated classifiers

<figure><img src="image_10.png" /></figure>

(b) Calibrated classifiers

<figure><img src="image_11.png" /><figcaption>**Figure 9.** a) Without calibration, estimation of NTV is not trivial even for boosting models. b) Calibrated classifiers are able to recover the true Normalized Total Variation for all datasets where it is available.</figcaption></figure>
[20] Fontana, Mark Alan, Lyman, Stephen, Sarker, Gourab K, Padgett, Douglas E, and MacLean, Catherine H. "Can machine learning algorithms predict which patients will achieve minimally clinically important differences from total joint arthroplasty?" In: *Clinical orthopaedics and related research* 477.6 (2019), p. 1267.

[21] Gao, Zijun, Hastie, Trevor, and Tibshirani, Robert. "Assessment of heterogeneous treatment effect estimation accuracy via matching". In: *Statistics in Medicine* 17 (2021).

[22] Gretton, Arthur, Borgwardt, Karsten M, Rasch, Malte J, Schölkopf, Bernhard, and Smola, Alexander. "A kernel two-sample test". In: *The Journal of Machine Learning Research* 13.1 (2012), pp. 723-773.

[23] Grose, Elysia, Wilson, Samuel, Barkun, Jeffrey, Bertens, Kimberly, Martel, Guillaume, Balaa, Fady, and Khalil, Jad Abou. "Use of Propensity Score Methodology in Contemporary High-Impact Surgical Literature". In: *Journal of the American College of Surgeons* 230.1 (Jan. 1, 2020), 101-112.e2.

[24] Gruber, Susan and van der Laan, Mark J. "tmle: An R Package for Targeted Maximum Likelihood Estimation". In: *Journal of Statistical Software* 51.13 (2012), pp. 1-35.

[25] Gutierrez, Pierre and Gerardy, Jean-Yves. "Causal Inference and Uplift Modeling A review of the literature". In: *Proceedings of The 3rd International Conference on Predictive Applications and APIs*. Proceedings of Machine Learning Research 67 (2016), p. 14.

[26] Hernán, MA and Robins, JM. *Causal Inference: What If*. CRC Boca Raton, FL, 2020.

[27] Hill, Jennifer L. "Bayesian Nonparametric Modeling for Causal Inference". In: *Journal of Computational and Graphical Statistics* 1 (Jan. 1, 2011), pp. 217-240.

[28] Hoogland, Jeroen, IntHout, Joanna, Belias, Michail, Rovers, Maroeska M, Riley, Richard D, E. Harrell Jr, Frank, Moons, Karel GM, Debray, Thomas PA, and Reitsma, Johannes B. "A tutorial on individualized treatment effect prediction from randomized trials with a binary endpoint". In: *Statistics in medicine* 40.26 (2021), pp. 5961-5981.

[29] Horng, Steven, Sontag, David A, Halpern, Yoni, Jernite, Yacine, Shapiro, Nathan I, and Nathanson, Larry A. "Creating an automated trigger for sepsis clinical decision support at emergency department triage using machine learning". In: *PloS one* 12.4 (2017), e0174708.

[30] Howe, Chanelle J., Cole, Stephen R., Westreich, Daniel J., Greenland, Sander, Napravnik, Sonia, and Eron, Joseph J. "Splines for trend analysis and continuous confounder control". In: *Epidemiology (Cambridge, Mass.)* 22.6 (Nov. 2011), pp. 874-875.

[31] Imbens, Guido W. and Rubin, Donald B. *Causal inference in statistics, social, and biomedical sciences*. Cambridge University Press, 2015.

[32] Ionides, Edward L. "Truncated Importance Sampling". In: *Journal of Computational and Graphical Statistics* 17.2 (June 2008), pp. 295-311.

[33] Jesson, Andrew, Mindermann, Sören, Shalit, Uri, and Gal, Yarin. "Identifying Causal-Effect Inference Failure with Uncertainty-Aware Models". In: *Advances in Neural Information Processing Systems* 33 (Oct. 22, 2020), pp. 11637-11649.

[34] Johansson, Fredrik D, Shalit, Uri, Kallus, Nathan, and Sontag, David. "Generalization bounds and representation learning for estimation of potential outcomes and causal effects". In: *The Journal of Machine Learning Research* 23.1 (2022), pp. 7489-7538.
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Good correlation between overlap measured as normalized Total Variation and Maximum Mean Discrepancy (200 sampled Gaussian datasets)</figcaption></figure>

<figure><img src="image_13.png" /></figure>
[35] Kendall, M. G. "A new measure of rank correlation". In: *Biometrika* 30.1-2 (June 1938), pp. 81-93.

[36] Khojaste-Sarakhsi, M, Haghighi, Seyedhamidreza Shahabi, Ghomi, SMT Fatemi, and Marchiori, Elena. "Deep learning for Alzheimer's disease diagnosis: A survey". In: *Artificial Intelligence in Medicine* (2022), p. 102332.

[37] Künzel, Sören R., Sekhon, Jasjeet S., Bickel, Peter J., and Yu, Bin. "Metalearners for estimating heterogeneous treatment effects using machine learning". In: *Proceedings of the National Academy of Sciences* 116.10 (Mar. 5, 2019), pp. 4156-4165.

[38] Laan, Mark J van der, Laan, MJ, and Robins, JM. *Unified methods for censored longitudinal data and causality*. Springer Science & Business Media, 2003.

[39] Laan, Mark J. van der, Polley, Eric C., and Hubbard, Alan E. "Super Learner". In: *Statistical Applications in Genetics and Molecular Biology* 6 (Sept. 16, 2007).

[40] Laan, Mark J. van der and Rose, Sherri. *Targeted Learning*. Springer Series in Statistics. 2011.

[41] Lamont, Andrea, Lyons, Michael D, Jaki, Thomas, Stuart, Elizabeth, Feaster, Daniel J, Tharmaratnam, Kukatharmini, Oberski, Daniel, Ishwaran, Hemant, Wilson, Dawn K, and Van Horn, M Lee. "Identification of predicted individual treatment effects in randomized clinical trials". In: *Statistical methods in medical research* 27.1 (2018), pp. 142-157.

[42] Loiseau, Nicolas, Trichelair, Paul, He, Maxime, Andreux, Mathieu, Zaslavskiy, Mikhail, Wainrib, Gilles, and Blum, Michael G. B. "External control arm analysis: an evaluation of propensity score approaches, G-computation, and doubly debiased machine learning". In: *BMC Medical Research Methodology* 22 (2022).

[43] Louizos, Christos, Shalit, Uri, Mooij, Joris, Sontag, David, Zemel, Richard, and Welling, Max. "Causal Effect Inference with Deep Latent-Variable Models". In: *Advances in neural information processing systems* (Nov. 6, 2017).

[44] MacDorman, M. F. and Atkinson, J. O. "Infant mortality statistics from the linked birth/infant death data set-1995 period data". eng. In: *Monthly Vital Statistics Report* 46.6 Suppl 2 (Feb. 1998), pp. 1-22.

[45] Minderer, Matthias, Djolonga, Josip, Romijnders, Rob, Hubis, Frances, Zhai, Xiaohua, Houlsby, Neil, Tran, Dustin, and Lucic, Mario. "Revisiting the Calibration of Modern Neural Networks". In: *Advances in Neural Information Processing Systems* 34 (2021), pp. 15682-15694.

[46] Mooney, Stephen J and Pejaver, Vikas. "Big data in public health: terminology, machine learning, and privacy". In: *Annual review of public health* 39 (2018), p. 95.

[47] Moons, Karel GM, Royston, Patrick, Vergouwe, Yvonne, Grobbee, Diederick E, and Altman, Douglas G. "Prognosis and prognostic research: what, why, and how?" In: *Bmj* 338 (2009).

[48] Naimi, Ashley I, Mishler, Alan E, and Kennedy, Edward H. "Challenges in Obtaining Valid Causal Effect Estimates with Machine Learning Algorithms". In: *American Journal of Epidemiology* (2021).

[49] Naimi, Ashley I and Whitcomb, Brian W. "Defining and Identifying Average Treatment Effects". In: *American Journal of Epidemiology* (2023).

[50] Nassif, Ali Bou, Talib, Manar Abu, Nasir, Qassim, Afadar, Yaman, and Elgendy, Omar. "Breast cancer detection using artificial intelligence techniques: A systematic literature review". In: *Artificial Intelligence in Medicine* (2022), p. 102276.
<figure><img src="image_14.png" /><figcaption>Figure 12. Maximal value of Inverse Propensity Weights increases exponentially with the overlap as measure by Normalized Total Variation.</figcaption></figure>

• Adding a Gaussian noise, $\epsilon \sim N(0, \sigma(x; a))$, we construct the potential outcomes: $y(a) = (1 - \omega)\mu_0(x) + a\omega\tau(x) + \epsilon(x, a)$ with $\omega$ the treatment effect relative size with respect to the baseline response.

We generated 1000 instances of this dataset with uniformly random overlap parameters $\theta \in [0, 2.5]$.

### Details on the semi-simulated datasets

ACIC 2018 [45]: The initial intervention was a child's birth weight ($A = 1$ if weight < 2.5kg), and outcome was the child's IQ after a follow-up period. The study contained $N = 4802$ data points with $D = 55$ features (5 binary, 27 count data, and 23 continuous). They simulated 77 different setups varying parameters for treatment and response models, overlap, and interactions between treatment and covariates<fnref n="6" />. We used 10 different seeds for each setup, totaling 770 dataset instances.

ACIC 2018 [66]: Starting from data from the Linked Births and Infant Deaths Database (LBIDD) [67] with $D = 177$ covariates, treatment and outcome models are simulated with complex models to reflect different scenarios. The data do not provide the true propensity scores, so we evaluate only feasible metrics, which do not require this nuisance parameter. We used all 432 datasets<fnref n="7" /> of size $N = 5000$.

Twins [68]: It is an augmentation of real data on twin births and mortality rates [69]. There are $N = 11984$ samples (pairs of twins), and $D = 50$ covariates<fnref n="8" />. The outcome is the mortality and the treatment is the weight of the heavier twin at birth. This is a "true" counterfactual dataset [89] in the sense that we have both potential outcomes with each twin. They simulate the treatment with a sigmoid model based on GESTAT10 (number of gestation weeks before birth) and $x$ the 45 other covariates:

$$ 
\begin{gathered}
t_i | \mathbf{x}_i, z_i \sim \text{Bern} \left( \sigma \left( w_0^\top \mathbf{x} + w_h(\mathbf{z}/10 - 0.1) \right) \right) \quad (19) \\
\text{with } w_0 \sim N(0, 0.1 \cdot I), w_h \sim N(5, 0.1)
\end{gathered}
$$

<footnote n="6">6 Original R code available at https://github.com/vdorie/aciccomp/tree/master/2016 to generate 77 simulations settings.</footnote>

<footnote n="7">7 Using the scaling part of the data, from github.com/IBM-HRL-MLHLS/IBM-Causal-Inference-Benchmarking-Framework</footnote>

<footnote n="8">8 We obtained the dataset from https://github.com/AMLab-Amsterdam/CEVAE/tree/master/datasets/TWINS</footnote>
[51] Niculescu-Mizil, Alexandru and Caruana, Rich. "Predicting good probabilities with supervised learning". en. In: *Proceedings of the 22nd international conference on Machine learning - ICML '05*. ACM Press, 2005, pp. 625–632.

[52] Nie, Xinkun and Wager, Stefan. "Quasi-Oracle Estimation of Heterogeneous Treatment Effects". In: *Biometrika* 108.2 (Dec. 13, 2017), pp. 299–319.

[53] Niswander, Kenneth R. and Stroke, United States National Institute of Neurological Diseases and. *The Women and Their Pregnancies: The Collaborative Perinatal Study of the National Institute of Neurological Diseases and Stroke*. Google-Books-ID: A0bdVhlhDQkC. National Institute of Health, 1972. 562 pp.

[54] Pedregosa, Fabian, Varoquaux, Gaël, Gramfort, Alexandre, Michel, Vincent, Thirion, Bertrand, Grisel, Olivier, Blondel, Mathieu, Prettenhofer, Peter, Weiss, Ron, Dubourg, Vincent, Vanderplas, Jake, Passos, Alexandre, Cournapeau, David, Brucher, Matthieu, Perrot, Matthieu, and Duchesnay, Édouard. "Scikit-learn: Machine Learning in Python". In: *Journal of Machine Learning Research* 12.85 (2011), pp. 2825–2830.

[55] Perez-Lebel, Alexandre, Morvan, Marine Le, and Varoquaux, Gaël. "Beyond calibration: estimating the grouping loss of modern neural networks". In: *arXiv preprint arXiv:2210.16315* (2022).

[56] Perperoglou, Aris, Sauerbrei, Willi, Abrahamowicz, Michal, and Schmid, Matthias. "A review of spline function procedures in R". In: *BMC Medical Research Methodology* 19.1 (Mar. 6, 2019), p. 46.

[57] Platt, John C. and Platt, John C. "Probabilistic Outputs for Support Vector Machines and Comparisons to Regularized Likelihood Methods". In: *Advances in Large Margin Classifiers* (1999), pp. 61–74.

[58] Poldrack, Russell A, Huckins, Grace, and Varoquaux, Gael. "Establishment of best practices for evidence for prediction: a review". In: *JAMA psychiatry* 77.5 (2020), pp. 534–540.

[59] Powers, Scott, Qian, Junyang, Jung, Kenneth, Schuler, Alejandro, Shah, Nigam H., Hastie, Trevor, and Tibshirani, Robert. "Some methods for heterogeneous treatment effect estimation in high dimensions". In: *Statistics in Medicine* 37.11 (2018), pp. 1767–1787.

[60] Rahimi, Ali and Recht, Benjamin. "Random Features for Large-Scale Kernel Machines". In: *Advances in Neural Information Processing Systems*. Vol. 20. 2008.

[61] Rajkomar, Alvin, Dean, Jeffrey, and Kohane, Isaac. "Machine learning in medicine". In: *New England Journal of Medicine* 380.14 (2019), pp. 1347–1358.

[62] Robins, James. "A new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect". In: *Mathematical Modelling* 7.9 (Jan. 1, 1986), pp. 1393–1512.

[63] Robins, James M. and Greenland, Sander. "The role of model selection in causal inference from non experimental data". In: *American Journal of Epidemiology* 123.3 (Mar. 1986), pp. 392–402.

[64] Robinson, P. M. "Root-N-Consistent Semiparametric Regression". In: *Econometrica* 4 (1988), pp. 931–954.

[65] Rolling, Craig A. and Yang, Yuhong. "Model selection for estimating treatment effects". In: *Journal of the Royal Statistical Society: Series B (Statistical Methodology)* 76.4 (Sept. 2014), pp. 749–769.
**Algorithm 2 Data simulation for the simulated dataset caussim**

Let the parameters of a given simulation be:

• N the number of samples.

• D the dimension of the basis expansion.

• p<sub>A</sub> the proportion of treated individuals.

• θ the overlap parameter.

• $\mathbb{P}(x|A = a) = W \cdot \mathcal{N} \left( \begin{bmatrix} (1-2a)\theta \\ 0 \end{bmatrix}; \begin{bmatrix} \sigma_0 & 0 \\ 0 & \sigma_1 \end{bmatrix} \right)$ the rotated gaussian for each population covariate. W is a random rotation matrix and $\sigma_0^2 = 2; \sigma_1^2 = 5$ are scaling parameters.

• $e(x) = \mathbb{P}(A = 1 | X = x)$ the oracle propensity score, obtained analytically from $\mathbb{P}(x|A = a)$.

• $\mathbb{P}(X) = p_A \mathbb{P}(X|A=1) + (1-p_A) \mathbb{P}(X|A=0)$ the mixture of gaussian for the whole population covariates.

• $b_1, \dots, b_D \sim \mathbb{P}(X)$, the basis sampled from the gaussian mixture.

• $z(x) = [\text{RBF}_\gamma(x, b_d)]_{d=1..D} \cdot Z^T \in \mathbb{R}^D$ the Nystroem expansion. $\text{RBF}_\gamma$ is the Gaussian kernel $K(x, y) = \exp(-\gamma ||x - y||^2)$ and Z is the normalization constant of the kernel basis, computed as the root inverse of the basis kernel $Z = [K(b_1, b_j)]_{i,j=1..D}^{-1/2}$.

• $\beta_\mu$ and $\beta_\tau$, the linear coefficients on top of the basis expansion. They are sampled from two normal distributions with mean 0 and unit variances $\mathcal{N}(0, I_{D+1})$.

• $\omega$ the treatment effect relative size with respect to the baseline response.

• $\varepsilon(x, a) \sim \mathcal{N}(0, \sigma_y)$ an exogen gaussian noise with $\sigma_y$ the scale of the noise.

Generation process for one sample $(x, a, e(x), y(0), y(1), y)$:

1: $x \sim \mathbb{P}(X)$ ▷ Sample from rotated gaussian mixture

2: $z(x)$ ▷ Compute the Nystroem expansion

3: $\mu_0(x) = [z(x);1] \cdot \beta_\mu^T$ ▷ Compute the non-treated response surface, ie. the baseline

4: $\tau(x) = [z(x);1] \cdot \beta_\tau^T$ ▷ Compute the conditional treatment effect (CATE)

5: $y(a) = (1-\omega)\mu_0(x) + a \omega \tau(x) + \varepsilon(x, a)$ ▷ Compute the potential outcomes

6: $a \sim \text{Bernoulli}(e(x))$ ▷ Compute the treatment status

7: $y = a y(1) + (1-a) y(0)$ ▷ Compute the observed outcome

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Example of the simulation setup in the input space with two knots – ie. basis functions. The top panel shows the observations in feature space, while the bottom panel displays the two response surfaces on a 1D cut along the black lines drawn on the top panel.</figcaption></figure>

We add a non-constant slope in the sigmoid to control the overlap between treated and control populations. We sampled uniformly 1000 different overlap parameters between 0 and 2.5, totaling 1000 dataset instances. Unlike the previous datasets, only the overlap varies for these instances. The response surfaces are set by the original outcomes.
[66] Rosenbaum, Paul R and Rubin, Donald B. "The central role of the propensity score in observational studies for causal effects". In: *Biometrika* 70 (1983), pp. 41–55.

[67] Rubin, Donald B. "Causal Inference Using Potential Outcomes". In: *Journal of the American Statistical Association* 100.469 (Mar. 1, 2005), pp. 322–331.

[68] Saito, Yuta and Yasui, Shota. "Counterfactual Cross-Validation: Stable Model Selection Procedure for Causal Inference Models". In: *International Conference on Machine Learning*. PMLR. 2020, pp. 8398–8407.

[69] Schulam, Peter and Saria, Suchi. "Reliable Decision Support using Counterfactual Models". In: *Advances in neural information processing systems* 30 (Mar. 30, 2017).

[70] Schuler, Alejandro, Baiocchi, Michael, Tibshirani, Robert, and Shah, Nigam. "A comparison of methods for model selection when estimating individual treatment effects". In: *arXiv:1804.05146 [cs, stat]* (June 13, 2018).

[71] Schuler, Megan S. and Rose, Sherri. "Targeted Maximum Likelihood Estimation for Causal Inference in Observational Studies". In: *American Journal of Epidemiology* 185.1 (Jan. 1, 2017), pp. 65–73.

[72] Shalit, Uri, Johansson, Fredrik D, and Sontag, David. "Estimating individual treatment effect: generalization bounds and algorithms". In: *International Conference on Machine Learning*. PMLR. 2017, pp. 3076–3085.

[73] Shen, Li, Margolies, Laurie R, Rothstein, Joseph H, Fluder, Eugene, McBride, Russell, and Sieh, Weiva. "Deep learning to improve breast cancer detection on screening mammography". In: *Scientific reports* 9.1 (2019), p. 12495.

[74] Shen, Lingjie, Geleijnse, Gijs, and Kaptein, Maurits. "RCTrep: An R Package for the Validation of Estimates of Average Treatment Effects". In: *Journal of Statistical Software* (2023).

[75] Shimoni, Yishai, Karavani, Ehud, Ravid, Sivan, Bak, Peter, Ng, Tan Hung, Alford, Sharon Hensley, Meade, Denise, and Goldschmidt, Yaara. "An Evaluation Toolkit to Guide Model Selection and Cohort Definition in Causal Inference". In: *arXiv preprint arXiv:1906.00442* (2019).

[76] Shimoni, Yishai, Yanover, Chen, Karavani, Ehud, and Goldschmидт, Yaara. "Benchmarking Framework for Performance-Evaluation of Causal Inference Analysis". In: *arXiv:1802.05046 [cs, stat]* (Mar. 20, 2018).

[77] Simon, Gregory E, Johnson, Eric, Lawrence, Jean M, Rossom, Rebecca C, Ahmedani, Brian, Lynch, Frances L, Beck, Arne, Waitzfelder, Beth, Ziebell, Rebecca, Penfold, Robert B, et al. "Predicting suicide attempts and suicide deaths following outpatient visits using electronic health records". In: *American Journal of Psychiatry* 175.10 (2018), pp. 951–960.

[78] Snowden, Jonathan M., Rose, Sherri, and Mortimer, Kathleen M. "Implementation of G-computation on a simulated data set: demonstration of a causal inference technique". eng. In: *American Journal of Epidemiology* 173.7 (Apr. 2011), pp. 731–738.

[79] Spasic, Irena, Nenadic, Goran, et al. "Clinical text data in machine learning: systematic review". In: *JMIR medical informatics* 8.3 (2020), e17984.

[80] Sperrin, Matthew, Jenkins, David, Martin, Glen P, and Peek, Niels. "Explicit causal reasoning is needed to prevent prognostic models being victims of their own success". In: *Journal of the American Medical Informatics Association* 26.12 (2019), pp. 1675–1676.
## Model selection procedures

**Nuisances estimation.** The nuisances are estimated with a stacked regressor inspired by the Super Learner framework, [90]). We select hyper-parameters with randomized search on a validation set $V$ and keep them fix for model selection. The searc grid is detailed in Table 3. All implementations come from scikit-learn [79]. As extreme inverse propensity weights induce high variance, we use clipping [91, 92] to bound $min(\varepsilon, 1 - \varepsilon)$ away from 0 with a fixed $\eta = 10^{-10}$, ensuring strict overlap for numerical stability.

<table><caption>**Table 3.** Hyper-parameters grid used for nuisance models</caption><thead><tr><th>Model</th><th>Estimator</th><th>Hyper-parameters grid</th></tr></thead><tbody><tr><td>Outcome, m</td><td>StackedRegressor<br>(HistGradientBoostingRegressor, ridge)</td><td>ridge regularization: [0.0001, 0.001, 0.01, 0.1, 1, 10, 100]<br>HistGradientBoostingRegressor learning rate: [0.01, 0.1, 1]<br>HistGradientBoostingRegressor max leaf nodes: [10, 20, 30, 50]</td></tr><tr><td>Treatment, e</td><td>StackedClassifier<br>(HistGradientBoostingClassifier, LogisticRegression)</td><td>LogisticRegression C: [0.0001, 0.001, 0.01, 0.1, 1, 10, 100]<br>HistGradientBoostingClassifier learning rate: [0.01, 0.1, 1]<br>HistGradientBoostingClassifier max leaf nodes: [10, 20, 30, 50]</td></tr></tbody></table>

## Additional Results

**Definition of the Kendall's tau, $\kappa$.** The Kendall's tau is a widely used statistics to measure the rank correlation between two set of observations. It measures the number of concordant pairs minus the discordant pairs normalized by the total number of pairs. It takes values in the $[-1, 1]$ range.

$$ \kappa = \frac{(\text{number of concordant pairs}) - (\text{number of discordant pairs})}{(\text{number of pairs})} \qquad (20) $$

**Values of relative $\kappa(\ell, \tau$-risk) compared to the mean over all metrics Kendall's as shown in the boxplots of Figure 3.**

**Figure 14 - Results measured in relative Kendall's for feasible and semi-oracle risks.** Because of extreme propensity scores in the denominator and bayes error residuals in the numerator, the semi-oracle $U$-risk has poor performances at bad overlap. Estimating these propensity scores in the is feasible $U$-risk reduces the variance since clipping is performed.

**Figure 15 - Results measured in absolute Kendall's.**

**Figure 16 - Results measured as distance to the oracle tau-risk.** To see practical gain in term of $\tau$-risk, we plot the results as the normalized distance between the estimator selected by the oracle $\tau$-risk and the estimator selected by each causal metric.

Then, $\widehat{R}$-risk* is more efficient than all other metrics. The gain are substantial for every datasets.

<figure><img src="image_16.png" /><figcaption>**Figure 14.** The $R$-risk is the best metric on average: Relative Kendall's $\tau$ agreement with $\tau$-risk. Strong and Weak overlap correspond to the first and last tertiles of the overlap distribution measured with Normalized Total Variation eq. 17.</figcaption></figure>
[81] Sriperumbudur, Bharath K., Fukumizu, Kenji, Gretton, Arthur, Schölkopf, Bernhard, and Lanckriet, Gert R. G. "On integral probability metrics, \phi-divergences and binary classification". In: *arXiv:0901.2698 [cs, math]* (Oct. 12, 2009).

[82] Steyerberg, Ewout W. *Clinical prediction models*. Springer, 2019.

[83] Su, Xiaogang, Peña, Annette T, Liu, Lei, and Levine, Richard A. "Random forests of interaction trees for estimating individualized treatment effects in randomized trials". In: *Statistics in medicine* 37.17 (2018), pp. 2547–2560.

[84] Swaminathan, Adith and Joachims, Thorsten. "Counterfactual risk minimization: Learning from logged bandit feedback". In: *International Conference on Machine Learning*. PMLR. 2015, pp. 814–823.

[85] VanderWeele, Tyler J. "Principles of confounder selection". In: *European journal of epidemiology* 34 (2019), pp. 211–219.

[86] Varoquaux, Gaël and Colliot, Olivier. *Evaluating machine learning models and their diagnostic value*. 2022.

[87] Wager, Stefan and Athey, Susan. "Estimation and Inference of Heterogeneous Treatment Effects using Random Forests". In: *Journal of the American Statistical Association* 113.523 (July 3, 2018), pp. 1228–1242.

[88] Wang, Hanyin, Li, Yikuan, Khan, Seema A, and Luo, Yuan. "Prediction of breast cancer distant recurrence using natural language processing and knowledge-guided convolutional neural network". In: *Artificial intelligence in medicine* 110 (2020), p. 101977.

[89] Wendling, T., Jung, K., Callahan, A., Schuler, A., Shah, N. H., and Gallego, B. "Comparing methods for estimation of heterogeneous treatment effects using observational data from health care databases". In: *Statistics in Medicine* 23 (2018), pp. 3309–3324.

[90] Wyatt, Jeremy C and Altman, Douglas G. "Commentary: Prognostic models: clinically useful or quickly forgotten?" In: *Bmj* 311.7019 (1995), pp. 1539–1541.

[91] Yala, Adam, Lehman, Constance, Schuster, Tal, Portnoi, Tally, and Barzilay, Regina. "A deep learning mammography-based model for improved breast cancer risk prediction". In: *Radiology* 292.1 (2019), pp. 60–66.

[92] Zadrozny, Bianca and Elkan, Charles. "Obtaining calibrated probability estimates from decision trees and naive Bayesian classifiers". In: (2001), p. 8.

[93] Zhang, Zhenwei and Sejdić, Ervin. "Radiological images and machine learning: trends, perspectives, and prospects". In: *Computers in biology and medicine* 108 (2019), pp. 354–370.
<table><caption>Table 4. Values of relative $κ(ℓ, τ-risk)$ compared to the mean over all metrics Kendall's as shown in the boxplots of Figure 3</caption><thead><tr><th rowspan="2">Metric</th><th rowspan="2">Dataset</th><th colspan="2">Strong Overlap</th><th colspan="2">Weak Overlap</th></tr><tr><th>Median</th><th>IQR</th><th>Median</th><th>IQR</th></tr></thead><tbody><tr><td rowspan="4">μ-risk</td><td>Twins (N=11 984)</td><td>-0.32</td><td>0.12</td><td>-0.19</td><td>0.12</td></tr><tr><td>ACIC 2016 (N=4 802)</td><td>-0.03</td><td>0.13</td><td>0.11</td><td>0.19</td></tr><tr><td>Caussim (N=5 000)</td><td>-0.40</td><td>0.55</td><td>-0.16</td><td>0.31</td></tr><tr><td>ACIC 2018 (N=5 000)</td><td>0.00</td><td>0.30</td><td>0.01</td><td>0.40</td></tr><tr><td rowspan="4">μ-risk<sub>IPW</sub></td><td>Twins (N=11 984)</td><td>-0.31</td><td>0.13</td><td>-0.17</td><td>0.12</td></tr><tr><td>ACIC 2016 (N=4 802)</td><td>-0.02</td><td>0.13</td><td>0.11</td><td>0.19</td></tr><tr><td>Caussim (N=5 000)</td><td>-0.34</td><td>0.50</td><td>0.09</td><td>0.31</td></tr><tr><td>ACIC 2018 (N=5 000)</td><td>0.00</td><td>0.30</td><td>-0.01</td><td>0.43</td></tr><tr><td rowspan="4">μ-risk<sub>IPW</sub><sup>*</sup></td><td>Twins (N=11 984)</td><td>-0.32</td><td>0.13</td><td>-0.17</td><td>0.13</td></tr><tr><td>ACIC 2016 (N=4 802)</td><td>-0.02</td><td>0.13</td><td>0.11</td><td>0.21</td></tr><tr><td>Caussim (N=5 000)</td><td>-0.33</td><td>0.54</td><td>0.26</td><td>0.27</td></tr><tr><td>Twins (N=11 984)</td><td>0.13</td><td>0.12</td><td>0.27</td><td>0.12</td></tr><tr><td rowspan="4">τ-risk<sub>IPW</sub></td><td>ACIC 2016 (N=4 802)</td><td>-0.07</td><td>0.18</td><td>0.05</td><td>0.31</td></tr><tr><td>Caussim (N=5 000)</td><td>-0.19</td><td>0.43</td><td>-0.14</td><td>0.18</td></tr><tr><td>ACIC 2018 (N=5 000)</td><td>-0.16</td><td>0.40</td><td>-0.11</td><td>0.66</td></tr><tr><td>Twins (N=11 984)</td><td>0.12</td><td>0.14</td><td>0.20</td><td>0.16</td></tr><tr><td rowspan="4">τ-risk<sub>IPW</sub><sup>*</sup></td><td>ACIC 2016 (N=4 802)</td><td>-0.03</td><td>0.16</td><td>-0.09</td><td>0.43</td></tr><tr><td>Caussim (N=5 000)</td><td>-0.15</td><td>0.46</td><td>-0.17</td><td>0.19</td></tr><tr><td>Twins (N=11 984)</td><td>0.13</td><td>0.12</td><td>0.02</td><td>0.25</td></tr><tr><td>ACIC 2016 (N=4 802)</td><td>0.04</td><td>0.11</td><td>0.11</td><td>0.26</td></tr><tr><td rowspan="3">U-risk</td><td>Caussim (N=5 000)</td><td>0.04</td><td>0.43</td><td>-0.04</td><td>0.17</td></tr><tr><td>ACIC 2018 (N=5 000)</td><td>0.12</td><td>0.26</td><td>-0.02</td><td>0.50</td></tr><tr><td>Twins (N=11 984)</td><td>0.25</td><td>0.08</td><td>-0.41</td><td>0.45</td></tr><tr><td rowspan="3">U-risk<sup>*</sup></td><td>ACIC 2016 (N=4 802)</td><td>0.08</td><td>0.13</td><td>-0.59</td><td>0.57</td></tr><tr><td>Caussim (N=5 000)</td><td>0.46</td><td>0.12</td><td>0.02</td><td>0.44</td></tr><tr><td>Twins (N=11 984)</td><td>0.15</td><td>0.10</td><td>0.25</td><td>0.18</td></tr><tr><td rowspan="3">R-risk</td><td>ACIC 2016 (N=4 802)</td><td>0.07</td><td>0.12</td><td>0.22</td><td>0.15</td></tr><tr><td>Caussim (N=5 000)</td><td>0.34</td><td>0.26</td><td>0.13</td><td>0.21</td></tr><tr><td>ACIC 2018 (N=5 000)</td><td>0.13</td><td>0.27</td><td>0.21</td><td>0.47</td></tr><tr><td rowspan="3">R-risk<sup>*</sup></td><td>Twins (N=11 984)</td><td>0.25</td><td>0.10</td><td>0.32</td><td>0.15</td></tr><tr><td>ACIC 2016 (N=4 802)</td><td>0.12</td><td>0.12</td><td>0.25</td><td>0.15</td></tr><tr><td>Caussim (N=5 000)</td><td>0.47</td><td>0.11</td><td>0.16</td><td>0.14</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>**Figure 7:** Different outcome models lead to different estimation errors on the Average Treatment Effects, on 77 classic simulations with known true causal effect¹⁸. The different models are ridge regression and random forests with different hyper-parameters (details A). The different configurations are plotted as a function of increasing difference between treated and untreated population –see section. There is no systematic best performer; data-driven model selection is important.</figcaption></figure>

## A Variability of ATE estimation on ACIC 2016

Figure 7 shows ATE estimations for six different models used in g-computation estimators on the 76 configurations of the ACIC 2016 dataset. Outcome models are fitted on half of the data and inference is done on the other half –ie. train/test with a split ratio of 0.5. For each configuration, and each model, this train test split was repeated ten times, yielding non parametric variance estimates¹¹. Figure 7 shows large variations obtained across different outcome estimators on semi-synthetic datasets¹⁸. Flexible models such as random forests are doing well in most settings except when treated and untreated populations differ noticeably, in which case a linear model (ridge) is to be preferred. However random forests with different hyper-parameters (max depth= 2) yield poor estimates. A simple rule of thumb such as preferring flexible models does not work in general; model selection is needed.

Outcome models are implemented with scikit-learn⁵⁴ and the following hyper-parameters:

<table><caption>**Table 2:** Hyper-parameters grid used for ACIC 2016 ATE variability</caption><thead><tr><td>Outcome Model</td><td>Hyper-parameters grid</td></tr></thead><tbody><tr><td>Random Forests</td><td>Max depth: [2, 10]</td></tr><tr><td>Ridge regression without treatment interaction</td><td>Ridge regularization: [0.1]</td></tr><tr><td>Ridge regression with treatment interaction</td><td>Ridge regularization: [0.1]</td></tr></tbody></table>
(a) Caussim

<figure><img src="image_17.png" /></figure>

(b) ACIC 2016

<figure><img src="image_18.png" /></figure>

(c) ACIC 2018

<figure><img src="image_19.png" /></figure>

(d) TWINS

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

**Figure 15.** Agreement with τ-risk ranking of methods function of overlap violation. The lines represent medians, estimated with a lowess. The transparent bands denote the 5% and 95% confidence intervals.
## B Prior work: model selection for outcome modeling (g-computation)

A natural way to select a predictive model for causal inference would be an error measure between a causal quantity such as the CATE and models' estimate. But such error is not a "feasible" risk: it cannot be computed solely from observed data and requires oracle knowledge.

**Simulation studies of causal model selection** Using eight simulations setups from$^{59}$, where the oracle CATE is known, Schuler et al. [70] compare four causal risks, concluding that for CATE estimation the best model-selection risk is the so-called R-risk$^{52}$-def. 5, below. Their empirical results are clear for randomized treatment allocation but less convincing for observational settings where both simple Mean Squared Error-MSE, μ-risk(f) def. 1- and reweighted MSE-μ-riskIPW def. 2- appear to perform better than R-risk on half of the simulations. Another work$^1$ studied empirically both MSE and reweighted MSE risks on the semi-synthetic ACIC 2016 datasets$^{18}$, but did not include the R-risk. We complete these prior empirical work by studying a wider variety of data generative processes and varying the influence of overlap, an important parameter of the data generation process which makes a given causal metric appropriate$^{16}$. We also study how to best adapt cross-validation procedures to causal metrics which themselves come with models to estimate.

**Theoretical studies of causal model selection** Several theoretical works have proposed causal model selection procedures that are consistent: select the best model in a family given asymptotically large data. These work rely on introducing a CATE estimator in the testing procedure: matching$^{65}$, an IPW estimate$^{25}$, a doubly robust estimator$^{68}$, or debiasing the error with influence functions$^1$. However, for theoretical guarantees to hold, the test-set correction needs to converge to the oracle: it needs to be flexible enough -well-posed- and asymptotic data. From a practical perspective, meeting such requirements implies having a good CATE estimate, thus having solved the original problem of causal model selection.

**Statistical guarantees on causal estimation procedures** Much work in causal inference has focused on procedures that guarantee asymptotically consistent estimators, such as Targeted Machine Learning Estimation (TMLE)$^{40,71}$ or Double Machine Learning$^{14}$. Here also, theories require asymptotic regimes and models to be well-specified.

By contrast, Johansson et al. [34] studies causal estimation without assuming that estimators are well specified. They derive an upper bound on the oracle error to the CATE (τ-risk) that involves the error on the outcome and the similarity of the distributions of treated and control patients. However, they use this upper bound for model optimization, and do not give insights on model selection. In addition, for hyperparameter selection, they rely on a plugin estimate of the τ-risk built with counterfactual nearest neighbors, which has been shown ineffective$^{70}$.

## C Causal assumptions

We assume the following four assumptions, referred as strong ignorability and necessary to assure identifiability of the causal estimands with observational data$^{67}$:

**Assumption 1 (Unconfoundedness)**

$$ \{Y(0), Y(1)\} \perp A | X $$
<figure><img src="image_22.png" /><figcaption>(b) ACIC 2016</figcaption></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>(c) ACIC 2018</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(d) TWINS</figcaption></figure>

<figure><img src="image_26.png" /></figure>

**Figure 16.** Metric performances by normalized tau-risk distance to the best method selected with τ-risk. All nuisances are learned with the same estimator stacking gradient boosting and ridge regression. Dotted and plain lines corresponds to 60% lowess quantile estimates. This choice of quantile allows to see better the oracle metrics lines for which outliers with a value of 0 distord the curves.
This condition –also called ignorability– is equivalent to the conditional independence on e(X)⁶⁶: {Y(0), Y(1)} ⊥ A|e(X).

**Assumption 2 (Overlap, also known as Positivity))**

$$ \eta < e(x) < 1 - \eta \quad \forall x \in \mathcal{X} \text{ and some } \eta > 0 $$

The treatment is not perfectly predictable. Or with different words, every patient has a chance to be treated and not to be treated. For a given set of covariates, we need examples of both to recover the ATE.

As noted by¹⁶, the choice of covariates X can be viewed as a trade-off between these two central assumptions. A bigger covariates set generally reinforces the ignorability assumption. In the contrary, overlap can be weakened by large X because of the potential inclusion of instruments: variables only linked to the treatment which could lead to arbitrarily small propensity scores.

**Assumption 3 (Consistency) The observed outcome is the potential outcome of the assigned treatment:**

$$ Y = AY(1) + (1 - A)Y(0) $$

Here, we assume that the intervention A has been well defined. This assumption focuses on the design of the experiment. It clearly states the link between the observed outcome and the potential outcomes through the intervention²⁶.

**Assumption 4 (Generalization) The training data on which we build the estimator and the test data on which we make the estimation are drawn from the same distribution D*, also known as the “no covariate shift” assumption³³.**

# D Definitions of feasible risks

**Definition 1 (Factual μ-risk)** ⁷² This is the usual Mean Squared Error on the target y. It is what is typically meant by “generalization error” in supervised learning:

$$ \mu\text{-risk}(f) = \mathbb{E} \left[ (Y - f(X; A))^2 \right] $$

**Definition 2 (μ-risk<sub>IPW</sub><sup>*</sup>)** ³⁸ Let the inverse propensity weighting function w(x, a) = a/(e(x)) + 1/(1-e(x)), we define the semi-oracle Inverse Propensity Weighting risk,

$$ \mu\text{-risk}_{IPW}^{*}(f) = \mathbb{E} \left[ \left( \frac{A}{e(X)} + \frac{1-A}{1-e(X)} \right) (Y - f(X; A))^2 \right] $$

**Definition 3 (τ-risk<sub>IPW</sub><sup>*</sup>)** ⁸⁷ The CATE τ(x) can be estimated with a regression against inverse propensity weighted outcomes⁴,²⁵,⁸⁷, the τ-risk<sub>IPW</sub>.

$$ \tau\text{-risk}_{IPW}^{*}(f) = \mathbb{E} \left[ \left( Y \frac{A - e(X)}{e(X)(1 - e(X))} - \tau_f(X) \right)^2 \right] $$
*Figure 17 - Stacked models for the nuisances is more efficient.* For each metrics the benefit of using a stacked model of linear and boosting estimators for nuisances compared to a linear model. The evaluation measure is Kendall's tau relative to the oracle R-risk* to have a stable reference between experiments. Thus, we do not include in this analysis the ACIC 2018 dataset since R-risk* is not available due to the lack of the true propensity score.

*Figure 18 Low population overlap hinders model selection for all metrics.*

*Figure 19 - Stacked models for the nuisances is more efficient.* For each metrics the benefit of using a stacked model of linear and boosting estimators for nuisances compared to a linear model. The evaluation measure is Kendall's tau relative to the oracle R-risk* to have a stable reference between experiments. Thus, we do not include in this analysis the ACIC 2018 dataset since R-risk* is not available due to the lack of the true propensity score.

*Figure 20 - Flexible models are performant in recovering nuisances even in linear setups.*

*Selecting different seeds and parameters is crucial to draw conclusions.* One strength of our study is the various number of different simulated and semi-simulated datasets. We are convinced that the usual practice of using only a small number of generation processes does not allow to draw statistically significant conclusions.

Figure 21 illustrate the dependence of the results on the generation process for caussim simulations. We highlighted the different trajectories induced by three different seeds for data generation and three different treatment ratio instead of 1000 different seeds. The result curves are relatively stable from one setup to another for R-risk, but vary strongly for μ-risk and μ-risk<sub>IPW</sub>.
**Definition 4 (U-risk*)**$^{37,52}$ Based on the Robinson decomposition –eq. 7, the U-learner uses the A – e(X) term in the denominator. The derived risk is:

$$U-risk^*(f) = \mathbb{E}\left[\left(\frac{Y - m(X)}{A - e(X)} - \tau_f(X)\right)^2\right]$$

Note that extreme propensity weights in the denominator term might inflate errors in the numerator due to imperfect estimation of the mean outcome m.

**Definition 5 (R-risk*)**$^{52,70}$ The R-risk also uses two nuisance m and e:

$$R-risk^*(f) = \mathbb{E}\left[(Y - m(X)) - (A - e(X))\tau_f(X)\right]^2$$

It is also based on the Robinson decomposition –eq. 7.

# E Proofs: Links between feasible and oracle risks

## E.1 Reformulation of the R-risk as reweighted τ-risk

**Proposition 1 (R-risk as reweighted τ-risk) Proof 1** We consider the R-decomposition:$^{64}$,

$$y(a) = m(x) + (a - e(x))\tau(x) + \varepsilon(x; a) \quad (11)$$

Where $\mathbb{E}[\varepsilon(X; A)|X, A] = 0$ We can use it as plug in the R-risk formula:

$$\begin{align*}
R-risk(f) &= \int_{Y \times X \times A} [(y - m(x)) - (a - e(x))\tau_f(x)]^2 p(y; x; a) dydxda \\
&= \int_{Y \times X \times A} [(a - e(x))\tau(x) + \varepsilon(x; a) - (a - e(x))\tau_f(x)]^2 p(y; x; a) dydxda \\
&= \int_{X \times A} (a - e(x))^2 (\tau(x) - \tau_f(x))^2 p(x; a) dxda \\
&\quad + 2 \int_{Y \times X \times A} (a - e(x)) (\tau(x) - \tau_f(x)) \int_{Y} \varepsilon(x; a) p(y | x; a) dyp(x; a) dxda \\
&\quad + \int_{X \times A} \int_{Y} \varepsilon^2(x; a) p(y | x; a) dyp(x; a) dxda
\end{align*}$$

The first term can be decomposed on control and treated populations to force e(x) to appear:

$$\begin{align*}
& \int_X (\tau(x) - \tau_f(x))^2 [e(x)^2 p(x; 0) + (1 - e(x))^2 p(x; 1)] dx \\
&= \int_X (\tau(x) - \tau_f(x))^2 [e(x)^2 (1 - e(x)) p(x) + (1 - e(x))^2 e(x) p(x)] dx \\
&= \int_X (\tau(x) - \tau_f(x))^2 (1 - e(x)) e(x) [1 - e(x) + e(x)] p(x) dx \\
&= \int_X (\tau(x) - \tau_f(x))^2 (1 - e(x)) e(x) p(x) dx.
\end{align*}$$

The second term is null since, $\mathbb{E}[\varepsilon(x, a)|X, A] = 0$.

The third term corresponds to the modulated residuals : $\tilde{\sigma}_B^2(0) + \tilde{\sigma}_B^2(1)$
(a) Caussim

<figure><img src="image_27.png" /></figure>

(b) ACIC 2016

<figure><img src="image_28.png" /></figure>

(c) Twins

<figure><img src="image_29.png" /><figcaption>Figure 17. Results are similar between the Shared nuisances/candidate set and the Separated nuisances set procedure. The experience has not been run on the full metrics for Caussim due to computation costs.</figcaption></figure>
## E.2 Interesting special cases

**Randomization special case** If the treatment is randomized as in RCTs, $p(A = 1 | X = x) = p(A = 1) = p_A$, thus $\mu$-risk$_{IPW}$ takes a simpler form:

$$ \mu\text{-risk}_{IPW} = \mathbb{E}_{(Y,X,A)\sim\mathcal{D}} \left[ \left( \frac{A}{p_A} + \frac{1-A}{1-p_A} \right) (Y - f(X; A))^2 \right] $$

However, we still can have large differences between $\tau$-risk and $\mu$-risk$_{IPW}$ coming from heterogeneous errors between populations as shown experimentally in Schuler et al. [70] and our results below.

Concerning the R-risk, replacing $e(x)$ by its randomized value $p_A$ in Proposition 1 yields the oracle $\tau$-risk up to multiplicative and additive constants:

$$ \text{R-risk} = p_A (1 - p_A) \tau\text{-risk} + (1 - p_A) \sigma_B^2(0) + p_A \sigma_B^2(1) $$

Thus, selecting estimators with R-risk* in randomized setting controls the $\tau$-risk. This explains the strong performances of R-risk in randomized setups$^{70}$ and is a strong argument to use it to estimate heterogeneity in RCTs.

**Oracle Bayes predictor** If we have access to the oracle Bayes predictor for the outcome ie. $f(x, a) = \mu(x, a)$, then all risks are equivalent up to the residual variance:

$$ \tau\text{-risk}(\mu) = \mathbb{E}_{X\sim p(X)}[(\tau(X) - \tau_{\mu}(X))^2] = 0 \quad (12) $$

$$ \begin{aligned} \mu\text{-risk}(\mu) &= \mathbb{E}_{(Y,X,A)\sim p(Y;X;A)}[(Y - \mu_A(X))^2] \\ &= \int_{X,A} \varepsilon(x, a)^2 p(a | x) p(x) dx da \le \sigma_B^2(0) + \sigma_B^2(1) \end{aligned} \quad (13) $$

$$ \mu\text{-risk}_{IPW}(\mu) = \sigma_B^2(0) + \sigma_B^2(1) \quad \text{from Lemma ??} \quad (14) $$

$$ \text{R-risk}(\mu) = \tilde{\sigma}_B^2(0) + \tilde{\sigma}_B^2(1) \le \sigma_B^2(0) + \sigma_B^2(1) \quad \text{from Proposition 1} \quad (15) $$

Thus, differences between causal risks only matter in finite sample regimes. Universally consistent learners converge to the Bayes risk in asymptotic regimes, making all model selection risks equivalent. In practice however, choices must be made in non-asymptotic regimes.

## F Measuring overlap

**Motivation of the Normalized Total Variation** Overlap is often assessed by comparing visually population distributions as in Figure 1 or computing standardized difference on each feature$^{6,7}$. While these methods are useful to decide if positivity holds, they do not yield a single measure. Rather, we compute the divergence between the population covariate distributions $\mathbb{P}(X|A = 0)$ and $\mathbb{P}(X|A = 1)^{16,34}$. Computing overlap when working only on samples of the observed distribution, outside of simulation, requires a sophisticated estimator of discrepancy between distributions, as two data points never have the same exact set of features. Maximum Mean Discrepancy$^{22}$ is typically used in the context of causal
<figure><img src="image_30.png" /><figcaption>**Figure 18.** Low population overlap hinders causal model selection for all metrics: Kendall's τ agreement with τ-risk. Strong, medium and Weak overlap correspond to the tertiles of the overlap distribution measured with Normalized Total Variation eq. 17.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>**Figure 19.** Learning the nuisances with stacked models (linear and gradient boosting) is important for successful model selection with R-risk. For Twins dataset, there is no improvement for stacked models compared to linear models because of the linearity of the propensity model.</figcaption></figure>
inference<sup>72,34</sup>. However it needs a kernel, typically Gaussian, to extrapolate across neighboring observations. We prefer avoiding the need to specify such a kernel, as it must be adapted to the data which is tricky with categorical or non-Gaussian features, a common situation for medical data.

For simulated and some semi-simulated data, we have access to the probability of treatment for each data point, which sample both densities in the same data point. Thus, we can directly use distribution discrepancy measures and rely on the Normalized Total Variation (NTV) distance to measure the overlap between the treated and control propensities. This is the empirical measure of the total variation distance<sup>81</sup> between the distributions, $TV(P(X|A = 1), P(X|A = 0))$. As we have both distribution sampled on the same points, we can rewrite it a sole function of the propensity score, a low dimensional score more tractable than the full distribution $P(X|A)$:

$$ \widehat{NTV}(e, 1-e) = \frac{1}{2N} \sum_{i=1}^{N} \left| \frac{e(x_i)}{p_A} - \frac{1-e(x_i)}{1-p_A} \right| \quad (16) $$

Formally, we can rewrite NTV as the Total Variation distance between the two population distributions. For a population $O = (Y(A), X, A) \sim D$:

$$ \begin{aligned} NTV(O) &= \frac{1}{2N} \sum_{i=1}^{N} \left| \frac{e(x_i)}{p_A} - \frac{1-e(x_i)}{1-p_A} \right| \\ &= \frac{1}{2N} \sum_{i=1}^{N} \left| \frac{P(A=1|X=x_i)}{p_A} - \frac{P(A=0|X=x_i)}{1-p_A} \right| \end{aligned} $$

Thus NTV approximates the following quantity in expectation over the data distribution $D$:

$$ \begin{aligned} NTV(D) &= \int_X \left| \frac{p(A=1|X=x)}{p_A} - \frac{p(A=0|X=x)}{1-p_A} \right| p(x) dx \\ &= \int_X \left| \frac{p(A=1, X=x)}{p_A} - \frac{p(A=0, X=x)}{1-p_A} \right| dx \\ &= \int_X \left| p(X=x|A=1) - p(X=x|A=0) \right| dx \end{aligned} $$

For countable sets, this expression corresponds to the Total Variation distance between treated and control populations covariate distributions: $TV(p_0(x), p_1(x))$.

**Measuring overlap without the oracle propensity scores:** For ACIC 2018, or for non-simulated data, the true propensity scores are not known. To measure overlap, we rely on flexible estimations of the Normalized Total Variation, using gradient boosting trees to approximate the propensity score. Empirical arguments for this plug-in approach is given in Figure 8.

**Empirical arguments** We show empirically that NTV is an appropriate measure of overlap by:

* Comparing the NTV distance with the MMD for Caussim which is gaussian distributed in Figure 10,

* Verifying that setups with penalized overlap from ACIC 2016 have a higher total variation distance than unpenalized setups in Figure 9.
<figure><img src="image_32.png" /><figcaption>**Figure 20.** Flexible models are performant in recovering nuisances in the downsampled Twins dataset. The propensity score is linear in this setup, making it particularly challenging for flexible models compared to linear methods.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>**Figure 21.** Kendall correlation coefficients for each causal metric. Each (color, shape) pair indicates a different (treatment ratio, seed) of the generation process.</figcaption></figure>
* Verifying that the Inverse Propensity Weights extrema (the inverse of the $\nu$ overlap constant appearing in the overlap Assumption 2) positevely correlates with NTV for Caussim, ACIC 2016 and Twins in Figure 11. Even if the same value of the maximum IPW could lead to different values of NTV, we expect both measures to be correlated : the higher the extrem propensity weights, the higher the NTV.

**Estimating NTV in practice** Finally, we verify that approximating the NTV distance with a learned plug-in estimates of $e(x)$ is reasonable. We used either a logistic regression or a gradient boosting classifier to learn the propensity models for the three datasets where we have access to the ground truth propensity scores: Caussim, Twins and ACIC 2016. We respectively sampled 1000, 1000 and 770 instances of these datasets with different seeds and overlap settings. We first run a hyperparameter search with cross-validation on the train set, then select the best estimator. We refit on the train set this estimator with or without calibration by cross validation and finally estimate the normalized TV with the obtained model. This training procedure reflects the one described in Algorithm 1 where nuisance models are fitted only on the train set.

The hyper parameters are : learning rate $\in [1e-3, 1e-2, 1e-1, 1]$, minimum samples leaf $\in [2, 10, 50, 100, 200]$ for boosting and L2 regularization $\in [1e-3, 1e-2, 1e-1, 1]$ for logistic regression.

Results in Figure 8 comparing bias to the true normalized Total Variation of each dataset instances versus growing true NTV indicate that calibration of the propensity model is crucial to recover a good approximation of the NTV.

# G Experiments

## G.1 Details on the data generation process

We use Gaussian-distributed covariates and random basis expansion based on Radial Basis Function kernels. A random basis of RBF kernel enables modeling non-linear and complex relationships between covariates in a similar way to the well known spline expansion. The estimators of the response function are learned with a linear model on another random basis (which can be seen as a stochastic approximation of the full data kernel$^{60}$). We carefully control the amount of overlap between treated and control populations, a crucial assumption for causal inference. Figure 12 illustrates 2D examples of the simulation.

* The raw features for both populations are drawn from a mixture of Gaussians: $P(X) = p_A P(X|A = 1) + (1 - p_A)P(X|A = 0)$ where $P(x|A = a)$ is a rotated Gaussian:

$$P(x|A = a) = W \cdot \mathcal{N} \left( \begin{bmatrix} (1-2a)\theta \\ 0 \end{bmatrix}; \begin{bmatrix} \sigma_0 & 0 \\ 0 & \sigma_1 \end{bmatrix} \right) \quad (17)$$

with $\theta$ a parameter controlling overlap (bigger yields poorer overlap), $W$ a random rotation matrix and $\sigma_0^2 = 2; \sigma_1^2 = 5$.

This generation process allows to analytically compute the oracle propensity scores $e(x)$, to simply control for overlap with the parameter $\theta$, the distance between the two Gaussian main axes and to visualize response surfaces.

* A basis expansion of the raw features increases the problem dimension. Using Radial Basis Func-
<figure><img src="image_34.png" /><figcaption>Figure 22. a) For CATE, a train/test ratio of 0.9/0.1 appears a good trade-off. b) For ATE, there is a small signal pointing also to 0.9/0.1 (K=10). for ATE. Experiences on 10 replications of all 78 instances of the ACIC 2016 data.</figcaption></figure>

## A.8 Data split choices

### Use 90% of the data to estimate outcome models, 10% to select them

The analyst faces a compromise: given a finite data sample, should she allocate more data to estimate the outcome model, thus improving the quality of the outcome model but leaving little data for model selection. Or, she could choose a bigger test set for model selection and effect estimation. For causal model selection, there is no established practice (as reviewed in A.8).

We investigate such tradeoff varying the ratio between train and test data size. For this, we first split out 30% of the data as a holdout set $V$ on which we use the oracle response functions to derive silver-standard estimates of causal quantities. We then use the standard estimation procedure on the remaining 70% of the data, splitting it into train $\tau$ and test $S$ of varying sizes. We finally measure the error between this estimate and the silver standard.

We consider two different analytic goals: estimating a average treatment effect –a single number used for policy making– and a CATE –a full model of the treatment effect as a function of covariates $X$. Given that the latter is a much more complex object than the former, the optimal train/test ratio might vary. To measure errors, we use for the ATE the relative absolute ATE bias between the ATE computed with the selected outcome model on the test set, and the true ATE as evaluated on the holdout set $V$. For the CATE, we compare the $\tau$-risk of the best selected model applied on the holdout set $V$. We explore this trade-off for the ACIC 2016 dataset and the R-risk.

Figure 22 shows that a train/test ratio of 0.9/0.1 (K=10) or 0.8/0.2 (K=5) appears best to estimate CATE and ATE.

### Heterogeneity in practices for data split

Splitting the data is common when using machine learning for causal inference, but practices vary widely in terms of the fraction of data to allocate to train models, outcomes and nuisances, and to evaluate them.

Before even model selection, data splitting is often required for estimation of the treatment effect, ATE or CATE, for instance to compute the nuisances required to optimize the outcome model (as the R-risk, Definition 5). The most frequent choice is use 80% of the data to fit the models, and 20% to evaluate them. For instance, for CATE estimation, the R-learner has been introduced using K-folds with $K = 5$ and $K = 10$: 80% of the data (4 folds) to train the nuisances and the remaining fold to minimize the corresponding R-loss [42]. Yet, it has been implemented with $K=5$ in causallib [93] or $K=3$ in econML [94]. Likewise, for ATE estimation, [43] introduce doubly-robust machine learning, recommending $K=5$ based on an empirical comparison $K=2$. However, subsequent works use doubly robust ML with varying choices of $K$: [95] use $K=3$, [96] use $K=2$. In the econML implementation, $K$ is set to 3 [94]. [97] evaluate various machine-learning approaches –including R-learners– using $K=5$ and 10, drawing inspiration from the TMLE literature which sets $K=5$ in the TMLE package [98].

Causal model selection has been much less discussed. The only study that we are aware of, [50], use a different data split: a 2-folds train/test procedure, training the nuisances on the first half of the data, and using the second half to estimate the R-risk and select the best treatment effect model.

## A.9 Sensitivity of the results to different ratio of causal effect to baseline response

Going beyond overlap, we study the effect of another important parameter of the data generation process: the ratio of causal effect to the baseline response.

Here, we detail the experimental setup of the results presented in Figure 7.

For each dataset of our simulation (caussim), we measure this ratio as the absolute mean difference between the causal effect and the baseline response, measured empirically on each dataset instance $\Delta_{\mu} = \frac{1}{N} \sum_{i=1}^{N} \left| \frac{\mu_{1}(x_{i}) - \mu_{0}(x_{i})}{\mu_{0}(x_{i}) + \mu_{1}(x_{i}) - \frac{1}{N} \sum_{j=1}^{N} \mu_{0}(x_{j}) + \mu_{1}(x_{j})} \right|$. We vary a dedicated parameter ($\omega$ in Algorithm 2) in the simulation to force this ratio to cover a wide range of values. For 179 dataset instances, it ranges several order of magnitudes with the full distribution of values given in Table 5.

<table><caption>Table 5. Distribution of the causal effect ratio $\Delta_{\mu}$ for the experiment on causal effect ratio, results in Figure 7.</caption><thead><tr><th></th><th>count</th><th>mean</th><th>std</th><th>min</th><th>1%</th><th>10%</th><th>25%</th><th>50%</th><th>75%</th><th>90%</th><th>99%</th><th>max</th></tr></thead><tbody><tr><td>effect_ratio</td><td>179</td><td>16.95</td><td>24.92</td><td>0.043</td><td>0.054</td><td>0.32</td><td>1.34</td><td>9.12</td><td>22.55</td><td>39.97</td><td>101.12</td><td>206.53</td></tr></tbody></table>

We also studied another measure, focusing on the ratio of the variance of the causal effect to the variance of the mean shared response ($\mu_{0} + \mu_{1}$) of both population: $\frac{\text{var}(\mu_{1}(x_{i}) - \mu_{0}(x_{i}))}{\text{var}(\mu_{1}(x_{i}) + \mu_{0}(x_{i}))}$, covering a range from 0 to 11.08 with the full distribution described in Table 6. Figure 23 shows that for high values of the variance, R-risk is outperformed by the $\mu$-risk$_{IPW}$ and the $\tau$-risk$_{IPW}$. However, on average, the R-risk
<figure><img src="image_11.png" /><figcaption>(b) Calibrated classifiers</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8:** a) Without calibration, estimation of NTV is not trivial even for boosting models. b) Calibrated classifiers are able to recover the true Normalized Total Variation for all datasets where it is available.</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>**Figure 23.** R-risk is robust to a wide range of effect variance: Kendall's τ agreement with τ-risk. Strong, medium and Weak Causal effect variance are the tertiles of the ratio of causal effect variance to the mean outcome variance: Low [0; 0.04[, Medium [0.04; 1.04[, High [1.04; 11.08[.</figcaption></figure>

is still the better risk.

<table><caption>**Table 6.** Statistical summary of the causal effect variance.</caption><thead><tr><th></th><th>count</th><th>mean</th><th>std</th><th>min</th><th>1%</th><th>10%</th><th>25%</th><th>50%</th><th>75%</th><th>90%</th><th>99%</th><th>max</th></tr></thead><tbody><tr><td>effect-variance</td><td>179</td><td>1.013</td><td>1.781</td><td>1.680e-5</td><td>1.734e-5</td><td>4.875e-4</td><td>6.449e-3</td><td>0.789</td><td>1.160</td><td>1.807</td><td>9.641</td><td>11.079</td></tr></tbody></table>

## References

1. Karel GM Moons, Patrick Royston, Yvonne Vergouwe, Diederick E Grobbee, and Douglas G Altman. Prognosis and prognostic research: what, why, and how? *Bmj*, 338, 2009.

2. Ewout W Steyerberg. *Clinical prediction models*. Springer.

3. Andrew L Beam and Isaac S Kohane. Big data and machine learning in health care. *Jama*, 319(13):1317–1318, 2018.

4. Alvin Rajkomar, Jeffrey Dean, and Isaac Kohane. Machine learning in medicine. *New England Journal of Medicine*, 380(14):1347–1358, 2019.

5. M Khojaste-Sarakhsi, Seyedhamidreza Shahabi Haghighi, SMT Fatemi Ghomi, and Elena Marchiori. Deep learning for alzheimer's disease diagnosis: A survey. *Artificial Intelligence in Medicine*, page 102332, 2022.
<figure><img src="image_13.png" /><figcaption>**Figure 10:** Good correlation between overlap measured as normalized Total Variation and Maximum Mean Discrepancy (200 sampled Gaussian datasets)</figcaption></figure>

<figure><img src="image_14.png" /></figure>
6. Zhenwei Zhang and Ervin Sejdić. Radiological images and machine learning: trends, perspectives, and prospects. *Computers in biology and medicine*, 108:354-370, 2019.

7. Adam Yala, Constance Lehman, Tal Schuster, Tally Portnoi, and Regina Barzilay. A deep learning mammography-based model for improved breast cancer risk prediction. *Radiology*, 292(1):60-66, 2019.

8. Li Shen, Laurie R Margolies, Joseph H Rothstein, Eugene Fluder, Russell McBride, and Weiva Sieh. Deep learning to improve breast cancer detection on screening mammography. *Scientific reports*, 9(1):12495, 2019.

9. Ali Bou Nassif, Manar Abu Talib, Qassim Nasir, Yaman Afadar, and Omar Elgendy. Breast cancer detection using artificial intelligence techniques: A systematic literature review. *Artificial Intelligence in Medicine*, page 102276, 2022.

10. Stephen J Mooney and Vikas Pejaver. Big data in public health: terminology, machine learning, and privacy. *Annual review of public health*, 39:95, 2018.

11. Rishi J Desai, Shirley V Wang, Muthiah Vaduganathan, Thomas Evers, and Sebastian Schneeweiss. Comparison of machine learning methods with traditional models for use of administrative claims with electronic medical records to predict heart failure outcomes. *JAMA network open*, 3(1):e1918962-e1918962, 2020.

12. Gregory E Simon, Eric Johnson, Jean M Lawrence, Rebecca C Rossom, Brian Ahmedani, Frances L Lynch, Arne Beck, Beth Waitzfelder, Rebecca Ziebell, Robert B Penfold, et al. Predicting suicide attempts and suicide deaths following outpatient visits using electronic health records. *American Journal of Psychiatry*, 175(10):951-960, 2018.

13. Steven Horng, David A Sontag, Yoni Halpern, Yacine Jernite, Nathan I Shapiro, and Larry A Nathanson. Creating an automated trigger for sepsis clinical decision support at emergency department triage using machine learning. *PloS one*, 12(4):e0174708, 2017.

14. Hanyin Wang, Yikuan Li, Seema A Khan, and Yuan Luo. Prediction of breast cancer distant recurrence using natural language processing and knowledge-guided convolutional neural network. *Artificial intelligence in medicine*, 110:101977, 2020.

15. Irena Spasic, Goran Nenadic, et al. Clinical text data in machine learning: systematic review. *JMIR medical informatics*, 8(3):e17984, 2020.

16. Douglas G Altman, Yvonne Vergouwe, Patrick Royston, and Karel GM Moons. Prognosis and prognostic research: validating a prognostic model. *Bmj*, 338, 2009.

17. Russell A Poldrack, Grace Huckins, and Gael Varoquaux. Establishment of best practices for evidence for prediction: a review. *JAMA psychiatry*, 77(5):534-540, 2020.

18. Gaël Varoquaux and Olivier Colliot. Evaluating machine learning models and their diagnostic value, 2022.

19. Jeremy C Wyatt and Douglas G Altman. Commentary: Prognostic models: clinically useful or quickly forgotten? *Bmj*, 311(7019): 1539-1541, 1995.

20. Mark Alan Fontana, Stephen Lyman, Gourab K Sarker, Douglas E Padgett, and Catherine H MacLean. Can machine learning algorithms predict which patients will achieve minimally clinically important differences from total joint arthroplasty? *Clinical orthopaedics and related research*, 477(6):1267, 2019.

21. Jonathan M. Snowden, Sherri Rose, and Kathleen M. Mortimer. Implementation of G-computation on a simulated data set: demonstration of a causal inference technique. *American Journal of Epidemiology*, 173(7):731-738, April 2011.

22. Matthew Sperrin, David Jenkins, Glen P Martin, and Niels Peek. Explicit causal reasoning is needed to prevent prognostic models being victims of their own success. *Journal of the American Medical Informatics Association*, 26(12):1675-1676, 2019.

23. Tony Blakely, John Lynch, Koen Simons, Rebecca Bentley, and Sherri Rose. Reflection on modern methods: when worlds collide—prediction, machine learning and causal inference. *International journal of epidemiology*, 49(6):2058-2064, 2020.

24. MA Hernán and JM Robins. *Causal Inference: What If*. CRC Boca Raton, FL.

25. Tyler J VanderWeele. Principles of confounder selection. *European journal of epidemiology*, 34:211-219, 2019.

26. James M. Robins and Sander Greenland. The role of model selection in causal inference from non experimental data. *American Journal of Epidemiology*, 123(3):392-402.

27. T. Wendling, K. Jung, A. Callahan, A. Schuler, N. H. Shah, and B. Gallego. Comparing methods for estimation of heterogeneous treatment effects using observational data from health care databases. *Statistics in Medicine*, (23):3309-3324.

28. Paul R Rosenbaum and Donald B Rubin. The central role of the propensity score in observational studies for causal effects. *Biometrika*, 70:41-55.

29. Peter C. Austin and Elizabeth A. Stuart. Moving towards best practice when using inverse probability of treatment weighting (IPTW) using the propensity score to estimate causal treatment effects in observational studies. *Statistics in Medicine*, 34(28):3661-3679.

30. Sabrina Casucci, Li Lin, Sharon Hewner, and Alexander Nikolaev. Estimating the causal effects of chronic disease combinations on 30-day hospital readmissions based on observational medicaid data. *Journal of the American Medical Informatics Association*, 25(6): 670-678, 2018.

31. Elysia Grose, Samuel Wilson, Jeffrey Barkun, Kimberly Bertens, Guillaume Martel, Fady Balaa, and Jad Abou Khalil. Use of propensity score methodology in contemporary high-impact surgical literature. *Journal of the American College of Surgeons*, 230(1):101-112.e2.

32. Xiaogang Su, Annette T Peña, Lei Liu, and Richard A Levine. Random forests of interaction trees for estimating individualized treatment effects in randomized trials. *Statistics in medicine*, 37(17):2547-2560, 2018.

33. Andrea Lamont, Michael D Lyons, Thomas Jaki, Elizabeth Stuart, Daniel J Feaster, Kukatharmini Tharmaratnam, Daniel Oberski, Hemant Ishwaran, Dawn K Wilson, and M Lee Van Horn. Identification of predicted individual treatment effects in randomized clinical trials. *Statistical methods in medical research*, 27(1):142-157, 2018.

34. Jeroen Hoogland, Joanna IntHout, Michail Belias, Maroeska M Rovers, Richard D Riley, Frank E. Harrell Jr, Karel GM Moons, Thomas PA Debray, and Johannes B Reitsma. A tutorial on individualized treatment effect prediction from randomized trials with a binary endpoint. *Statistics in medicine*, 40(26):5961-5981, 2021.

35. Jennifer L. Hill. Bayesian nonparametric modeling for causal inference. *Journal of Computational and Graphical Statistics*, (1):217-240.

36. Mark J. van der Laan and Sherri Rose. *Targeted Learning*. Springer Series in Statistics.

37. Megan S. Schuler and Sherri Rose. Targeted maximum likelihood estimation for causal inference in observational studies. *American Journal of Epidemiology*, 185(1):65-73.

38. Scott Powers, Junyang Qian, Kenneth Jung, Alejandro Schuler, Nigam H. Shah, Trevor Hastie, and Robert Tibshirani. Some methods for heterogeneous treatment effect estimation in high dimensions. *Statistics in Medicine*, 37(11):1767-1787.

39. Stefan Wager and Susan Athey. Estimation and inference of heterogeneous treatment effects using random forests. *Journal of the*
<figure><img src="image_15.png" /><figcaption>**Figure 11:** Maximal value of Inverse Propensity Weights increases exponentially with the overlap as measure by Normalized Total Variation.</figcaption></figure>

tion (RBF) Nystroem transformation<fnref n="5" />, we expand the raw features into a transformed space. The basis expansion samples randomly a small number of representatives in the raw data. Then, it computes an approximation of the full N-dimensional kernel with these basis components, yielding the transformed features $z(x)$. The number of basis functions –ie. *knots*–, controls the complexity of the ground-truth response surfaces and treatment. We first use this process to draw the non-treated response surface $\mu_0$ and the causal effect $\tau$. We then draw the observations from a mixture of two Gaussians, for the treated and non-treated. We vary the separation between the two Gaussians to control the overlap between treated and non-treated populations, an important parameter for causal inference (related to $\eta$ in Proposition ??). Finally, we generate observed outcomes adding Gaussian noise.

More formally, we generate the basis following the original data distribution, $[b_1..b_D] \sim P(x)$, with $D=2$ in our simulations. Then, we compute an approximation of the full kernel of the data generation process $RBF(x, \cdot)$ with $x \sim P(x)$ with these representatives: $z(x) = [RBF_\gamma(x, b_d)]_{d=1..D} \cdot Z^T \in \mathbb{R}^D$ with $RBF_\gamma$ being the Gaussian kernel $K(x, y) = \exp(-\gamma\|x - y\|^2)$ and $Z$ the normalization constant of the kernel basis, computed as the root inverse of the basis kernel $Z =$

<footnote n="5">⁵We use the Sklearn implementation,⁵⁴</footnote>
*American Statistical Association*, 113(523):1228-1242.

40. Susan Athey, Julie Tibshirani, and Stefan Wager. Generalized random forests. *Annals of Statistics*, 47(2):1148-1178.

41. Sören R. Künzel, Jasjeet S. Sekhon, Peter J. Bickel, and Bin Yu. Metalearners for estimating heterogeneous treatment effects using machine learning. *Proceedings of the National Academy of Sciences*, 116(10):4156-4165.

42. Xinkun Nie and Stefan Wager. Quasi-oracle estimation of heterogeneous treatment effects. *Biometrika*, 108(2):299-319.

43. Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins. Double/debiased machine learning for treatment and structural parameters. *The Econometrics Journal*, page 71.

44. Gang Fang, Izabela E Annis, Jennifer Elston-Lafata, and Samuel Cykert. Applying machine learning to predict real-world individual treatment effects: insights from a virtual patient cohort. *Journal of the American Medical Informatics Association*, 26(10):977-988, 2019.

45. Vincent Dorie, Jennifer Hill, Uri Shalit, Marc Scott, and Dan Cervone. Automated versus do-it-yourself methods for causal inference: Lessons learned from a data analysis competition. *Statistical Science*, 34(1):43-68, 2019.

46. Florent Le Borgne, Arthur Chatton, Maxime Léger, Rémi Lenain, and Yohann Foucher. G-computation and machine learning for estimating the causal effects of binary exposure statuses on binary outcomes. *Scientific reports*, 11(1):1435, 2021.

47. Jinma Ren, Paul Cislo, Joseph C Cappelleri, Patrick Hlavacek, and Marco DiBonaventura. Comparing g-computation, propensity score-based weighting, and targeted maximum likelihood estimation for analyzing externally controlled trials with both measured and unmeasured confounders: a simulation study. *BMC Medical Research Methodology*, 23(1):18, 2023.

48. Jean Kaddour, Aengus Lynch, Qi Liu, Matt J Kusner, and Ricardo Silva. Causal machine learning: A survey and open problems. arxiv 2022. arXiv preprint arXiv:2206.15475, 2022.

49. Victor Chernozhukov, Christian Hansen, Nathan Kallus, Martin Spindler, and Vasilis Syrgkanis. Applied causal inference powered by ml and ai. arXiv preprint arXiv:2403.02467, 2024.

50. Alejandro Schuler, Michael Baiocchi, Robert Tibshirani, and Nigam Shah. A comparison of methods for model selection when estimating individual treatment effects. arXiv:1804.05146 [cs, stat].

51. Ahmed Alaa and Mihaela Van Der Schaar. Validating causal inference models via influence functions. *International Conference on Machine Learning*, pages 191-201.

52. Mary E. Charlson, Peter Pompei, Kathy L. Ales, and C.Ronald MacKenzie. A new method of classifying prognostic comorbidity in longitudinal studies: Development and validation. *Journal of Chronic Diseases*, 40(5):373-383.

53. James Robins. A new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect. *Mathematical Modelling*, 7(9):1393-1512.

54. Ashley I Naimi and Brian W Whitcomb. Defining and identifying average treatment effects. *American Journal of Epidemiology*, 2023.

55. Guido W. Imbens and Donald B. Rubin. *Causal inference in statistics, social, and biomedical sciences*. Cambridge University Press.

56. Donald B Rubin. Causal inference using potential outcomes. *Journal of the American Statistical Association*, 100(469):322-331.

57. P. M. Robinson. Root-n-consistent semiparametric regression. *Econometrica*, (4):931-954.

58. Mark J van der Laan, MJ Laan, and JM Robins. *Unified methods for censored longitudinal data and causality*. Springer Science & Business Media.

59. Peter Schulam and Suchi Saria. Reliable decision support using counterfactual models. *Advances in neural information processing systems*, 30.

60. Uri Shalit, Fredrik D Johansson, and David Sontag. Estimating individual treatment effect: generalization bounds and algorithms. In *International Conference on Machine Learning*, pages 3076-3085. PMLR, 2017.

61. Chanelle J. Howe, Stephen R. Cole, Daniel J. Westreich, Sander Greenland, Sonia Napravnik, and Joseph J. Eron. Splines for trend analysis and continuous confounder control. 22(6):874-875.

62. Aris Perperoglou, Willi Sauerbrei, Michal Abrahamowicz, and Matthias Schmid. A review of spline function procedures in r. 19(1):46.

63. Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In *Advances in Neural Information Processing Systems*, volume 20.

64. Lingjie Shen, Gijs Geleijnse, and Maurits Kaptein. Rctrep: An r package for the validation of estimates of average treatment effects. *Journal of Statistical Software*, 2023.

65. Kenneth R. Niswander and United States National Institute of Neurologialac Diseases and Stroke. *The Women and Their Pregnancies: The Collaborative Perinatal Study of the National Institute of Neurological Diseases and Stroke*. National Institute of Health. Google-Books-ID: AobdVhlhDQkc.

66. Yishai Shimoni, Chen Yanover, Ehud Karavani, and Yaara Goldschmiedt. Benchmarking framework for performance-evaluation of causal inference analysis. arXiv:1802.05046 [cs, stat].

67. M. F. MacDorman and J. O. Atkinson. Infant mortality statistics from the linked birth/infant death data set-1995 period data. *Monthly Vital Statistics Report*, 46(6 Suppl 2):1-22, February 1998.

68. Christos Louizos, Uri Shalit, Joris Mooij, David Sontag, Richard Zemel, and Max Welling. Causal effect inference with deep latent-variable models. *Advances in neural information processing systems*.

69. Douglas Almond, Kenneth Y. Chay, and David S. Lee. The Costs of Low Birth Weight. *The Quarterly Journal of Economics*, 120(3), 2005.

70. M. G. Kendall. A new measure of rank correlation. *Biometrika*, 30(1-2):81-93, June 1938.

71. Susan Athey and Guido Imbens. Recursive partitioning for heterogeneous causal effects. *Proceedings of the National Academy of Sciences*, 113(27):7353-7360, 2016.

72. Pierre Gutierrez and Jean-Yves Gerardy. Causal inference and uplift modeling a review of the literature. *Proceedings of The 3rd International Conference on Predictive Applications and APIs*, (67):14.

73. John C. Platt and John C. Platt. Probabilistic Outputs for Support Vector Machines and Comparisons to Regularized Likelihood Methods. *Advances in Large Margin Classifiers*, pages 61-74, 1999.

74. Bianca Zadrozny and Charles Elkan. Obtaining calibrated probability estimates from decision trees and naive Bayesian classifiers. page 8.

75. Alexandru Niculescu-Mizil and Rich Caruana. Predicting good probabilities with supervised learning. In *Proceedings of the 22nd international conference on Machine learning - ICML '05*, pages 625-632. ACM Press, 2005.

76. Matthias Minderer, Josip Djolonga, Rob Romijnders, Frances Hubis, Xiaohua Zhai, Neil Houlsby, Dustin Tran, and Mario Lucic. Revisiting the Calibration of Modern Neural Networks. *Advances in Neural Information Processing Systems*, 34:15682-15694, 2021.
<figure><img src="image_16.png" /><figcaption>**Figure 12:** Example of the simulation setup in the input space with two knots -ie. basis functions. The top panel shows the observations in feature space, while the bottom panel displays the two response surfaces on a 1D cut along the black lines drawn on the top panel.</figcaption></figure>

$$[K(b_i, b_j)]_{i,j \in 1..D}^{-1/2}$$

* Functions $\mu_0$, $\tau$ are distinct linear functions of the transformed features:

$$\mu_0(x) = [z(x); 1] \cdot \beta_{\mu}^T$$

$$\tau(x) = [z(x); 1] \cdot \beta_{\tau}^T$$

* Adding a Gaussian noise, $\varepsilon \sim N(0, \sigma(x; a))$, we construct the potential outcomes: $y(a) = \mu_0(x) + a \tau(x) + \varepsilon(x, a)$

We generated 1000 instances of this dataset with uniformly random overlap parameters $\theta \in [0, 2.5]$.

## G.2 Details on the semi-simulated datasets

**ACIC 2018** ¹⁸: The initial intervention was a child's birth weight ($A = 1$ if weight < 2.5kg), and outcome was the child's IQ after a follow-up period. The study contained $N = 4802$ data points with $D = 55$ features (5 binary, 27 count data, and 23 continuous). They simulated 77 different setups
77. Alexandre Perez-Lebel, Marine Le Morvan, and Gaël Varoquaux. Beyond calibration: estimating the grouping loss of modern neural networks. *arXiv preprint arXiv:2210.16315*, 2022.

78. Xavier Bouthillier, Pierre Delaunay, Mirko Bronzi, Assya Trofimov, Brennan Nichyporuk, Justin Szeto, Nazanin Mohammadi Sepahvand, Edward Raff, Kanika Madan, Vikram Voleti, Samira Ebrahimi Kahou, Vincent Michalski, Tal Arbel, Chris Pal, Gael Varoquaux, and Pascal Vincent. Accounting for variance in machine learning benchmarks. *Proceedings of Machine Learning and Systems*, 3:747–769.

79. Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jake Vanderplas, Alexandre Passos, David Cournapeau, Matthieu Brucher, Matthieu Perrot, and Édouard Duchesnay. Scikit-learn: Machine learning in python. *Journal of Machine Learning Research*, 12(85):2825–2830.

80. Alexander D'Amour, Peng Ding, Avi Feller, Lihua Lei, and Jasjeet Sekhon. Overlap in observational studies with high-dimensional covariates. *Journal of Econometrics*, 221(2):644–654, 2021.

81. Craig A. Rolling and Yuhong Yang. Model selection for estimating treatment effects. 76(4):749–769.

82. Yuta Saito and Shota Yasui. Counterfactual cross-validation: Stable model selection procedure for causal inference models. In *International Conference on Machine Learning*, pages 8398–8407. PMLR, 2020.

83. Fredrik D Johansson, Uri Shalit, Nathan Kallus, and David Sontag. Generalization bounds and representation learning for estimation of potential outcomes and causal effects. *The Journal of Machine Learning Research*, 23(1):7489–7538, 2022.

84. Dylan J Foster and Vasilis Syrgkanis. Orthogonal statistical learning. *The Annals of Statistics*, 51(3):879–908, 2023.

85. Andrew Jesson, Sören Mindermann, Uri Shalit, and Yarin Gal. Identifying causal-effect inference failure with uncertainty-aware models. *Advances in Neural Information Processing Systems*, 33:11637–11649.

86. Peter C. Austin. An Introduction to Propensity Score Methods for Reducing the Effects of Confounding in Observational Studies. *Multivariate Behavioral Research*, (3):399–424, May 2011.

87. Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf, and Alexander Smola. A kernel two-sample test. *The Journal of Machine Learning Research*, 13(1):723–773.

88. Bharath K. Sriperumbudur, Kenji Fukumizu, Arthur Gretton, Bernhard Schölkopf, and Gert R. G. Lanckriet. On integral probability metrics, \phi-divergences and binary classification. *arXiv preprint arXiv:0901.2698 [cs, math]*.

89. Alicia Curth, David Svensson, and James Weatherall. Really doing great at estimating CATE? a critical look at ML benchmarking practices in treatment effect estimation. *Neurips Process 2021*, page 14.

90. Mark J. van der Laan, Eric C. Polley, and Alan E. Hubbard. Super learner. 6.

91. Adith Swaminathan and Thorsten Joachims. Counterfactual risk minimization: Learning from logged bandit feedback. In *International Conference on Machine Learning*, pages 814–823. PMLR, 2015.

92. Edward L Ionides. Truncated importance sampling. *Journal of Computational and Graphical Statistics*, 17(2):295–311.

93. Yishai Shimoni, Ehud Karavani, Sivan Ravid, Peter Bak, Tan Hung Ng, Sharon Hensley Alford, Denise Meade, and Yaara Goldschmidt. An evaluation toolkit to guide model selection and cohort definition in causal inference. *arXiv preprint arXiv:1906.00442*, 2019.

94. Keith Battocchi, Eleanor Dillon, Maggie Hei, Greg Lewis, Paul Oka, Miruna Oprescu, and Vasilis Syrgkanis. EconML: A Python Package for ML-Based Heterogeneous Treatment Effects Estimation. https://github.com/microsoft/EconML.

95. Nicolas Loiseau, Paul Trichelair, Maxime He, Mathieu Andreux, Mikhail Zaslavskiy, Gilles Wainrib, and Michael G. B. Blum. External control arm analysis: an evaluation of propensity score approaches, G-computation, and doubly debiased machine learning. *BMC Medical Research Methodology*, 22, 2022.

96. Zijun Gao, Trevor Hastie, and Robert Tibshirani. Assessment of heterogeneous treatment effect estimation accuracy via matching. *Statistics in Medicine*, (17), 2021.

97. Ashley I Naimi, Alan E Mishler, and Edward H Kennedy. Challenges in obtaining valid causal effect estimates with machine learning algorithms. *American Journal of Epidemiology*, 2021.

98. Susan Gruber and Mark J. van der Laan. tmle: An R package for targeted maximum likelihood estimation. *Journal of Statistical Software*, 51(13):1–35, 2012.
varying parameters for treatment and response models, overlap, and interactions between treatment and covariates<fnref n="6" />. We used 10 different seeds for each setup, totaling 770 dataset instances.

**ACIC 2018**⁷⁶: Starting from data from the Linked Births and Infant Deaths Database (LBIDD)⁴⁴ with $D = 177$ covariates, treatment and outcome models are simulated with complex models to reflect different scenarios. The data do not provide the true propensity scores, so we evaluate only feasible metrics, which do not require this nuisance parameter. We used all 432 datasets<fnref n="7" /> of size $N = 5000$.

**Twins**⁴³: It is an augmentation of real data on twin births and mortality rates². There are $N = 11984$ samples (pairs of twins), and $D = 50$ covariates<fnref n="8" />. The outcome is the mortality and the treatment is the weight of the heavier twin at birth. This is a "true" counterfactual dataset¹⁵ in the sense that we have both potential outcomes with each twin. They simulate the treatment with a sigmoid model based on GESTAT10 (number of gestation weeks before birth) and x the 45 other covariates:

$$ 
\begin{gathered}
\mathbf{t}_i | \mathbf{x}_i, \mathbf{z}_i \sim \text{Bern} \left( \sigma \left( w_o^\top \mathbf{x} + w_h (\mathbf{z}/10 - 0.1) \right) \right) \quad (18) \\
\text{with } w_o \sim \mathcal{N}(0, 0.1 \cdot I), \quad w_h \sim \mathcal{N}(5, 0.1)
\end{gathered}
$$

We add a non-constant slope in the sigmoid to control the overlap between treated and control populations. We sampled uniformly 1000 different overlap parameters between 0 and 2.5, totaling 1000 dataset instances. Unlike the previous datasets, only the overlap varies for these instances. The response surfaces are set by the original outcomes.

## G.3 Model selection procedures

**Nuisances estimation** The nuisances are estimated with a stacked regressor inspired by the Super Learner framework,³⁹). We select hyper-parameters with randomized search on a validation set $V$ and keep them fix for model selection. The searc grid is detailed in Table 3. All implementations come from scikit-learn⁵⁴. As extreme inverse propensity weights induce high variance, we use clipping⁸⁴,³² to bound $min(\tilde{e}, 1 - \tilde{e})$ away from 0 with a fixed $\eta = 10^{-10}$, ensuring strict overlap for numerical stability.

<table><caption>**Table 3:** Hyper-parameters grid used for nuisance models</caption><thead><tr><th>Model</th><th>Estimator</th><th>Hyper-parameters grid</th></tr></thead><tbody><tr><td>Outcome, m</td><td>StackedRegressor<br>(HistGradientBoostingRegressor, ridge)</td><td>ridge regularization: [0.0001, 0.001, 0.01, 0.1, 1]<br>HistGradientBoostingRegressor learning rate: [0.0001, 0.001, 0.01, 0.1, 1]</td></tr><tr><td>Treatment, e</td><td>StackedClassifier<br>(HistGradientBoostingClassifier, LogisticRegression)</td><td>LogisticRegression C: [0.0001, 0.001, 0.01, 0.1, 1]<br>HistGradientBoostingClassifier learning rate: [0.0001, 0.001, 0.01, 0.1, 1]<br>HistGradientBoostingClassifier max leaf nodes: 10</td></tr></tbody></table>

<footnote n="6">⁶Original R code available at https://github.com/vdorie/aciccomp/tree/master/2016 to generate 77 simulations settings.</footnote>

<footnote n="7">⁷Using the scaling part of the data, from github.com/IBM-HRL-MLHLS/IBM-Causal-Inference-Benchmarking-Framework</footnote>

<footnote n="8">⁸We obtained the dataset from https://github.com/AMLab-Amsterdam/CEVAE/tree/master/datasets/TWINS</footnote>
## G.4 Additional Results

**Definition of the Kendall's tau, $\kappa$** The Kendall's tau is a widely used statistics to measure the rank correlation between two set of observations. It measures the number of concordant pairs minus the discordant pairs normalized by the total number of pairs. It takes values in the [-1, 1] range.

$$ \kappa = \frac{(\text{number of concordant pairs}) - (\text{number of discordant pairs})}{(\text{number of pairs})} \qquad (19) $$

**Values of relative $\kappa(l, \tau$-risk) compared to the mean over all metrics Kendall's as shown in the boxplots of Figure 3**

**Figure 13 - Results measured in relative Kendall's for feasible and semi-oracle risks** Because of extreme propensity scores in the denominator and bayes error residuals in the numerator, the semi-oracle $U$-risk has poor performances at bad overlap. Estimating these propensity scores in the is feasible $U$-risk reduces the variance since clipping is performed.

**Figure 14 - Results measured in absolute Kendall's**

**Figure 15 - Results measured as distance to the oracle tau-risk** To see practical gain in term of $\tau$-risk, we plot the results as the normalized distance between the estimator selected by the oracle $\tau$-risk and the estimator selected by each causal metric.

Then, $\tilde{R}$-risk* is more efficient than all other metrics. The gain are substantial for every datasets.

**Figure 16 - Stacked models for the nuisances is more efficient** For each metrics the benefit of using a stacked model of linear and boosting estimators for nuisances compared to a linear model. The evaluation measure is Kendall's tau relative to the oracle $R$-risk* to have a stable reference between experiments.

<figure><img src="image_17.png" /><figcaption>**Figure 13: The R-risk is the best metric:** Relative Kendall's $\tau$ agreement with $\tau$-risk. Strong and Weak overlap correspond to the first and last tertiles of the overlap distribution measured with Normalized Total Variation eq. 16.</figcaption></figure>
<table><caption>**Table 4:** Values of relative $κ(l, τ-risk)$ compared to the mean over all metrics Kendall's as shown in the boxplots of Figure 3</caption><thead><tr><th rowspan="2">Metric</th><th rowspan="2">Dataset</th><th colspan="2">Strong Overlap</th><th colspan="2">Weak Overlap</th></tr><tr><th>Median</th><th>IQR</th><th>Median</th><th>IQR</th></tr></thead><tbody><tr><td rowspan="4">μ̂-risk</td><td>Twins (N=11 984)</td><td>-0.32</td><td>0.12</td><td>-0.19</td><td>0.12</td></tr><tr><td>ACIC 2016 (N=4 802)</td><td>-0.03</td><td>0.13</td><td>0.11</td><td>0.19</td></tr><tr><td>Caussim (N=5 000)</td><td>-0.40</td><td>0.55</td><td>-0.16</td><td>0.31</td></tr><tr><td>ACIC 2018 (N=5 000)</td><td>0.00</td><td>0.30</td><td>0.01</td><td>0.40</td></tr><tr><td rowspan="4">μ̂-risk<sub>IPW</sub></td><td>Twins (N= 11 984)</td><td>-0.31</td><td>0.13</td><td>-0.17</td><td>0.12</td></tr><tr><td>ACIC 2016 (N=4 802)</td><td>-0.02</td><td>0.13</td><td>0.11</td><td>0.19</td></tr><tr><td>Caussim (N=5 000)</td><td>-0.34</td><td>0.50</td><td>0.09</td><td>0.31</td></tr><tr><td>ACIC 2018 (N=5 000)</td><td>0.00</td><td>0.30</td><td>-0.01</td><td>0.43</td></tr><tr><td rowspan="4">μ̂-risk<sup>*</sup><sub>IPW</sub></td><td>Twins (N= 11 984)</td><td>-0.32</td><td>0.13</td><td>-0.17</td><td>0.13</td></tr><tr><td>ACIC 2016 (N=4 802)</td><td>-0.02</td><td>0.13</td><td>0.11</td><td>0.21</td></tr><tr><td>Caussim (N=5 000)</td><td>-0.33</td><td>0.54</td><td>0.26</td><td>0.27</td></tr><tr><td>Twins (N= 11 984)</td><td>0.13</td><td>0.12</td><td>0.27</td><td>0.12</td></tr><tr><td rowspan="4">τ̂-risk<sub>IPW</sub></td><td>ACIC 2016 (N=4 802)</td><td>-0.07</td><td>0.18</td><td>0.05</td><td>0.31</td></tr><tr><td>Caussim (N=5 000)</td><td>-0.19</td><td>0.43</td><td>-0.14</td><td>0.18</td></tr><tr><td>ACIC 2018 (N=5 000)</td><td>-0.16</td><td>0.40</td><td>-0.11</td><td>0.66</td></tr><tr><td>Twins (N= 11 984)</td><td>0.12</td><td>0.14</td><td>0.20</td><td>0.16</td></tr><tr><td rowspan="4">τ̂-risk<sup>*</sup><sub>IPW</sub></td><td>ACIC 2016 (N=4 802)</td><td>-0.03</td><td>0.16</td><td>-0.09</td><td>0.43</td></tr><tr><td>Caussim (N=5 000)</td><td>-0.15</td><td>0.46</td><td>-0.17</td><td>0.19</td></tr><tr><td>Twins (N= 11 984)</td><td>0.13</td><td>0.12</td><td>0.02</td><td>0.25</td></tr><tr><td>ACIC 2016 (N=4 802)</td><td>0.04</td><td>0.11</td><td>0.11</td><td>0.26</td></tr><tr><td rowspan="4">Û-risk</td><td>Caussim (N=5 000)</td><td>0.04</td><td>0.43</td><td>-0.04</td><td>0.17</td></tr><tr><td>ACIC 2018 (N=5 000)</td><td>0.12</td><td>0.26</td><td>-0.02</td><td>0.50</td></tr><tr><td>Twins (N= 11 984)</td><td>0.25</td><td>0.08</td><td>-0.41</td><td>0.45</td></tr><tr><td>ACIC 2016 (N=4 802)</td><td>0.08</td><td>0.13</td><td>-0.59</td><td>0.57</td></tr><tr><td rowspan="4">Û-risk<sup>*</sup></td><td>Caussim (N=5 000)</td><td>0.46</td><td>0.12</td><td>0.02</td><td>0.44</td></tr><tr><td>Twins (N= 11 984)</td><td>0.15</td><td>0.10</td><td>0.25</td><td>0.18</td></tr><tr><td>ACIC 2016 (N=4 802)</td><td>0.07</td><td>0.12</td><td>0.22</td><td>0.15</td></tr><tr><td>Caussim (N=5 000)</td><td>0.34</td><td>0.26</td><td>0.13</td><td>0.21</td></tr><tr><td rowspan="4">R̂-risk</td><td>ACIC 2018 (N=5 000)</td><td>0.13</td><td>0.27</td><td>0.21</td><td>0.47</td></tr><tr><td>Twins (N= 11 984)</td><td>0.25</td><td>0.10</td><td>0.32</td><td>0.15</td></tr><tr><td>ACIC 2016 (N=4 802)</td><td>0.12</td><td>0.12</td><td>0.25</td><td>0.15</td></tr><tr><td>Caussim (N=5 000)</td><td>0.47</td><td>0.11</td><td>0.16</td><td>0.14</td></tr></tbody></table>
<figure><img src="image_18.png" /><figcaption>**Figure 14:** Agreement with τ-risk ranking of methods function of overlap violation. The lines represent medians, estimated with a lowess. The transparent bands denote the 5% and 95% confidence intervals.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 15:** Metric performances by normalized tau-risk distance to the best method selected with τ-risk. All nui- sances are learned with the same estimator stacking gradient boosting and ridge regression. Dotted and plain lines corresponds to 60% lowess quantile estimates. This choice of quantile allows to see better the oracle metrics lines for which outliers with a value of 0 distord the curves.</figcaption></figure>
Thus, we do not include in this analysis the ACIC 2018 dataset since *R*-risk* is not available due to the lack of the true propensity score.

**Figure 17** *Low population overlap hinders model selection for all metrics*

**Figure 18** - Stacked models for the nuisances is more efficient For each metrics the benefit of using a stacked model of linear and boosting estimators for nuisances compared to a linear model. The evaluation measure is Kendall’s tau relative to the oracle *R*-risk* to have a stable reference between experiments. Thus, we do not include in this analysis the ACIC 2018 dataset since *R*-risk* is not available due to the lack of the true propensity score.

**Figure 19** - Flexible models are performant in recovering nuisances even in linear setups

**Selecting different seeds and parameters is crucial to draw conclusions** One strength of our study is the various number of different simulated and semi-simulated datasets. We are convinced that the usual practice of using only a small number of generation processes does not allow to draw statistically significant conclusions.

Figure 20 illustrate the dependence of the results on the generation process for caussim simulations. We highlighted the different trajectories induced by three different seeds for data generation and three different treatment ratio instead of 1000 different seeds. The result curves are relatively stable from one setup to another for *R*-risk, but vary strongly for *μ*-risk and *μ*-riskIPW.
(a) Caussim

<figure><img src="image_20.png" /></figure>

(b) ACIC 2016

<figure><img src="image_21.png" /></figure>

(c) Twins

<figure><img src="image_22.png" /><figcaption>**Figure 16:** Results are similar between the Shared nuisances/candidate set and the Separated nuisances set procedure. The experience has not been run on the full metrics for Caussim due to computation costs.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 17:** Low population overlap hinders causal model selection for all metrics: Kendall's τ agreement with τ-risk. Strong, medium and Weak overlap correspond to the tertiles of the overlap distribution measured with Normalized Total Variation eq. 16.</figcaption></figure>
(a) Twins

<figure><img src="image_24.png" /></figure>

(b) Twins downsampled

<figure><img src="image_25.png" /></figure>

(c) Caussim

<figure><img src="image_26.png" /></figure>

(d) ACIC 2016

<figure><img src="image_27.png" /></figure>

**Figure 18:** Learning the nuisances with *stacked models* (linear and gradient boosting) is important for successful model selection with R-risk. For Twins dataset, there is no improvement for *stacked models* compared to *linear models* because of the linearity of the propensity model.

<figure><img src="image_28.png" /></figure>

**Figure 19:** Flexible models are performant in recovering nuisances in the downsampled Twins dataset. The propensity score is linear in this setup, making it particularly challenging for flexible models compared to linear methods.
<figure><img src="image_29.png" /><figcaption>**Figure 21:** a) For CATE, a train/test ratio of 0.9/0.1 appears a good trade-off. b) For ATE, there is a small signal pointing also to 0.9/0.1 (K=10). For ATE, experiences on 10 replications of all 78 instances of the ACIC 2016 data.</figcaption></figure>

# H Data split choices

## H.1 Use 90% of the data to estimate outcome models, 10% to select them

The analyst faces a compromise: given a finite data sample, should she allocate more data to estimate the outcome model, thus improving the quality of the outcome model but leaving little data for model selection. Or, she could choose a bigger test set for model selection and effect estimation. For causal model selection, there is no established practice (as reviewed in H.2).

We investigate such tradeoff varying the ratio between train and test data size. For this, we first split out 30% of the data as a holdout set $V$ on which we use the oracle response functions to derive silver-standard estimates of causal quantities. We then use the standard estimation procedure on the remaining 70% of the data, splitting it into train $T$ and test $S$ of varying sizes. We finally measure the error between this estimate and the silver standard.

We consider two different analytic goals: estimating a average treatment effect –a single number used for policy making– and a CATE –a full model of the treatment effect as a function of covariates $X$. Given that the latter is a much more complex object than the former, the optimal train/test ratio might vary. To measure errors, we use for the ATE the relative absolute ATE bias between the ATE computed with the selected outcome model on the test set, and the true ATE as evaluated on the holdout set $V$. For the CATE, we compare the $\tau$-risk of the best selected model applied on the holdout set $V$. We explore this trade-off for the ACIC 2016 dataset and the R-risk.

Figure 21 shows that a train/test ratio of 0.9/0.1 (K=10) or 0.8/0.2 (K=5) appears best to estimate CATE and ATE.
## H.2 Heterogeneity in practices for data split

Splitting the data is common when using machine learning for causal inference, but practices vary widely in terms of the fraction of data to allocate to train models, outcomes and nuisances, and to evaluate them.

Before even model selection, data splitting is often required for estimation of the treatment effect, ATE or CATE, for instance to compute the nuisances required to optimize the outcome model (as the R-risk, Definition 5). The most frequent choice is use 80% of the data to fit the models, and 20% to evaluate them. For instance, for CATE estimation, the R-learner has been introduced using K-folds with $K = 5$ and $K = 10$: 80% of the data (4 folds) to train the nuisances and the remaining fold to minimize the corresponding R-loss⁵². Yet, it has been implemented with $K=5$ in causallib⁷⁵ or $K=3$ in econML⁸. Likewise, for ATE estimation, Chernozhukov et al. [14] introduce doubly-robust machine learning, recommending $K=5$ based on an empirical comparison $K=2$. However, subsequent works use doubly robust ML with varying choices of $K$: Loiseau et al. [42] use $K=3$, Gao et al. [21] use $K=2$. In the econML implementation, $K$ is set to 3⁸. Naimi et al. [48] evaluate various machine-learning approaches –including R-learners– using $K=5$ and 10, drawing inspiration from the TMLE literature which sets $K=5$ in the TMLE package²⁴.

Causal model selection has been much less discussed. The only study that we are aware of, Schuler et al. [70], use a different data split: a 2-folds train/test procedure, training the nuisances on the first half of the data, and using the second half to estimate the R-risk and select the best treatment effect model.