# Optimal Transport for Counterfactual Estimation:
A Method for Causal Inference

Arthur Charpentiera∗, Emmanuel Flachaireb & Ewen Gallicb
a Université du Québec à Montréal (UQAM), Montréal (Québec), Canada
b Aix Marseille Univ, CNRS, AMSE, Marseille, France

∗Corresponding author<EMAIL_ADDRESS>
###### Abstract

Many problems ask a question that can be formulated as a causal question: what
would have happened if…? For example, would the person have had surgery if he
or she had been Black? To address this kind of questions, calculating an
average treatment effect (ATE) is often uninformative, because one would like
to know how much impact a variable (such as skin color) has on a specific
individual, characterized by certain covariates. Trying to calculate a
conditional ATE (CATE) seems more appropriate. In causal inference, the
propensity score approach assumes that the treatment is influenced by
$\boldsymbol{x}$, a collection of covariates. Here, we will have the dual
view: doing an intervention, or changing the treatment (even just
hypothetically, in a thought experiment, for example by asking what would have
happened if a person had been Black) can have an impact on the values of
$\boldsymbol{x}$. We will see here that optimal transport allows us to change
certain characteristics that are influenced by the variable we are trying to
quantify the effect of. We propose here a mutatis mutandis version of the
CATE, which will be done simply in dimension one by saying that the CATE must
be computed relative to a level of probability, associated to the proportion
of ${x}$ (a single covariate) in the control population, and by looking for
the equivalent quantile in the test population. In higher dimension, it will
be necessary to go through transport, and an application will be proposed on
the impact of some variables on the probability of having an unnatural birth
(the fact that the mother smokes, or that the mother is Black).

Keywords Causality; Conditional Average Treatment Effects (CATE);
Counterfactual; Mutatis Mutandis; Optimal Transport; Quantiles

##  1 Introduction

###  1.1 From intervention to counterfactuals

In Pearl and Mackenzie (2018), a “ladder of causation” is introduced, to
describe the three levels of causal reasoning. The first level, named
“association”, discusses associations (not to use the word “correlation”)
between variables. Questions such as “is variable $X$ associated with variable
$Y$?” can be answered at this level. Econometric models are usually simply
based on such associations. The second level is labelled “intervention”.
Reasoning on this level answers questions of the form “if I make the
intervention $T$, how will this affect the level of the outcome $Y$?” For
example, the question “would a patient heal faster at home or at the hospital,
after some surgery?” is a standard question on this second level of the ladder
of causation. This kind of reasoning invokes causality and can be used to
investigate more questions than the reasoning of the first level. The third
level of the “ladder of causation” is labelled “counterfactuals” and involves
answering questions which ask what might have been, had circumstances been
different. Counterfactual modeling implies that, to each individual in the
control space, described through variables $\boldsymbol{x}$ and $y$, we will
associate a counterfactual version of that individual in the hypothetical
space. More formally, we will use notations of causal inference to answer
counterfactual questions, such as “would that person have had surgery if she
had been Afro-American?”

###  1.2 Causal inference framework

Consider, as in Rubin (1974) or Hernán and Robins (2010), the following
framework: let $t$ denote some binary treatment, $t\in\\{0,1\\}$, with
respectively, the control and the treatment. Let $\boldsymbol{x}$ be some
covariates, $y$ the observed outcome, with $y_{T\leftarrow 1}^{\star}$ and
$y_{T\leftarrow 0}^{\star}$ the potential outcomes (also denoted $y(1)$ and
$y(0)$ in Imbens and Rubin (2015) or Imai (2018), or $y^{1}$ and $y^{0}$ in
Morgan and Winship (2014) or Cunningham (2021), even $y_{t=1}$ and $y_{t=0}$
in Pearl and Mackenzie (2018)), realized either under treatment condition
($t=1$) or under control condition ($t=0$). Note that the observed outcome is
$y=y_{T\leftarrow t}^{\star}$, or $y=t\cdot y_{T\leftarrow
1}^{\star}+(1-t)\cdot y_{T\leftarrow 0}^{\star}$. An illustration is reported
in Table 1.

| Treatment | Outcome | Age | Gender | Height | Weight
---|---|---|---|---|---|---
| $t_{i}$ | $y_{i}$ | $y_{i,T\leftarrow 1}^{\star}$ | $y_{i,T\leftarrow o}^{\star}$ | $x_{1,i}$ | $x_{2,i}$ | $x_{3,i}$ | $x_{4,i}$
1 | 1 | 121 | 121 | ? | 37 | F | 160 | 56
2 | 0 | 109 | ? | 109 | 28 | F | 156 | 54
3 | 1 | 162 | 162 | ? | 53 | M | 190 | 87
Table 1: Potential outcome framework of causal inference, with one binary
treatment $t_{i}$, the observed outcome variable $y_{i}$ and the two potential
outcomes $y_{i,T\leftarrow 1}^{\star}$ and $y_{i,T\leftarrow 0}^{\star}$, as
well as some covariates $\boldsymbol{x}_{i}$. One of the two potential
outcomes is observed, and the other is missing, indicated by the question mark
in the table.

We will use the term “treatment” (and letter $t$) even if interventions are
not possible, so it is no per se a “treatment”. In this article, we try to
answer a hypothetical question, like most questions asked at the third level
of the “ladder of causality”. For instance, in a context of quantifying
discrimination, the “treatment” will denote the sensitive attribute, as in
Charpentier (2023), such as the race of an individual, e.g., “what would have
been the outcome if that person had been Afro-American?” Since our approach
proposes an improvement on the metrics used in causal inference literature, we
will use similar notations.

There will be a significant impact of treatment $t$ on $y$ if
$y^{\star}_{T\leftarrow 0}\neq y^{\star}_{T\leftarrow 1}$. More specifically,
the causal effect for individual $i$ is ${\tau_{i}=y^{\star}_{i,T\leftarrow
1}-y^{\star}_{i,T\leftarrow 0}}$. The average treatment effect (ATE) can the
be defined as follows:

$\tau=\text{ATE}=\mathbb{E}\big{[}Y^{\star}_{i,T\leftarrow
1}-Y^{\star}_{i,T\leftarrow 0}\big{]}.$

Its empirical counterpart, the sample average treatment effect (SATE) writes:

$\widehat{\tau}=\text{SATE}=\frac{1}{n}\sum_{i=1}^{n}y^{\star}_{i,T\leftarrow
1}-y^{\star}_{i,T\leftarrow 0}.$

Unfortunately, the latter is not directly observable, since one of the two is
always missing, but some techniques can be used to provide some robust
estimate of that quantity (we will present some of them in the next section).

Lastly, in the context of possibly heterogeneous effects, captured through
covariates $\boldsymbol{x}$ (that can be a subset of the entire set of
covariates), the conditional average treatment effect (CATE) is defined as the
functional

$\tau(\boldsymbol{x})=\text{CATE}(\boldsymbol{x})=\mathbb{E}\big{[}Y^{\star}_{T\leftarrow
1}-Y^{\star}_{T\leftarrow 0}\big{|}\boldsymbol{X}=\boldsymbol{x}\big{]}$

that can be written

$\tau(\boldsymbol{x})=\text{CATE}(\boldsymbol{x})=\mathbb{E}\big{[}Y^{\star}_{T\leftarrow
1}\big{|}\boldsymbol{X}=\boldsymbol{x}\big{]}-\mathbb{E}\big{[}Y^{\star}_{T\leftarrow
0}\big{|}\boldsymbol{X}=\boldsymbol{x}\big{]},$

as introduced in Hahn (1998) and Heckman et al. (1998). More recently, Hitsch
and Misra (2018) used that measure to quantify heterogeneous treatment effects
to evaluate optimal targeting policies, as well as Powers et al. (2018) and
Fan et al. (2022). Wager and Athey (2018), Athey and Wager (2019) and Athey et
al. (2019) suggested to use random forests to estimate this quantity, inspired
by Davis and Heller (2017). See also Künzel et al. (2019) or Hsu et al. (2022)
for additional discussion on that quantity.

A classical assumption is that $(t_{i},y_{i},\boldsymbol{x}_{i})$ is a random
sample of size $n$ from some joint random vector $(T,Y,\boldsymbol{X})$.
Rosenbaum and Rubin (1983) suggested a strong “ignorable treatment assignment”
assumption defined as a conditional independence between
$(Y^{\star}_{T\leftarrow 0},Y^{\star}_{T\leftarrow 1})$ and $T$, conditional
on the covariates $\boldsymbol{X}$.

###  1.3 Agenda

In Section 2, and more specifically in Section 2.1, we will discuss further
the (possible) connection between covariates $\boldsymbol{x}$, treatment $t$
and the outcome $y$. Following our example on discrimination, the treatment
variable $t$ (such as skin color) is an “exogenous variable”, in the sense
that it cannot be influenced either by covariates $\boldsymbol{x}$ or by the
outcome $y$. Using the terminology from directed acyclic graphs (DAGs), $t$
will have no parent, so in a sense, it will be easier to pretend that an
hypothetical intervention on $t$ is possible. In most applications, $t$ will
have an impact on the outcome $y$, but not only. More precisely, it is
possible that $t$ might influence some covariates $\boldsymbol{x}$, and those
covariates can, in turn, impact the outcome $y$. In Section 2.2, we suggest an
extension from the standard ceteris paribus $\text{CATE}(\boldsymbol{x})$
defined as the difference $\mathbb{E}\big{[}Y^{*}_{T\leftarrow
1}\big{|}\boldsymbol{x}\big{]}-\mathbb{E}\big{[}Y^{*}_{T\leftarrow
0}\big{|}\boldsymbol{x}\big{]}$, to some mutatis mutandis
$\text{CATE}(\boldsymbol{x})$ defined as the difference
$\mathbb{E}\big{[}Y^{*}_{T\leftarrow 1}\big{|}\boldsymbol{x}_{T\leftarrow
1}\big{]}-\mathbb{E}\big{[}Y^{*}_{T\leftarrow 0}\big{|}\boldsymbol{x}\big{]}$,
where, if $\boldsymbol{x}$ is considered with respect to the control group,
the counterfactual in the treated population should be based on a different
version of $\boldsymbol{x}$, in the treated space. As discussed in Section
2.3, the classical tool used in econometrics is the propensity score, based on
$\mathbb{P}[T=1|\boldsymbol{X}=\boldsymbol{x}]$, that is usually considered to
take into account the association that exists between the treatment and the
covariates. At the second stage of the “ladder of causation” –the
intervention– we consider the fact that $\boldsymbol{x}$ might influence $t$.
When answering the question “would a patient heal faster at home or at the
hospital, after some surgery?”, it might be relevant to assume that the
propensity score can be used to correct for the bias we have in the data,
since some patient have been healing at the hospital, not by choice, but
because of some $\boldsymbol{x}$. At the third stage of the ladder –the
counterfactuals– some sort of dual version should be considered, since $t$ is
not influenced by $\boldsymbol{x}$, quite the opposite: some $\boldsymbol{x}$
might be influenced by $t$. A simple toy example, based on a Gaussian
structural equation model (SEM), is presented in Section 2.4, while in Section
2.5, we briefly present real data that we will use in the next sections to
illustrate various algorithms, based on births in the United States. The
variable of interest $y$ is a binary variable, indicating whether a birth was
natural, or not. The covariates $\boldsymbol{x}$ considered here will be the
weight of the newborn, and the weight gain of the mother. And various
“treatments” are considered: whether the mother is Afro-American, or not;
whether the mother is a smoker, or not; whether the baby is a girl, or not
(results for the last two are reported in Appendix A.2).

In Section 3, we will focus on the case where only one covariate $x$ is
considered. We will start with classical matching techniques in Section 3.1,
used to match each point in $(y_{i},x_{i},t_{i}=0)$ –in the control group–
with another one in $(y_{j},x_{j},t_{j}=1)$ –in the treated group– when the
two groups have the same size. In Section 3.2, we will suggest on “optimal”
matching algorithm, to associate individual $i$ (in the control group) to $j$
(in the treated group), that we will denote $j_{i}^{\star}$. Then, in Section
3.3, we will discuss the case where the two groups have different sizes, that
will be called optimal “coupling”. In Section 3.4, we will define an
estimator, the mutatis mutandis CATE,
$\widehat{m}_{1}\big{(}\widehat{\mathcal{T}}(x)\big{)}-\widehat{m}_{0}\big{(}x\big{)}$,
where $\widehat{\mathcal{T}}(x)=\widehat{F}_{1}^{-1}\circ\widehat{F}_{0}(x)$,
with $\widehat{F}_{0}$ and $\widehat{F}_{1}$ denoting the empirical
distribution functions of $x$ conditional on $t=0$ and $t=1$, respectively. We
will use quantiles to optimally “transport” $\boldsymbol{x}$’s from the
control group to the treated group, formally through the $\mathcal{T}$
mapping. Finally, in Section 3.5, we will illustrate this on probability to
have a non-natural baby delivery, on our dataset.

In Section 4, we will extend our previous approach to the case where several
covariates $\boldsymbol{x}$ are considered. Formally, we will use optimal
transport techniques to get a proper counterfactual of $\boldsymbol{x}$, not
in the control group, but in the treated group. In Section 4.1, we will define
the optimal transport problem for any number of dimensions and then, in
Section 4.2, we will explain how to optimally associate each observation
$\boldsymbol{x}_{i}$ in the control group (when $t=0$) with a single
counterfactual observation $\boldsymbol{x}_{j}$ in the treated group (when
$t=1$) when the two groups have the same size. This can be related to the
Gaussian SEM discussed in Section 2.4. In Section 4.3, we will see the
extension when the two groups have different sizes. Unfortunately, those
approach do not provide an explicit mapping $\mathcal{T}$, but simply a
matching of a single individual $\boldsymbol{x}_{i}$ (in the control group) to
a weighted sum of multiple $\boldsymbol{x}_{j}$ (in the treated group). As we
will see in Section 4.4, it will be possible to get explicit formulation for
the mapping $\mathcal{T}$ (from the space of covariates in the control group
to the space of covariates in the treated group) when we assume that
$\boldsymbol{X}$ conditional on $T$ has Gaussian distributions. In Section
4.5, those techniques will be further discussed in the context of the
application to non-natural birth111See https://github.com/3wen/counterfactual-
estimation-optimal-transport for more details..

##  2 Ceteris Paribus vs. Mutatis Mutandis

Before introducing another concept of CATE, we will formalize a little bit
more the connections between the “treatment” $t$, the outcome $y$ and the
covariates $\boldsymbol{x}$.

###  2.1 Exogeneity, endogeneity and causal graphs

As discussed earlier, when presenting the second stage of the “ladder of
causation”, $t$ is a treatment. For example, in epidemiology, $t$ may be a
treatment given to patients, possibly resulting from an intervention. At the
third level, the treatment would be more a thought experiment (the
“gedankenexperiment” in Mach (1893)), to answer a question such as “what if
$t$ had taken another value?”, without being able to make an experiment.
Chisholm (1946) introduced the idea of “contrary-to-fact conditional”, coined
as “counterfactual” in Goodman (1947). A classical example would be when
$t\in\\{\text{smoker},\text{non-smoker}\\}$, since it is not ethically
possible to force someone to smoke, but it can also be used on inherent
variables, such as the gender or the race of a person, that cannot be changed
in a real experiment, to quantify possible discrimination.

Covariates $\boldsymbol{x}$ are available variables that have an impact on the
outcome $y$. It is necessary here to distinguish two kinds of covariates, with
variables that are influenced by the value of $t$, that might be seen as
”endogenous”, and those that are not influenced by the value of $t$, that
might be seen as “exogenous”. For example, the weight of the baby $x$ is an
endogenous variable with respect to the variable indicating whether the mother
is a smoker or not. Using a terminology used on causal graphs, “endogeneous”
covariates $x$ are mediator variables (between $t$ and $y$), while
“exogeneous” ones are variables colliding with $t$ on $y$, sometimes called
collider variables (see Figure 1).

The Markov assumption, on causal networks, states that each variable is
conditionally independent of its non-descendants, given its parents. In Figure
1, in the ‘cofounder’ case (with the fork $t\to x$ and $t\to y$), and in the
‘mediator’ case (with the chain $t\to x\to y$), $y$ is independent of $t$,
conditional on $x$. But in the “colider” case (with $x\to y$ and $t\to y$),
while $x$ and $t$ are independent, they become conditionally dependent,
conditional on $y$. We will not discuss here the construction of the causal
graphs, that is supposed to be given (see, e.g., Vowels et al. (2022) for a
survey on techniques used to discover causal structures).

$y$$\boldsymbol{x}^{m}$$t$$\boldsymbol{x}^{c}$treatmentoutcome mediator
(endogeneous) colider (exogeneous) $\boldsymbol{x}^{p}$confounding(noise,
proxy) Figure 1: Distinction of covariates, with confounding variables that
will not influence $y$ on the top right, and two sets of explanatory variables
that will influence $y$, that are influenced, or not, by “treatment” $t$, with
mediators and colliders, at the bottom left.

###  2.2 Impact of a treatment $t$ on $y$ and $\boldsymbol{x}$, and CATE

$y$$\boldsymbol{x}^{m}$$t$$\boldsymbol{x}^{c}$$y^{\prime}$$\boldsymbol{x}^{\prime
m}$$t^{\prime}$$\boldsymbol{x}^{c}$$y^{\prime}=h(t^{\prime},\boldsymbol{x}^{\prime
m},\boldsymbol{x}^{c})$DAGintervention on the treatment Figure 2: A causal
graph on the left, and the impact of an intervention on the treatment $t$ on
the right.

Consider some treatment $t$. Let $\boldsymbol{x}^{m}$ denote the set of
mediator variables and $\boldsymbol{x}^{c}$ denote the set of collider
variables, as in Figure 2. Following the SEM terminology used in causal
inference, consider data generated according to the equations on the left
below (real world), prior to intervention on $t$. The right hand equations
describe the data generating process with an intervention on $t$ (denoted
$do(t)$ in Pearl and Mackenzie (2018)):

real world | with intervention ($do(t)$)
---|---
$\displaystyle{\begin{cases}T=h_{t}(U_{t})\\\ \boldsymbol{X}^{m}=h_{m}(T,\boldsymbol{U}_{m})\\\ \boldsymbol{X}^{c}=h_{c}(\boldsymbol{U}_{c})\\\ Y=h_{y}(T,\boldsymbol{X}^{m},\boldsymbol{X}^{c},U_{y})\\\ \end{cases}}$ | $\displaystyle{\begin{cases}T=t\\\ \boldsymbol{X}^{m}_{T\leftarrow t}=h_{m}(t,\boldsymbol{U}_{m})\\\ \boldsymbol{X}^{c}=h_{c}(\boldsymbol{U}_{c})\\\ Y_{T\leftarrow t}=h_{y}(t,\boldsymbol{X}^{m}_{T\leftarrow t},\boldsymbol{X}^{c},U_{y})\\\ \end{cases}}$

Consider some independent noise variables
$\\{U_{t},\boldsymbol{U}_{m},\boldsymbol{U}_{c},U_{y}\\}$ (that can be assumed
to be centered Gaussian to be close to the econometric literature). In the
“real world”, $T$ is a function of $U_{t}$, and $U_{t}$ only, through some
$h_{t}:\mathbb{R}\to\\{0,1\\}$ function,
$h_{t}(u)=\boldsymbol{1}(u>\text{threshold})$. Then we have two possible
explanatory variables: mediator (endogenous) and collider (exogenous). If
$\boldsymbol{X}^{c}$ are functions of the noise $\boldsymbol{U}_{c}$ only
(through function $h_{c}$), $\boldsymbol{X}^{m}$ are functions of the noise
$\boldsymbol{U}_{m}$ and the treatment $T$ (through function $h_{c}$). And
finally, the outcome $Y$ is function of $\boldsymbol{X}^{c}$ and
$\boldsymbol{X}^{m}$, also possibly $T$, and some idiosyncratic noise $U_{y}$.

In a ceteris paribus approach, $\text{CATE}(x)$ is equal to
$\mathbb{E}\big{[}Y^{*}_{T\leftarrow
1}\big{|}{x}\big{]}-\mathbb{E}\big{[}Y^{*}_{T\leftarrow 0}\big{|}{x}\big{]}$.
In a mutatis mutandis version, we should not consider $x$, but a version of
$x$ that should be influenced by the treatment $t$, denoted ${x}_{T\leftarrow
1}$. In a general setting, we have the following definition:

######  Definition 2.1.

The mutatis mutandis CATE is

$\text{CATE}(\boldsymbol{x})=\mathbb{E}\big{[}Y^{*}_{T\leftarrow
1}\big{|}\boldsymbol{x}_{T\leftarrow
1}\big{]}-\mathbb{E}\big{[}Y^{*}_{T\leftarrow 0}\big{|}\boldsymbol{x}\big{]}$

(we might denote $\boldsymbol{x}_{T\leftarrow 0}$ instead of $x$ to avoid
confusion for the second term).

More specifically, when we ask the question “what would have been the
probability to have a non-natural delivery for a baby with weight $x$ if the
mother had been smoking?”, we have to take into account the fact that if the
mother had been smoking, the weight of the baby would have been impacted. The
original weight $x$, associated with a non-Black mother, would become
$\boldsymbol{x}_{T\leftarrow 1}$ (instead of $x$) if we seek a counterfactual
version of $x$ in the treated population.

###  2.3 Propensity score weighting

The classical approach in causal inference is based on the idea that $T$ is
not really exogenous, and can be influenced by $\boldsymbol{x}$. Therefore,
the average treatment effect $\text{ATE}=\mathbb{E}[Y^{\star}_{T\leftarrow
1}-Y^{\star}_{T\leftarrow 0}]$, that can be written

$\text{ATE}=\mathbb{E}\left[\frac{TY}{p(\boldsymbol{X})}-\frac{(1-T)Y}{1-p(\boldsymbol{X})}\right]$

would be estimated by

$\text{SATE}=\frac{1}{n}\sum_{i=1}^{n}\frac{t_{i}y_{i}}{\widehat{p}(\boldsymbol{x}_{i})}-\frac{(1-t_{i})y_{i}}{1-\widehat{p}(\boldsymbol{x}_{i})},$

where $p(\boldsymbol{x})$ is a “propensity score” defined as
$p(\boldsymbol{x})=\mathbb{P}[T=1|\boldsymbol{X}=\boldsymbol{x}]$, that can be
estimated using, for instance, a logistic regression

$\widehat{p}(\boldsymbol{x})=\frac{\exp[\boldsymbol{x}^{\top}\widehat{\boldsymbol{\beta}}]}{1+\exp[\boldsymbol{x}^{\top}\widehat{\boldsymbol{\beta}}]}.$

Thus, the SATE can be seen as the difference between two weighted averages of
$y_{i}$’s. As discussed in Abrevaya et al. (2015), it can be used to estimate
$\text{CATE}(x)$, on a subset of features, with a local estimate of the
average

$\text{CATE}(x)=\frac{1}{\sum
K_{h}(x_{i}-x)}\sum\left(\frac{t_{i}y_{i}}{\widehat{p}(\boldsymbol{x}_{i})}-\frac{(1-t_{i})y_{i}}{1-\widehat{p}(\boldsymbol{x}_{i})}\right)K_{h}(x_{i}-x),$

using some kernel function $K_{h}$. A $k$-nearest neighbors estimate can also
be considered:

$\text{CATE}(x)=\frac{1}{k}\sum_{i\in\mathcal{V}_{k}(x)}\left(\frac{t_{i}y_{i}}{\widehat{p}(\boldsymbol{x}_{i})}-\frac{(1-t_{i})y_{i}}{1-\widehat{p}(\boldsymbol{x}_{i})}\right),$

where $i\in\mathcal{V}_{k}(x)$ when $x_{i}$ is among the $k$-nearest neighbors
of $x$. If $y$ is binary (as the example we will use later on), the ATE is a
difference between two probabilities, and logtistic regressions can be used to
properly estimate $\mathbb{E}[Y^{\star}_{T\leftarrow
t}|\boldsymbol{X}=\boldsymbol{x}]$, with weights in the regressions, that
would be either the inverse of $1-\widehat{p}(\boldsymbol{x}_{i})$ if
$t_{i}=0$ or the inverse of $\widehat{p}(\boldsymbol{x}_{i})$ if $t_{i}=1$, as
in Li et al. (2018).

###  2.4 A toy (Gaussian) example

To illustrate our approach, as an alternative to the use of a propensity
score, consider the following toy example, with three explanatory variables,
two endogenous (and correlated) ones, and an exogenous one, with some linear
model (a Gaussian structural equation model, SEM):

$\displaystyle{\begin{cases}T=\boldsymbol{1}(U_{t}<0),\leavevmode\nobreak\
U_{t}\sim\mathcal{N}(0,1)\\\
\boldsymbol{X}^{m}=\boldsymbol{\mu}_{T}+\boldsymbol{\Sigma}_{T}^{1/2}\boldsymbol{U}_{m},\leavevmode\nobreak\
\boldsymbol{U}_{m}\sim\mathcal{N}(\boldsymbol{0},\mathbb{I})\\\
{X}^{c}=\mu+\sigma{U}_{c},\leavevmode\nobreak\ U_{c}\sim\mathcal{N}(0,1)\\\
Y=\alpha+(\boldsymbol{\beta}_{m},\beta_{c})(\boldsymbol{X}^{m},X^{c})^{\top}+\gamma
T+U_{y},\leavevmode\nobreak\ U_{y}\sim\mathcal{N}(0,1)\\\ \end{cases}}$ (1)

where all the noises $(U_{t},\boldsymbol{U}_{m},U_{c},U_{y})$ are assumed to
be centered, and independent. Here $\boldsymbol{\Sigma}_{0}^{1/2}$ is Cholesky
decomposition of $\boldsymbol{\Sigma}_{0}$, so that $\boldsymbol{X}^{m}$
conditional on $T=t$ has distribution
$\mathcal{N}(\boldsymbol{\mu}_{t},\boldsymbol{\Sigma}_{t})$. Treatment $T$ is
a binary variable, well-balanced since $\mathbb{P}(T=0)=\mathbb{P}(T=1)$.
Conditional on $T=t$, the mediator (endogenous) variables $\boldsymbol{X}^{m}$
have a Gaussian distribution, with mean $\boldsymbol{\mu}_{t}$ and variance
matrix $\boldsymbol{\Sigma}_{t}$. A collider variable $X^{c}$ is supposed to
be independent of the other ones. And finally, $Y$ is a Gaussian variable
where the average is a linear combination of $\boldsymbol{X}^{m}$ and $X^{c}$,
plus $\gamma$ when $T=1$. In Figure 3, the left-hand panel shows a scatter
plot of $\boldsymbol{x}^{m}=(x^{m}_{1},x^{m}_{2})$ with blue points when $t=0$
and red points when $t=1$. The right-hand panel shows $(x^{m}_{1},t)$ on a
scatter plot, with the two conditional densities, as well as the logistic
regression of $t$ against $x_{1}^{m}$ (that could be seen as the propensity
score).

Figure 3: Scatter plot of $\boldsymbol{x}^{m}=(x^{m}_{1},x^{m}_{2})$ with blue
points when $t=0$, and red points when $t=1$, on the left, and the logistic
regression of $t$ against $x_{1}^{m}$ on the right. Toy dataset generated from
Equation (1).

The two interventions yield

$do(T=0)$ | $do(T=1)$
---|---
$\displaystyle{\begin{cases}T\leftarrow 0\\\ \boldsymbol{X}^{m}=\boldsymbol{\mu}_{0}+\boldsymbol{\Sigma}_{0}^{1/2}\boldsymbol{U}_{m}\\\ {X}^{c}=\mu+\sigma{U}_{c}\\\ Y=\alpha+(\boldsymbol{\beta}_{m},\beta_{c})(\boldsymbol{X}^{m},X^{c})^{\top}+U_{y}\\\ \end{cases}}$ | $\displaystyle{\begin{cases}T\leftarrow 1\\\ \boldsymbol{X}^{m}=\boldsymbol{\mu}_{1}+\boldsymbol{\Sigma}_{1}^{1/2}\boldsymbol{U}_{m}^{\prime}\\\ {X}^{c}=\mu+\sigma{U}_{c}^{\prime}\\\ Y=\alpha+(\boldsymbol{\beta}_{m},\beta_{c})(\boldsymbol{X}^{m},X^{c})^{\top}+\gamma+U_{y}^{\prime}\\\ \end{cases}}$

more precisely, in that model with three covariates,
$\boldsymbol{X}^{m}=(X_{1}^{m},X_{2}^{m})$, and since

$\boldsymbol{\Sigma}_{t}=\begin{pmatrix}\sigma_{t1}^{2}&r_{t}\sigma_{t1}\sigma_{t2}\\\
r_{t}\sigma_{t1}\sigma_{t2}&\sigma_{t2}^{2}\\\ \end{pmatrix}\text{ and
}\boldsymbol{\Sigma}_{t}^{1/2}=\begin{pmatrix}\sigma_{t1}&0\\\
\sigma_{t2}r_{t}&\sigma_{t2}\sqrt{1-r_{t}^{2}}\\\ \end{pmatrix}$

we can write

$do(T=0)$ | $do(T=1)$
---|---
$\displaystyle{\begin{cases}T\leftarrow 0\\\ X_{1}^{m}=\mu_{01}+\sigma_{01}U_{1}^{m}\\\ X_{2}^{m}=\mu_{02}+\sigma_{02}(r_{0}U_{1}^{m}+\sqrt{1-r_{0}^{2}}U_{2}^{m})\\\ {X}^{c}=\mu+\sigma{U}_{c}\\\ Y=\alpha+\beta_{1}^{m}X_{1}^{m}+\beta_{2}^{m}X_{2}^{m}+\beta^{c}X^{c}+U_{y}\\\ \end{cases}}$ | $\displaystyle{\begin{cases}T\leftarrow 1\\\ X_{1}^{m}=\mu_{11}+\sigma_{11}U_{1}^{m^{\prime}}\\\ X_{2}^{m}=\mu_{12}+\sigma_{12}(r_{1}U_{1}^{m^{\prime}}+\sqrt{1-r_{1}^{2}}U_{2}^{m^{\prime}})\\\ {X}^{c}=\mu+\sigma{U}_{c}^{\prime}\\\ Y=\alpha+\beta_{1}^{m}X_{1}^{m}+\beta_{2}^{m}X_{2}^{m}+\beta^{c}X^{c}+\gamma+U_{y}^{\prime}\\\ \end{cases}}$

and therefore

$\displaystyle{\begin{cases}Y_{T\leftarrow
0}=\alpha+\beta_{1}^{m}x_{1}+\beta_{2}^{m}\big{(}\mu_{02}+\sigma_{02}(r_{0}\sigma_{01}^{-1}[x_{1}-\mu_{01}]+\sqrt{1-r_{0}^{2}}U_{2}^{m})\big{)}+\beta^{c}(\mu+\sigma{U}_{c})+U_{y}\\\
Y_{T\leftarrow
1}=\alpha+\beta_{1}^{m}x_{1}^{\prime}+\beta_{2}^{m}\big{(}\mu_{12}+\sigma_{12}(r_{1}\sigma_{11}^{-1}[x_{1}^{\prime}-\mu_{11}]+\sqrt{1-r_{1}^{2}}U_{2}^{m^{\prime}})\big{)}+\beta^{c}(\mu+\sigma{U}_{c}^{\prime})+\gamma+U_{y}^{\prime}\\\
\end{cases}}.$

Hence,

$\text{ATE}=\mathbb{E}[Y_{T\leftarrow 1}-Y_{T\leftarrow 0}]=\gamma.$

For conditional average treatment effects,

$\displaystyle{\begin{cases}\mathbb{E}[Y_{T\leftarrow
0}|X_{1}^{m}=x_{1}]=\alpha+\beta_{1}^{m}x_{1}+\beta_{2}^{m}\big{(}\mu_{02}+\sigma_{02}r_{0}\sigma_{01}^{-1}[x_{1}-\mu_{01}]\big{)}+\beta^{c}\mu\\\
\mathbb{E}[Y_{T\leftarrow
1}|X_{1}^{m}=x_{1}^{\prime}]=\alpha+\beta_{1}^{m}x_{1}^{\prime}+\beta_{2}^{m}\big{(}\mu_{12}+\sigma_{12}r_{1}\sigma_{11}^{-1}[x_{1}^{\prime}-\mu_{11}]\big{)}+\beta^{c}\mu+\gamma\\\
\end{cases}}.$

Ceteris paribus, we suppose that $x_{1}^{\prime}=x_{1}$, then

$\text{CATE}_{cp}(x_{1})=\mathbb{E}[Y_{T\leftarrow
1}|X_{1}^{m}=x_{1}]-\mathbb{E}[Y_{T\leftarrow
0}|X_{1}^{m}=x_{1}]=\text{ATE}+\delta x_{1}+\kappa,$

where

$\begin{cases}\kappa=\beta_{2}^{m}\big{(}\mu_{12}+\sigma_{02}r_{0}\sigma_{01}^{-1}\mu_{01}-\sigma_{12}r_{1}\sigma_{11}^{-1}\mu_{11}-\mu_{02}\big{)}\\\
\delta=\beta_{2}^{m}\big{(}\sigma_{12}r_{1}\sigma_{11}^{-1}-\sigma_{02}r_{0}\sigma_{01}^{-1}\big{)}\end{cases}.$

Mutatis mutandis, since $X_{1}^{m}=\mu_{01}+\sigma_{01}U_{1}^{m}$ when $T=0$
while $X_{1}^{m}=\mu_{11}+\sigma_{11}U_{1}^{m}$ when $t=1$, it is legitimate
to consider that $x_{1}^{\prime}=x_{1:T\leftarrow
1}=\mu_{11}+\sigma_{11}(\sigma_{01}^{-1}[x_{1}-\mu_{01}])$. Therefore, mutatis
mutandis,

$\text{CATE}_{mm}(x_{1})=\text{ATE}+\delta^{\prime}x_{1}+\kappa^{\prime},$

where

$\begin{cases}\kappa^{\prime}=\kappa+\beta_{1}^{m}[\mu_{11}-\sigma_{11}\sigma_{01}^{-1}\mu_{01}]\kappa+k\\\
\delta^{\prime}=\delta+\beta_{1}^{m}(\sigma_{11}\sigma_{01}^{-1}-1)=\delta+d\end{cases},$

so that we can also write

$\text{CATE}_{mm}(x_{1})=\text{CATE}_{cp}(x_{1})+\big{(}dx_{1}+k\big{)}.$

In Figure 4, the horizontal orange line is the true average treatment effect
(ATE). The green line is the true ceteris paribus CATE, while the blue line is
the true mutatis mutandis CATE, both function of $x_{1}^{m}$. The dashed and
erratic lines on the right-hand graph are estimations of the CATE function
using two techniques, described in the next section.

Figure 4: ATE, ceteribs paris $\text{CATE}_{cp}(x_{1})$ and mutatis mutandis
$\text{CATE}_{mm}(x_{1})$ on the left, with an estimate of mutatis mutandis
$\text{CATE}_{mm}(x_{1})$ on the right, from the toy dataset from example 1.
Numerical details are given in Appendix A.1.

###  2.5 Application on birth data

Let us now consider the dataset of all deliveries in the U.S. in
2013.222https://www.cdc.gov/nchs/data_access/Vitalstatsonline.htm Those data
have been intensively used to discuss the “low birth weight paradox”. As
explained in Wilcox (1993, 2001), low birth weight of babies $x$ is strongly
associated with increased neonatal mortality $y$. However, low birth weight
infants born to mothers who smoke $t=1$ usually have lower mortality rates
than low birth weight infants born to nonsmoking mothers $t=0$. Hernández-Díaz
et al. (2006) discussed the birth weight paradox based on causal directed
acyclic graphs as a conceptual framework. Multiple causal models have been
considered. Figure 5 illustrates four situations, using directed acyclic
graphs. In the first case (Figure 5a), birth weight $x$ has a direct effect on
mortality $y$, while smoking $t$ has not. It is also possible to consider a
second case where birth weight $x$, and possibly smoking $t$, have a direct
effect on mortality $y$ (Figure 5b). To increase the plausibility of this
scenario, some known common causes of lower birth weight and mortality,
denoted $z$, can be added (Figure 5c). In this third case, Hernández-Díaz et
al. (2006) claims that the variables $z$ might induce an association between
smoking and mortality, conditional on birth weight $x$. Lastly, a fourth
situation that combines the second and the third can be considered (Figure
5d).

Figure 5: Directed acyclic graphs for the birth weight paradox, when $y$ is
the mortality indicator, $x$ the birth weight and $t$ a smoking indicator. $z$
denotes some possible common causes of infant death, from Hernández-Díaz et
al. (2006).

Here, instead of focusing on newborn mortality (which is an unbalanced
variable, with less than $0.5\%$ mortality rate), we consider
$y=\boldsymbol{1}(\text{non-natural delivery})$. As can be seen in Table 2,
about a third of all deliveries can be considered as “un-natural” (or
“complicated”, involving a least a C-section). Among possible explanatory
variables, we consider the weight of the newborn infant $x_{1}$ and the weight
gain of the mother $x_{2}$. Conditional densities, of
$\boldsymbol{x}=\begin{bmatrix}x_{1}&x_{2}\end{bmatrix}$ given $y$ can be
visualized in Figure 8. To illustrate various techniques based on optimal
transport, we will consider $\text{CATE}(\boldsymbol{x})$,

$\tau(\boldsymbol{x})=\text{CATE}(\boldsymbol{x})=\mathbb{P}\big{[}Y^{\star}_{T\leftarrow
1}=1\big{|}\boldsymbol{X}=\boldsymbol{x}\big{]}-\mathbb{P}\big{[}Y^{\star}_{T\leftarrow
0}=1\big{|}\boldsymbol{X}=\boldsymbol{x}\big{]},$

for several possible “treatment” $t$, that can be visualized in Figure 6, with
either a smoker indicator (for the mother) or a variable indicating whether
the newborn is a boy or not. However, emphasis will be placed on a variable
indicating whether the mother is Black (Afro-American) or not. Conditional
densities of $\boldsymbol{x}$ given $t$ can be visualized in Figure 9. In a
nutshell, we want to address the following questions “what would have been the
probability of a non-natural delivery for a baby of weight $x_{1}$ whose
mother gained weight $x_{2}$ during pregnancy, if the mother had been Afro-
American?” or “if the mother had been smoking?”

| Variable of interest
---|---
|  $y=0$ (natural) |  $y=1$ (non-natural)
$n$ number of observations | 2,221,522 (65.70%) | 1,159,776 (34.30%)
$x_{1}$ weight of newborn | average 3,299 g. | average 3,231 g.
$x_{2}$ weight gain of mother | average 30.02 lbs. | average 31.16 lbs.
| “Treatment”
| $t=0$ | $t=1$
Afro-American variable | non-Black | 2,980,387 | (88.14%) | Black | 400,911 | (11.86%)
smoker variable | non-smoker | 2,959,847 | (91.54%) | smoker | 273,685 | (8.46%)
sex variable | baby boy | 1,730,837 | (51.18%) | baby girl | 1,650,461 | (48.82%)
Table 2: Statistics about the variable of interest $y$, indicating a non-
natural delivery, and two explanatory variables, the weight of the newborn
child ($x_{1}$) and the weight gain of the mother $(x_{2})$, on top; and
statistics about the “treatment” considered at the bottom. $t$$x_{1}$ weight
gain non-natural delivery newborn weight $y$$x_{2}$ smoker black mother
gender of infant Figure 6: Directed acyclic graphs to explain non-natural
deliveries, when $y=\boldsymbol{1}(\text{non-natural delivery})$, $x$ is
either the birth weight of the infant ($x_{1}$), or the weight gain of the
pregnant mother ($x_{2}$), and $t$ is either a smoker indicator (for the
mother), or an indicator that the mother is Black (Afro-American), or that the
baby is a boy.

##  3 Quantile based matching

In this section, we consider the simple case where $x$ is univariate. This
allows us to introduce properties that will be extended more formally in
higher dimension in the next section. Following the example of Section 2.4, we
will propose some techniques to generate a counterfactual version of
$(x,y,t=0)$, or $(x,y_{T\leftarrow 0}^{\star})$, that will be $(x_{T\leftarrow
1},y_{T\leftarrow 1}^{\star})$. In Section 3.1, we will discuss classical
matching techniques, used to match each point in $(y_{i},x_{i},t_{i}=0)$ –in
the control group– with another one in $(y_{j},x_{j},t_{j}=1)$ –in the treated
group– when the two groups have the same size. In Section 3.2, we will suggest
on optimal matching algorithm, to associate individual $i$ (in the control
group) to $j$ (in the treated group), or $j_{i}^{\star}$. Then, in Section
3.3, we will discuss the case where the two groups have different sizes, that
will be called optimal “coupling”. In Section 3.4, we will define an
estimator, the mutatis mutandis CATE,
$\widehat{m}_{1}\big{(}\widehat{\mathcal{T}}(x)\big{)}-\widehat{m}_{0}\big{(}x\big{)}$,
where $\widehat{\mathcal{T}}(x)=\widehat{F}_{1}^{-1}\circ\widehat{F}_{0}(x)$,
with $\widehat{F}_{0}$ and $\widehat{F}_{1}$ denoting the empirical
distribution functions of $x$ conditional on $t=0$ and $t=1$, respectively.
Thus, we will use quantiles to optimal “transport” $x$’s from the control
group to the treated group, formally through the $\mathcal{T}$ mapping.
Finally, in Section 3.5, we will illustrate this on the probability that a
non-natural baby delivery occurs.

###  3.1 Classical matching techniques

To estimate the average treatment effect
$\displaystyle{\tau=\mathbb{E}\big{[}Y^{\star}_{T\leftarrow
1}-Y^{\star}_{T\leftarrow 0}\big{]}}$, a standard technique is to consider
matching techniques to match each point in $(y_{i},x_{i},t_{i}=0)$ or
$(y_{i}^{(0)},x_{i}^{(0)})$ with another one in $(y_{j},x_{j},t_{j}=1)$, or
$(y_{j}^{(1)},x_{j}^{(1)})$. In this coupling approach, we assume that there
are $n$ treated and $n$ non-treated individuals. A treated individual $i$
($t_{i}=1$) is matched to someone in the non-treated group ($t_{j}=0$) that is
close enough for some distance on the set of covariates $\mathcal{X}$,
$j^{\star}_{i}=\displaystyle{\underset{j:t_{j}=0}{\text{argmin}}\\{d(x_{i}^{(0)},x_{j}^{(1)})\\}}$,
so that

$\widehat{\tau}=\frac{1}{n}\sum_{i=1}^{n}\big{(}y^{(1)}_{j^{\star}_{i}}-y_{i}^{(0)}\big{)}=\frac{1}{n}\sum_{i=1}^{n}y_{j^{\star}_{i}}^{(1)}-\frac{1}{n}\sum_{i=1}^{n}y^{(0)}_{i}=\overline{y}^{(1)}-\overline{y}^{(0)},$

since we simply consider a re-ordering of the treated population. But
interestingly, that approach provides a counterfactual version of
$(x_{i},y_{i})$ in the treated population,
$(x_{j^{\star}_{i}},y_{j^{\star}_{i}})$. An algorithm performing such a
matching would be Algorithm 1.

Algorithm 1 Counterfactual matching – “1:1 nearest neighbor matching”
(classical)

$\mathcal{D}\leftarrow\\{(y_{i},\boldsymbol{x}_{i},t_{i})\\}$

function Counterfactual1($\mathcal{D}$)

$\mathcal{D}_{0}\leftarrow$ subset of $\mathcal{D}$ when $t=0$ (size $n$)
shuffled, with indices $i$

$\mathcal{D}_{1}\leftarrow$ subset of $\mathcal{D}$ when $t=1$ (size $n$),
with indices $j$

for $i=1,2,\cdots,n$ do

$j^{\star}_{i}=\displaystyle{\underset{j:t_{j}=1}{\text{argmin}}\\{d(\boldsymbol{x}_{i},\boldsymbol{x}_{j})\\}}$
in $\mathcal{D}_{1}$,

$L_{i}\leftarrow(i,j^{\star}_{i},y_{j^{\star}_{i}}^{(1)}-y^{(0)}_{i})$

remove observation $j^{\star}_{i}$ from $\mathcal{D}_{1}$

end for

return matrix $L$ ($n\times 3$, with $L=(L_{i})$)

end function

This algorithm, introduced by Rubin (1973), is described in Stuart (2010)
under the name “1:1 nearest neighbor matching”, and properties are discussed
in Ho et al. (2007) or Dehejia and Wahba (1999) that focuses on the problem of
not removing selected observations (also called “Greedy Matching”).

Quite naturally, it is possible to define some local version of the previous
quantity using weights or some $k$ nearest neighbors approach, to derive an
estimate of the CATE $\widehat{\tau}(x)$, as in Algorithm 2

$\widehat{\tau}(x)\propto\sum_{i=1}^{n}\omega_{i}(x)\big{(}y_{j^{\star}_{i}}^{(1)}-y^{(0)}_{i}\big{)},$

where weight $\omega_{i}(x)$ are all the higher that $x_{i}$ is close to $x$,
either based on a $k$-nearest neighbors approach
($\omega_{i}(x)=\boldsymbol{1}(i\in V_{{x}}^{k})$, as in Algorithm 2) or based
on a kernel approach ($\omega_{i}(x)=K(|{x}-{x}_{i}|)$ for some kernel $K$).

Algorithm 2 Estimate SCATE (classical, with $k$-NN)

dataset $\mathcal{D}\leftarrow\\{(y_{i},\boldsymbol{x}_{i},t_{i})\\}$,

function scate1($\mathcal{D},k,\boldsymbol{x}$)

$L\leftarrow$ Counterfactual1$(\mathcal{D})$

$V_{\boldsymbol{x}}^{k}\leftarrow$ list of $k$ nearest neighbors of
$\boldsymbol{x}_{i}$’s in $\mathcal{D}_{0}$ close to $\boldsymbol{x}$

for $i\in V_{\boldsymbol{x}}^{k}$ do

$d_{i}\leftarrow L\$d(i)$

end for

return $\displaystyle{\frac{1}{k}\sum_{i\in V_{\boldsymbol{x}}^{k}}d_{i}}$

end function

Unfortunately, that matching mechanism can be very sensitive to the initial
permutation: individuals picked first will have a counterfactual in the
treated group close to them, but it might not be the case for the individuals
picked last. In the next section, we will consider some optimal matching among
individuals in the two populations.

###  3.2 Optimal matching

The matching procedure described previously is characterized by some $n\times
n$ permutation matrix, $P$, with entries in $\\{0,1\\}$, satisfying
$\mathbb{P}\boldsymbol{1}_{n}=\boldsymbol{1}_{n}$ and
$\mathbb{P}^{\star\top}\boldsymbol{1}_{n}=\boldsymbol{1}_{n}$, see Brualdi
(2006). Hence, there is a permutation $\sigma$ of $\\{1,\cdots,n\\}$ such that
$j_{i}^{\star}=\sigma(i)$, and $P$ is the matrix associated with $\sigma$
(that satisfies $\boldsymbol{e}_{i}P=\boldsymbol{e}_{\sigma(i)}$, where
$\boldsymbol{e}_{i}$’s denote the standard basis vector, i.e., a row vector of
length $n$ with $1$ in the $i$-th position and $0$ in every other position).
It is possible to seek an “optimal” permutation: if $C$ is the $n\times n$
matrix that quantifies the distance between individuals in the two groups,
$C_{i,j}=d(x_{i}^{(0)},x_{j}^{(1)})=\delta(x_{i}^{(0)}-x_{j}^{(1)})$, the
optimal matching is solution of

$\min_{P\in\mathcal{P}}\langle
P,C\rangle=\min_{P\in\mathcal{P}}\sum_{i,j}P_{i,j}C_{i,j},$

where $\mathcal{P}$ is the set of permutation matrices, and
$\langle\cdot,\cdot\rangle$ is the Frobenius dot-product. This is also called
Kantorovich’s optimal transport problem, from Kantorovich (1942). If $\delta$
is (strictly) convex –as is the standard Euclidean distance– it can be proven
that this optimal transport problem has a simple solution. Instead of using
$(y_{i}^{(0)},x_{i}^{(0)})$, let $r_{i}^{(0)}$ denote the rank of
$x_{i}^{(0)}$ in $\\{x_{1}^{(0)},\cdots,x_{n}^{(0)}\\}$. Similarly, let
$r_{i}^{(1)}$ denote the rank of $x_{i}^{(1)}$ in the treated dataset
$\\{x_{1}^{(1)},\cdots,x_{n}^{(1)}\\}$. The procedure then becomes simply a
matching based on ranks, in the sense that $j_{i}^{\star}$ satisfies
$r_{j_{i}^{\star}}^{(1)}=r_{i}^{(0)}$, as discussed in Chapter 2 of
Santambrogio (2015). Since ranks are defined on $\\{1,2,\cdots,n\\}$, vectors
$\boldsymbol{r}^{(0)}$ and $\boldsymbol{r}^{(1)}$ correspond to two
permutations of $\\{1,2,\cdots,n\\}$, that we can denote $\sigma_{0}$ and
$\sigma_{1}$, respectively. The optimal coupling is based on permutation
$\sigma=\sigma_{1}\circ\sigma_{0}^{-1}$ in the sense that $x_{i}^{(0)}$ is
associated to $x_{\sigma(i)}^{(1)}$. If the $\boldsymbol{x}^{(0)}$’s and the
$\boldsymbol{x}^{(1)}$’s are sorted, then $P=\mathbb{I}_{n}$, i.e.,
$x_{i}^{(0)}$ is coupled with $x_{i}^{(1)}$. Or, if $\widehat{F}_{0}$ and
$\widehat{F}_{1}$ are the cumulative distribution functions associated with
sample $\boldsymbol{x}^{(0)}$ and $\boldsymbol{x}^{(1)}$, we can see that if
$u\in(0,1)$ is such that $\widehat{F}_{0}^{-1}(u)=x_{i}^{(0)}$, then
$\widehat{F}_{1}^{-1}(u)=x_{i}^{(1)}$, with the exact same $i$.

###  3.3 Optimal coupling

The previous procedure can be extended in the case where the two groups do not
necessarily have the same size. If the two groups $({x}_{i},t_{i}=0)$ and
$({x}_{j},t_{j}=1)$ have different sizes, namely $n_{0}$ and $n_{1}$,
respectively, it is possible to define some matching using weights, and
weighted mean of individuals in the two groups.

In a very general setting, if $\boldsymbol{a}_{0}\in\mathbb{R}_{+}^{n_{0}}$
and $\boldsymbol{a}_{1}\in\mathbb{R}_{+}^{n_{1}}$ satisfy
$\boldsymbol{a}_{0}^{\top}\boldsymbol{}{1}_{n_{0}}=\boldsymbol{a}_{1}^{\top}\boldsymbol{}{1}_{n_{1}}$
(identical sums), define

$U(\boldsymbol{a}_{0},\boldsymbol{a}_{1})=\big{\\{}M\in\mathbb{R}_{+}^{n_{0}\times
n_{1}}:M\boldsymbol{1}_{n_{1}}=\boldsymbol{a}_{0}\text{ and
}{M}^{\top}\boldsymbol{1}_{n_{0}}=\boldsymbol{a}_{1}\big{\\}}.$

This set of matrices is a convex polytope (see Brualdi (2006)). The optimal
coupling is matrix $P^{\star}$ solution of

$\min_{P\in U(\boldsymbol{a}_{0},\boldsymbol{a}_{1})}\left\\{\langle
C,P\rangle\right\\},$

which is solved using linear programming, by casting matrix
$P\in\mathbb{R}_{+}^{n_{0}\times n_{1}}$ as a vector
$\boldsymbol{p}\in\mathbb{R}_{+}^{n_{0}n_{1}}$ such that
$\boldsymbol{p}_{i+n(j-1)}=P_{i,j}$, and similarly for the cost matrix $C$.
The constraint $P\in U(\boldsymbol{a}_{0},\boldsymbol{a}_{1})$ becomes
equivalently

$\begin{pmatrix}\boldsymbol{1}_{n_{0}}^{\top}\otimes\mathbb{I}_{n_{1}}\\\
\mathbb{I}_{n_{0}}\otimes\boldsymbol{1}_{n_{1}}^{\top}\end{pmatrix}\boldsymbol{p}=A\boldsymbol{p}=(\boldsymbol{a}_{0},\boldsymbol{a}_{1})^{\top}=\begin{pmatrix}\boldsymbol{a}_{0}\\\
\boldsymbol{a}_{1}\end{pmatrix},$

where $A$ is some $(n_{0}+n_{1})\times(n_{0}n_{1})$ matrix. The optimal
matching problem is then simply

$\min\left\\{\boldsymbol{c}^{\top}\boldsymbol{p}\right\\}\text{ subject to
}A\boldsymbol{p}=(\boldsymbol{a}_{0},\boldsymbol{a}_{1})^{\top}.$

In our case, let $U_{n_{0},n_{1}}$ denote
$U(\boldsymbol{1}_{0},\frac{n_{0}}{n_{1}}\boldsymbol{1}_{1})$

$P^{*}\in\underset{P\in U_{n_{0},n_{1}}}{\text{argmin}}\langle
P,C\rangle\text{ ou }\underset{P\in
U_{n_{0},n_{1}}}{\text{argmin}}\sum_{i=1}^{n_{0}}\sum_{j=1}^{n_{1}}P_{i,j}C_{i,j}.$
(2)

One can notice that this matrix optimisation problem does not depend on the
dimension of space, so it will easily be extended to the case where $x$ is
multivariate. Nevertheless, in the univariate setting, this approach can be
related to quantile functions.

###  3.4 From optimal matching to CATE

Let $F_{0}$ and $F_{1}$ denote the two conditional distributions of $X$, an
absolutely continuous variable, in the control group ($t=0$) and in the
treatment group ($t=1$), respectively. Then the optimal matching between the
two groups is based on transformation $\mathcal{T}:x_{0}\mapsto
x_{1}=F_{1}^{-1}\circ F_{0}(x_{0})$. From the probability integral transform
property: if $X_{0}\sim F_{0}$, then $F_{0}(X_{0})$ is uniform on the unit
interval $[0,1]$, and then $X_{1}=\mathcal{T}(X_{0})\sim F_{1}$.

######  Lemma 1.

If $X_{0}\sim F_{0}$, then $X_{1}=\mathcal{T}(X_{0})\sim F_{1}$, where
$\mathcal{T}:x_{0}\mapsto x_{1}=F_{1}^{-1}\circ F_{0}(x_{0})$.

######  Definition 3.1.

The mutatis mutandis quantile-based CATE is

$\text{QCATE}(u)=\mathbb{E}\big{[}Y^{*}_{T\leftarrow
1}\big{|}X=F_{1}^{-1}(u)\big{]}-\mathbb{E}\big{[}Y^{*}_{T\leftarrow
0}\big{|}X=F_{0}^{-1}(u)\big{]},$ (3)

where $F_{t}$ is the cumulative distribution function of $X$, conditional on
$T=t$, or

$\text{CATE}(x)=\mathbb{E}\big{[}Y^{*}_{T\leftarrow
1}\big{|}X=\mathcal{T}(x)\big{]}-\mathbb{E}\big{[}Y^{*}_{T\leftarrow
0}\big{|}X=x\big{]},\leavevmode\nobreak\ \mathcal{T}=F_{1}^{-1}\circ F_{0}$
(4)

where $x$ is considered with respect to the control group.

Thus, $\text{CATE}(x)=\text{QCATE}(F_{0}(x))$.

######  Definition 3.2.

Consider two models, $\widehat{m}_{0}(x)$ and $\widehat{m}_{1}(x)$, that
estimate, respectively, $\mathbb{E}[Y|X=x,T=0]$ and $\mathbb{E}[Y|X=x,T=1]$. A
natural estimator of the mutatis mutandis CATE is

$\text{SCATE}(x)=\widehat{m}_{1}\big{(}\widehat{\mathcal{T}}(x)\big{)}-\widehat{m}_{0}\big{(}x\big{)}$

where $\widehat{\mathcal{T}}(x)=\widehat{F}_{1}^{-1}\circ\widehat{F}_{0}(x)$,
$\widehat{F}_{0}$ with $\widehat{F}_{1}$ denoting the empirical distribution
functions of $x$ conditional on $t=0$ and $t=1$, respectively.

Note that a simple parametric transformation can be obtained, based on the
assumption that $X$ conditional on $T$ is Gaussian. More precisely, if
$X_{1}\overset{\mathcal{L}}{=}{X}|t=1\sim\mathcal{N}({\mu}_{1},{\sigma}_{1}^{2})$
and
$X_{0}\overset{\mathcal{L}}{=}{X}|t=0\sim\mathcal{N}({\mu}_{0},{\Sigma}_{0})$,

$\mu_{1}+\sigma_{1}\cdot\frac{X_{0}-\mu_{0}}{\sigma_{0}}\overset{\mathcal{L}}{=}X_{1}$

######  Definition 3.3.

Consider two models, $\widehat{m}_{0}(x)$ and $\widehat{m}_{1}(x)$, that
estimate respectively $\mathbb{E}[Y|X=x,T=0]$ and $\mathbb{E}[Y|X=x,T=0]$. A
Gaussian estimator of the mutatis mutandis CATE is

$\text{SCATE}_{\mathcal{N}}(x)=\widehat{m}_{1}\big{(}\widehat{\mathcal{T}}_{\mathcal{N}}(x)\big{)}-\widehat{m}_{0}\big{(}x\big{)}$

where
$\widehat{\mathcal{T}}_{\mathcal{N}}(x)=\overline{x}_{1}+s_{1}s_{0}^{-1}(x-\overline{x}_{0})$,
$\overline{x}_{0}$ and $\overline{x}_{0}$ being respectively the averages of
$x$ in the two sub-populations, and $s_{0}$ and $s_{1}$ the sample standard
deviations.

An algorithm to compute that estimator is Algorithm 6 (in higher dimension).

###  3.5 Application to non-natural deliveries

In Figure 7, we can visualize $x\mapsto\widehat{\mathcal{T}}(x)$ when $x$ is
either the weight of the newborn infant on the left, or the weight gain of the
mother on the right, when $t$ indicates whether the mother is Black or not.
The $x$-axis is the value of $x$ in the control group ($t=0$) and the $y$-axis
is the value of $x$ in the treated group ($t=1$). On the left, observe that
$x\mapsto\widehat{\mathcal{T}}(x)$ is almost linear, parallel to the first
diagonal, below. This corresponds to the fact that the distribution of $X$
conditional on $T=0$ and $T=1$ are similar, up to a translation (same standard
deviation but different mean if a Gaussian transport
$\widehat{\mathcal{T}}_{\mathcal{N}}$ was considered). On the right,
$x\mapsto\widehat{\mathcal{T}}(x)$ is single-crossing the first diagonal. This
corresponds to the fact that the distribution of $X$ conditional on $T=0$ and
$T=1$ have different variances.

Figure 7: Optimal transport (quantile based) when $X$ is the weight of the
newborn infant on the left, and the weight gain of the mother on the right,
when $T$ indicates whether the mother is Black or not in the middle. The solid
line depicts the transported values while the dashed line is the identity
line. See Figure 26 in Appendix A.2 for similar graphs when $T$ indicates
whether the mother is a smoker or not, or indicates the sex of the newborn.

In Figure 8, we can visualize the conditional distributions of $x$, when $y=0$
and $y=1$ (natural and non-natural deliveries, respectively), when $x$ is the
weight of the baby (on the left) and the weight gain of the mother (on the
right). In Figure 9, we can visualize the conditional distributions of $x$,
when $y=0$ and $y=1$, when $t=0$ and $t=1$, where $t$ denotes whether the
mother is Afro-American or not.

Figure 8: Distribution of the weight of the newborn infant (in grams) on the
left and distribution of the weight gain of the mother on the right,
conditional on the delivery mode, $Y=\boldsymbol{1}(\text{non-natural
delivery})$.

Figure 9: Distribution of the weight of the newborn infant (in grams) on the
left and distribution of the weight gain of the mother on the right, whether
the mother is Black or not. See Figure 17 in Appendix A.2 for similar graphs
when $T$ indicates whether the mother is a smoker or not, or indicates the sex
of the newborn.

In Figure 7, we can visualize the empirical optimal coupling function
$\widehat{\mathcal{T}}:x_{0}\mapsto
x_{1}=\widehat{F}_{1}^{-1}\circ\widehat{F}_{0}(x_{0})$, where
$\widehat{F}_{0}$ and $\widehat{F}_{1}$ denote the empirical distribution
functions of $x$ conditional on $t=0$ and $t=1$, respectively.

In Figures 10 and 11, we can visualize $\widehat{m}_{0}(x)$ and
$\widehat{m}_{1}\big{(}\widehat{\mathcal{T}}(x)\big{)}$ on the left, when $t$
indicates whether the mother is Afro-American or not, when $x$ the weight of
the newborn infant in Figure 10 and when $x$ is the weight gain of the mother
in Figure 11. On the right, we can visualize
$x\mapsto\text{CATE}(x)=\widehat{m}_{1}\big{(}\widehat{\mathcal{T}}(x)\big{)}-\widehat{m}_{0}\big{(}x\big{)}$
as a function of $x$. The light curve in the back is
$\widehat{m}_{1}\big{(}x\big{)}-\widehat{m}_{0}\big{(}x\big{)}$. Numerical
values are given in Table 3 when $x$ is the weight of the newborn, and Table 4
when $x$ is the weight gain of the pregnant mother. For instance, a baby
weighting $2500$g (7.46% quantile in the non-Black population) corresponds to
a baby weighting $2301$g if the mother had been Black. The probability to have
a non-natural delivery has then an additional $5.5\%$ compared with non-Black
mother, using the GAM-SCATE approach. Using a Gaussian transport, the
counterfactual in the Black population is a $2297$g baby, and the probability
to have a non-natural delivery has then an additional $5.60\%$ compared with a
non-Black mother, using the GAM-$\text{SCATE}_{\mathcal{N}}$ approach.
Similarly, a baby weighting $3500$g (64.13% quantile in the non-Black
population) corresponds to a baby weighting $3375$g had the mother been Black
(about $3.6\%$ less). The probability to have a non-natural delivery has then
an additional $4.42\%$ compared with non-Black mother, using the GAM-SCATE
approach. Using a Gaussian transport, estimates are similar.

$t$: mother is Afro-American
---
$x$ (newborn’s weight) | 2000 | 2500 | 3000 | 3500 | 4000 | 4500
$u$ | 2.67% | 7.46% | 25.13% | 64.13% | 91.73% | 98.87%
$\text{CATE}_{0}(x)$ (GAM) | 0.58% | 1.99% | 3.24% | 4.86% | 7.78% | 11.70%
$\widehat{\mathcal{T}}(x)$ | 1595 | 2301 | 2863 | 3375 | 3890 | 4415
$\text{SCATE}(x)$ (GAM) | 7.94% | 5.53% | 4.53% | 4.42% | 5.16% | 7.46%
$\widehat{\mathcal{T}}_{\mathcal{N}}(x)$ | 1758 | 2297 | 2836 | 3376 | 3915 | 4455
$\text{SCATE}_{\mathcal{N}}(x)$ (GAM) | 5.15% | 5.60% | 4.82% | 4.42% | 5.71% | 9.41%
$\text{SCATE}_{\mathcal{N}}(x)$ (kernel) | 6.98% | 6.64% | 4.34% | 4.53% | 5.34% | 7.13%
Table 3: Estimation of the conditional average treatment (CATE), on the
probability to have a non-natural birth ($y$), as a function of the weight of
the baby ($x$, in g.), when the mother is Afro-American. Several weights $x$
are considered, from $2$ to $4.5$kg. $u$ is the probability associated with
$x$, in the baseline population ($t=0$). $\text{CATE}_{0}$ is simply the
difference $\widehat{m}_{1}(x)-\widehat{m}_{0}(x)$, where both
$\widehat{m}_{0}$ and $\widehat{m}_{1}$ are GAMs. $\widehat{\mathcal{T}}(x)$
is the quantile based transport function
($\widehat{\mathcal{T}}(x)=\widehat{F}_{1}^{-1}\circ\widehat{F}_{0}(x)$),
while $\widehat{\mathcal{T}}_{\mathcal{N}}(x)$ is the Gaussian one. Thus,
$\text{SCATE}(x)$ is the mutatis mutandis CATE
$\text{SCATE}(x)=\widehat{m}_{1}\big{(}\widehat{\mathcal{T}}(x)\big{)}-\widehat{m}_{0}\big{(}x\big{)}$,
while
$\text{SCATE}_{\mathcal{N}}(x)=\widehat{m}_{1}\big{(}\widehat{\mathcal{T}_{\mathcal{N}}}(x)\big{)}-\widehat{m}_{0}\big{(}x\big{)}$,
where both $\widehat{m}_{0}$ and $\widehat{m}_{1}$ are GAMs. Finally, the last
estimate is obtained when $\widehat{m}_{0}$ and $\widehat{m}_{1}$ are simple
local averages, using kernels. See Table 5 in Appendix A.2 for similar table
when $T$ indicates whether the mother is a smoker or not, or indicates the sex
of the newborn. $t$: mother is Afro-American
---
$x$ (weight gain of the mother) | 5 | 15 | 25 | 35 | 45 | 55
$u$ | 4.57% | 14.34% | 37.15% | 66.81% | 86.34% | 94.94
$\text{CATE}_{0}(x)$ (GAM) | 3.79% | 4.79% | 5.06% | 4.82% | 4.18% | 3.26%
$\widehat{\mathcal{T}}(x)$ | 1 | 12 | 24 | 35 | 47 | 58
$\text{CATE}(x)$ (GAM) | 5.25% | 5.25% | 5.04% | 4.82% | 4.69% | 4.19%
$\widehat{\mathcal{T}}_{\mathcal{N}}(x)$ | 1 | 12 | 23 | 34 | 46 | 57
$\text{CATE}_{\mathcal{N}}(x)$ (GAM) | 5.22% | 5.21% | 5.03% | 4.74% | 4.33% | 3.78%
$\text{CATE}_{\mathcal{N}}(x)$ (kernel) | 3.78% | 5.49% | 5.31% | 4.49% | 4.12% | 3.61%
Table 4: Estimation of the conditional average treatment (CATE), on the
probability to have a non-natural birth ($y$), as a function of the weight
gain of the mother ($x$, in lbs), when the mother is Afro-American. Several
weight gains $x$ are considered, from $5$ to $55$lbs. $u$ is the probability
associated with $x$, in the baseline population ($t=0$). $\text{CATE}_{0}$ is
simply the difference $\widehat{m}_{1}(x)-\widehat{m}_{0}(x)$, where both
$\widehat{m}_{0}$ and $\widehat{m}_{1}$ are GAMs. $\widehat{\mathcal{T}}(x)$
is the quantile based transport function
($\widehat{\mathcal{T}}(x)=\widehat{F}_{1}^{-1}\circ\widehat{F}_{0}(x)$),
while $\widehat{\mathcal{T}}_{\mathcal{N}}(x)$ is the Gaussian one. Thus,
$\text{SCATE}(x)$ is the mutatis mutandis CATE
$\text{SCATE}(x)=\widehat{m}_{1}\big{(}\widehat{\mathcal{T}}(x)\big{)}-\widehat{m}_{0}\big{(}x\big{)}$,
while
$\text{SCATE}_{\mathcal{N}}(x)=\widehat{m}_{1}\big{(}\widehat{\mathcal{T}_{\mathcal{N}}}(x)\big{)}-\widehat{m}_{0}\big{(}x\big{)}$,
where both $\widehat{m}_{0}$ and $\widehat{m}_{1}$ are GAMs. Finally, the last
estimate is obtained when $\widehat{m}_{0}$ and $\widehat{m}_{1}$ are simple
local averages, using kernels. See Table 6 in Appendix A.2 for similar table
when $T$ indicates whether the mother is a smoker or not, or indicates the sex
of the newborn.

In Figure 12, as previously, $\widehat{m}_{0}(x)$ and
$\widehat{m}_{1}\big{(}\widehat{\mathcal{T}}_{\mathcal{N}}(x)\big{)}$ can be
visualized on the left, when $t$ indicates whether the mother is Afro-American
or not, and when $x$ is the gain weight of the mother. On the right, we can
visualize
$x\mapsto\text{CATE}(x)=\widehat{m}_{1}\big{(}\widehat{\mathcal{T}}(x)\big{)}-\widehat{m}_{0}\big{(}x\big{)}$
as a function of $x$. Numerical values are given in Table 3 when $x$ is the
weight of the newborn, and Table 4 when $x$ is the weight gain of the pregnant
mother.

In Figure 13, some local kernels are used to estimate $\widehat{m}_{0}(x)$ and
$\widehat{m}_{1}\big{(}\widehat{\mathcal{T}}_{\mathcal{N}}(x)\big{)}$ on the
left. Numerical values are given in Table 3 when $x$ is the weight of the
newborn, and Table 4 when $x$ is the weight gain of the pregnant mother.

Figure 10: On the left, evolution of $x\mapsto\mathbb{E}[Y|X_{T\leftarrow
t}=x,T=t]$, estimated using a logistic GAM model, when
$Y=\boldsymbol{1}(\text{non-natural delivery})$, and $X$ is the weight of the
newborn infant, respectively when $T$ indicates whether the mother is Black or
not. On the right, evolution of $x\mapsto\text{SCATE}[Y|X=x]$. See Figure 18
in Appendix A.2 for similar graphs when $T$ indicates whether the mother is a
smoker or not, or indicates the sex of the newborn.

Figure 11: On the left, evolution of $x\mapsto\mathbb{E}[Y|X_{T\leftarrow
t}=x,T=t]$, estimated using a logistic GAM model, when
$Y=\boldsymbol{1}(\text{non-natural delivery})$, and $X$ is the weight gain of
the mother, respectively when $T$ indicates whether the mother is Black or
not. On the right, evolution of $x\mapsto\text{SCATE}[Y|X=x]$. See Figure 19
in Appendix A.2 for similar graphs when $T$ indicates whether the mother is a
smoker or not, or indicates the sex of the newborn.

Figure 12: On the left, evolution of $x\mapsto\mathbb{E}[Y|X_{T\leftarrow
t}=x,T=t]$, estimated using a logistic GAM model, when
$Y=\boldsymbol{1}(\text{non-natural delivery})$, and $X$ is the weight gain of
the mother, respectively when $T$ indicates whether the mother is Black or
not. On the right, evolution of $x\mapsto\text{SCATE}_{\mathcal{N}}[Y|X=x]$.
See Figure 20 in Appendix A.2 for similar graphs when $T$ indicates whether
the mother is a smoker or not.

Figure 13: On the left, evolution of $x\mapsto\mathbb{E}[Y|X_{T\leftarrow
t}=x,T=t]$, estimated using a kernel based local average, when
$Y=\boldsymbol{1}(\text{non-natural delivery})$, and $X$ is the weight of the
newborn infant, respectively when $T$ indicates whether the mother is a smoker
or not (on top), when the mother is Black or not in the middle, and the sex of
the infant below. On the right, evolution of
$x\mapsto\text{SCATE}_{\mathcal{N}}[Y|X=x]$ with an without transport, based
on a Gaussian transport. See Figure 22 in Appendix A.2 for similar graphs when
$T$ indicates whether the mother is a smoker or not, or indicates the sex of
the newborn.

##  4 Optimal transport based matching

In this section, we will extend what was derived in the previous section.
Heuristically, optimal matching of margins components of $\boldsymbol{x}$ will
probably not work, and the mapping should be multivariate. We will therefore
use optimal transport techniques to get a proper counterfactual of
$\boldsymbol{x}$, not in the control group, but in the treated group. In
Section 4.1, we will define properly the optimal transport problem (in any
dimension). Then, in Section 4.2, we will describe how to optimally associate
each observation $\boldsymbol{x}_{i}$ in the control group (when $t=0$) with a
single counterfactual observation $\boldsymbol{x}_{j}$ in the treated group
(when $t=1$), when two groups have the same size. This can be related to the
Gaussian SEM discussed in Section 2.4. In Section 4.3, we will present the
extension when the two groups have different sizes. In Section 4.4, we will
give an explicit formulation for $\mathcal{T}$ when we the distribution of
$\boldsymbol{X}$ conditional on $T$ is assumed to be Gaussian. The application
to non-natural deliveries will finally be discussed in Section 4.5.

###  4.1 Optimal transport

In the mathematical formulation of Monge (1781)’s problem, we want to push a
distribution from $\mathbb{P}_{0}$ to $\mathbb{P}_{1}$ (distributions on
$\mathbb{R}^{k}$, not necessarily in $\mathbb{R}$ as considered in the
previous section). Given
$\mathcal{T}:\mathbb{R}^{k}\rightarrow\mathbb{R}^{k}$, define the “push-
forward” measure,

$\mathbb{P}_{1}(A)=\mathcal{T}_{\\#}\mathbb{P}_{0}(A)=\mathbb{P}_{0}\big{(}\mathcal{T}^{-1}(A)\big{)},\leavevmode\nobreak\
\forall A\subset\mathbb{R}^{k}.$

For instance, when $k=1$, if $F$ is the cumulative distribution of a
univariate random variable $X$ under $\mathbb{P}$ (i.e.,
$F(x)=\mathbb{P}[X\leq x]$) then $\mathbb{Q}=F_{\\#}\mathbb{P}$ is the uniform
distribution on the unit interval $[0,1]$ as well as
$\mathbb{Q}^{\prime}=\overline{F}_{\\#}\mathbb{P}$, where $\overline{F}$ is
the survival function associated with $F$ (i.e.,
$\overline{F}(x)=\mathbb{P}[X>x]$). Similarly, or conversely, if $Q$ is the
quantile function associated with $F$ –$Q(u)=F^{-1}(u)$ for any $u\in(0,1)$–
then if $\mathbb{P}$ is the uniform distribution on the unit interval $[0,1]$,
$\mathbb{Q}=Q_{\\#}\mathbb{P}$ satisfies $\mathbb{Q}[X\leq x]=Q^{-1}(x)=F(x)$,
and similarly for $\overline{Q}$ where $\overline{Q}(u)=F^{-1}(1-u)$.

Observe that if $\mathbb{P}_{0}$ and $\mathbb{P}_{1}$ have densities $f_{0}$
and $f_{1}$, respectively, and if $T$ is continuously differentiable,
$\mathbb{P}_{1}=\mathcal{T}_{\\#}\mathbb{P}_{0}$ is any only if
$f_{0}(\boldsymbol{x})=f_{1}(\mathcal{T}(\boldsymbol{x}))\cdot|\det\nabla\mathcal{T}(\boldsymbol{x})|$,
for all $\boldsymbol{x}$. This non-linear function is a special case of the
so-called Monge-Ampère partial differential equations.

An optimal transport $\mathcal{T}^{\star}$ (in Brenier’s sense, from Brenier
(1991), see Villani (2009) or Galichon (2016)) from $\mathbb{P}_{0}$ towards
$\mathbb{P}_{1}$ will be solution of

$\mathcal{T}^{\star}\in\underset{\mathcal{T}:\mathcal{T}_{\\#}\mathbb{P}_{0}=\mathbb{P}_{1}}{\text{arginf}}\left\\{\int_{\mathbb{R}^{k}}\|\boldsymbol{x}-\mathcal{T}(\boldsymbol{x})\|^{2}d\mathbb{P}_{0}(\boldsymbol{x})\right\\},$

for a quadratic cost, or more generally,

$\mathcal{T}^{\star}\in\underset{\mathcal{T}:\mathcal{T}_{\\#}\mathbb{P}_{0}=\mathbb{P}_{1}}{\text{arginf}}\left\\{\int_{\mathbb{R}^{k}}\gamma(\boldsymbol{x},\mathcal{T}(\boldsymbol{x}))d\mathbb{P}_{0}(\boldsymbol{x})\right\\},$

for some cost function
$\gamma:\mathbb{R}^{k}\times\mathbb{R}^{k}\to\mathbb{R}_{+}$.

If $k=1$, and if the cost function $\gamma$ can be written
$\gamma(x,y)=h(|x-y|)$ for some strictly convex and positive function $h$,
then $T^{\star}$ is an increasing function, and more precisely, if
$F_{0}(x)=\mathbb{P}_{0}[X\leq x]$ and $F_{1}(x)=\mathbb{P}_{1}[X\leq x]$,
with $F_{0}$ absolutely continuous, then
$\mathcal{T}^{\star}(x)=F_{1}^{-1}\circ F_{0}(x)$ satisfies
$\mathcal{T}^{\star}_{\\#}\mathbb{P}_{0}=\mathbb{P}_{1}$ (since
$F_{1}(x)=F_{0}(T^{\star-1}(x))$ and $\mathcal{T}^{\star}$ is optimal. the
quadratic cost function (when $h(x)=x^{2}$) is a particular case. The case
where $h$ is concave was discussed in McCann (1999).

In higher dimension, for a quadratic cost, one can prove (see Villani (2003,
2009) or Galichon (2016)) that $\mathcal{T}^{\star}=\nabla\psi$ where $\psi$
is a convex function.

###  4.2 Empirical version of optimal matching

This transport can be seen as a matching between individuals in the two
groups, both of size $n$, $(\boldsymbol{x}_{i},t_{i}=0)$ and
$(\boldsymbol{x}_{j},t_{j}=1)$, instead of two distributions $\mathbb{P}_{0}$
and $\mathbb{P}_{1}$. If $C$ is a $n\times n$ matrix that quantifies the
distance between individuals in the two groups,
$C_{i,j}=d(\boldsymbol{x}_{i},\boldsymbol{x}_{j})$, the optimal matching is
solution of

$\min_{P\in\mathcal{P}}\langle
P,C\rangle=\min_{P\in\mathcal{P}}\sum_{i=1}^{n}\sum_{j=1}^{n}P_{i,j}C_{i,j},$

where $\mathcal{P}$ is the set of permutation matrices, and
$\langle\cdot,\cdot\rangle$ is the Frobenius dot-product. This is also called
Kantorovich’s optimal transport problem, from Kantorovich (1942).
Interestingly, there are some algorithms that can be used to find that optimal
coupling, or matching, which can, in turn, be used to get a counterfactual for
all individuals in each group.

###  4.3 Empirical version of optimal coupling

If the two groups $(\boldsymbol{x}_{i},t_{i}=0)$ and
$(\boldsymbol{x}_{j},t_{j}=1)$ have different sizes, namely $n_{0}$ and
$n_{1}$, respectively, it is possible to define some matching using weights.
In the coupling case, described previously, $P$ was some $n\times n$
permutation matrix. But here, as in Section 3.3 some $n_{0}\times n_{1}$
matrices will be involved, and similar problems are considered

$\min_{P\in U(\boldsymbol{a}_{0},\boldsymbol{a}_{1})}\langle
P,C\rangle=\min_{P\in
U(\boldsymbol{a}_{0},\boldsymbol{a}_{1})}\sum_{i=1}^{n_{0}}\sum_{j=1}^{n_{1}}P_{i,j}C_{i,j}.$
(5)

And again, assuming Gaussian distributions for $\boldsymbol{X}$ conditional on
$T$ will provide an explicit simple transport formula that can be used to get
an estimation of the mutatis mutandis CATE. This algorithm is given by
Algorithm 4, used to compute the Average Treatment Effect.

Algorithm 3 Counterfactual matching, with optimal matching

$\mathcal{D}\leftarrow\\{(y_{i},\boldsymbol{x}_{i},t_{i})\\}$

function Counterfactual2($\mathcal{D}$)

$\mathcal{D}_{0}\leftarrow$ subset of $\mathcal{D}$ when $t=0$ (size $n_{0}$),
with indices $i$

$\mathcal{D}_{1}\leftarrow$ subset of $\mathcal{D}$ when $t=1$ (size $n_{1}$),
with indices $j$

$C\leftarrow$ matrix $n_{0}\times n_{1}$,
$C_{i,j}=d(\boldsymbol{x}_{i},\boldsymbol{x}_{j})$ between points in
$\mathcal{D}_{0}$ and $\mathcal{D}_{1}$

$P^{*}\leftarrow$ solution of Problem (2)

return matrix $P^{*}$ ($n_{0}\times n_{1}$)

end function

Algorithm 4 Estimate SCATE (optimal matching based)

dataset $\mathcal{D}\leftarrow\\{(y_{i},\boldsymbol{x}_{i},t_{i})\\}$,

function scate2($\mathcal{D},k,\boldsymbol{x}$)

$P\leftarrow$ Counterfactual2$(\mathcal{D})$

$V_{\boldsymbol{x}}^{k}\leftarrow$ list of $k$ nearest neighbors of
$\boldsymbol{x}_{i}$’s in $\mathcal{D}_{0}$ close to $\boldsymbol{x}$

return $\displaystyle{\frac{1}{k}\sum_{i\in
V_{\boldsymbol{x}}^{k}}y_{i}^{1}-{P}_{i}^{\top}\boldsymbol{y}^{0}}$

end function

###  4.4 Counterfactuals for Gaussian covariates

In the general case, there are no simple construction and interpretation of
the optimal mapping $\mathcal{T}^{*}$, as the one we had in the univariate
case, based on quantiles. If it is possible, following Hallin et al. (2021),
to define multivariate quantiles (and therefore to extend concepts defined in
Section 3.4). But here, we will simply consider the multivariate Gaussian
case. Suppose that
$\boldsymbol{X}|t=1\sim\mathcal{N}(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{1})$
and
$\boldsymbol{X}|t=0\sim\mathcal{N}(\boldsymbol{\mu}_{0},\boldsymbol{\Sigma}_{0})$.
There is an explicit expression for the optimal transport, which is simply an
affine map (see Villani (2003) for more details). In the univariate case,
$x_{1}=\mathcal{T}^{*}_{\mathcal{N}}(x_{0})=\mu_{1}+\displaystyle{\frac{\sigma_{1}}{\sigma_{0}}(x_{0}-\mu_{0})}$,
while in the multivariate case, an analogous expression can be derived:

$\boldsymbol{x}_{1}=\mathcal{T}^{*}_{\mathcal{N}}(\boldsymbol{x}_{0})=\boldsymbol{\mu}_{1}+\boldsymbol{A}(\boldsymbol{x}_{0}-\boldsymbol{\mu}_{0}),$

where $\boldsymbol{A}$ is a symmetric positive matrix that satisfies
$\boldsymbol{A}\boldsymbol{\Sigma}_{0}\boldsymbol{A}=\boldsymbol{\Sigma}_{1}$,
which has a unique solution given by
$\boldsymbol{A}=\boldsymbol{\Sigma}_{0}^{-1/2}\big{(}\boldsymbol{\Sigma}_{0}^{1/2}\boldsymbol{\Sigma}_{1}\boldsymbol{\Sigma}_{0}^{1/2}\big{)}^{1/2}\boldsymbol{\Sigma}_{0}^{-1/2}$,
where $\boldsymbol{M}^{1/2}$ is the square root of the square (symmetric)
positive matrix $\boldsymbol{M}$ based on the Schur decomposition
($\boldsymbol{M}^{1/2}$ is a positive symmetric matrix), as described in
Higham (2008).

######  Definition 4.1.

Consider two models, $\widehat{m}_{0}(\boldsymbol{x})$ and
$\widehat{m}_{1}(\boldsymbol{x})$, that estimate, respectively,
$\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x},T=0]$ and
$\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x},T=1]$. A Gaussian estimator of the
mutatis mutandis CATE is

$\text{SCATE}_{\mathcal{N}}(\boldsymbol{x})=\widehat{m}_{1}\big{(}\widehat{\mathcal{T}}_{\mathcal{N}}(\boldsymbol{x})\big{)}-\widehat{m}_{0}\big{(}\boldsymbol{x}\big{)},$

where
$\widehat{\mathcal{T}}_{\mathcal{N}}(\boldsymbol{x})=\overline{\boldsymbol{x}}_{1}+\widehat{\boldsymbol{A}}(\boldsymbol{x}-\overline{\boldsymbol{x}}_{0})$,
with $\overline{\boldsymbol{x}}_{0}$ and $\overline{\boldsymbol{x}}_{1}$
being, respectively, the averages of $x$ in the two sub-populations, and
$\widehat{\boldsymbol{A}}=\widehat{\boldsymbol{\Sigma}}_{0}^{-1/2}\big{(}\widehat{\boldsymbol{\Sigma}}_{0}^{1/2}\widehat{\boldsymbol{\Sigma}}_{1}\widehat{\boldsymbol{\Sigma}}_{0}^{1/2}\big{)}^{1/2}\widehat{\boldsymbol{\Sigma}}_{0}^{-1/2}$
where $\widehat{\boldsymbol{\Sigma}}_{0}$ and
$\widehat{\boldsymbol{\Sigma}}_{1}$ denote the sample variance.

The algorithm to compute that estimate is Algorithm 6.

Algorithm 5 Optimal Gaussian Transport

dataset $\mathcal{D}\leftarrow\\{(y_{i},\boldsymbol{x}_{i},t_{i})\\}$,

function Tgaussian($\mathcal{D}$)

$\mathcal{D}_{0}\leftarrow$ subset of $\mathcal{D}$ when $t=0$

$\mathcal{D}_{1}\leftarrow$ subset of $\mathcal{D}$ when $t=1$

estimate moments of $\boldsymbol{x}$’s $\hat{\boldsymbol{\mu}}_{0}$,
$\hat{\boldsymbol{\mu}}_{1}$, $\hat{\boldsymbol{\Sigma}}_{0}$ and
$\hat{\boldsymbol{\Sigma}}_{1}$, in $\mathcal{D}_{0}$ and $\mathcal{D}_{1}$

$\hat{\boldsymbol{A}}\leftarrow\hat{\boldsymbol{\Sigma}}_{0}^{-1/2}\big{(}\hat{\boldsymbol{\Sigma}}_{0}^{1/2}\hat{\boldsymbol{\Sigma}}_{1}\hat{\boldsymbol{\Sigma}}_{0}^{1/2}\big{)}^{1/2}\hat{\boldsymbol{\Sigma}}_{0}^{-1/2}$

function T($\boldsymbol{x}$)

return
$\hat{\boldsymbol{\mu}}_{1}+\hat{\boldsymbol{A}}(\boldsymbol{x}-\hat{\boldsymbol{\mu}}_{0})$,

end function

return function T

end function

Algorithm 6 Parametric Estimate $\text{SCATE}_{\mathcal{N}}$ (Gaussian
transport)

dataset $\mathcal{D}\leftarrow\\{(y_{i},\boldsymbol{x}_{i},t_{i})\\}$,

$\mathcal{D}_{0}\leftarrow$ subset of $\mathcal{D}$ when $t=0$

$\mathcal{D}_{1}\leftarrow$ subset of $\mathcal{D}$ when $t=1$

$\widehat{m}_{0}\leftarrow$ model to predict $y$ based on $\boldsymbol{x}$,
trained on $\mathcal{D}_{0}$

$\widehat{m}_{1}\leftarrow$ model to predict $y$ based on $\boldsymbol{x}$,
trained on $\mathcal{D}_{1}$

$T\leftarrow$ Tgaussian($\mathcal{D}$)

function Scate3($\widehat{m}_{0},\widehat{m}_{1},T,\boldsymbol{x}$)

return
$\displaystyle{\widehat{m}_{1}(T(\boldsymbol{x}))-\widehat{m}_{0}({\boldsymbol{x}})}$

end function

We should probably stress here that, in the very general case, we should
transport only endogenous variables $\boldsymbol{x}^{m}$ (or mediators) and
not exogenous ones $\boldsymbol{x}^{c}$ (or coliders), as discussed in Section
2.1 (and Figure 2).

###  4.5 Application to non-natural deliveries

Figure 14: Joint distributions of $\boldsymbol{X}$ (weight of the newborn
infant and weight gain of the mother), conditional on the treatment $T$, when
$T$ indicates whether the mother is Black or not on the left. On the right,
vector field associated with optimal Gaussian transport, in dimension two
(weight of the newborn infant and weight gain of the mother). Some numerical
values are given in Table 7. On the right, the origin of the arrow is
$\boldsymbol{x}$ in the control group (non-Black pregnant mother) and the
arrowhead is $\widehat{\mathcal{T}}_{\mathcal{N}}(\boldsymbol{x}$ in the
treated group (Black pregnant mother). See Figure 25 in Appendix A.2 for
similar graphs when $T$ indicates whether the pregnant mother is a smoker or
not.

The left-hand side of Figure 14 displays a scatter plot of
$\boldsymbol{x}=(x_{1},x_{2})$, where $x_{1}$ represents the weight of the
newborn infant while $x_{2}$ shows the weight gain of the mother, conditional
on the treatment $T$, when $T$ indicates whether the mother is Black or not
(see Figure 25 in Appendix A.2 for similar graphs when $T$ indicates whether
the mother is a smoker or not). The ellipses are the iso-density curves under
a Gaussian assumption, such that $95\%$ of the points lie in the ellipse. The
right-hand side of Figure 14, shows $\mathcal{T}_{\mathcal{N}}$ on the same
frame, $\boldsymbol{x}=(x_{1},x_{2})$, with, respectively, the weight of the
newborn infant on the $x$-axis and weight gain of the mother on the $y$-axis.
The origin of an arrow corresponds to $\boldsymbol{x}=(x_{1},x_{2})$, while
its end corresponds to $\widehat{\mathcal{T}}_{\mathcal{N}}(\boldsymbol{x})$.
Note that all the arrows point to the left. Regardless of the weight of the
mother, had the latter been Black, the weight of the newborn would have been
lower. Nevertheless, the length of the arrows varies according to the weight
of the newborn. For infants whose weight is relatively high, for example for
$x_{1}$ close to 4500g, had the mother been Black, the newborn’s weight would
have been almost the same. For newborns whose weight $x_{1}$ is much lower
than 4500g, had the mother been Black, the baby’s weight would have been much
smaller. Some numerical values are given in Table 7 in Appendix A.2. For
instance, if we consider a non-Black mother with a baby weighting 2584g, who
gained 10.8lbs, the counterfactual is a Black mother with a baby weighting
2392g, who gained 7.6lbs.

The top panel of Figures 15 shows the level curves of
$\widehat{m}_{0}:\boldsymbol{x}\mapsto\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x},T=0]$
(left-hand side) and
$\widehat{m}_{1}:\boldsymbol{x}\mapsto\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x},T=1]$
(right-hand side), when the treatment $T$ indicates whether a mother is Black
or not, estimated with logistic GAM models (cubic splines). The middle-level
panel displays curves of the ceteris paribus
$\boldsymbol{x}\mapsto\text{CATE}[\boldsymbol{x}]$ without any transport (on
the left), and $\boldsymbol{x}\mapsto\text{SCATE}[\boldsymbol{x}]$ mutatis
mutandis (on the right). Lastly, the bottom panel shows a positive/negative
distinction for the conditional average treatment effect (positive is red,
negative is blue). Figure 16 provides different results using more knots in
the cubic splines. We can observe that all mothers are more likely to get a
non-natural delivery would they be Black, whatever the weight of the baby (the
ceteris paribus approach would suggest that mothers with small babies, below
2.5kg would be less likely to get a non-natural delivery if they were Black).

Figure 15: On top, contours of
$\boldsymbol{x}\mapsto\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x},T=0]$ and
$\boldsymbol{x}\mapsto\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x},T=1]$ when
$T$ indicates whether a mother is Afro-American or not, estimated with
logistic GAM models (cubic splines). In the middle, contours of the ceteris
paribus $\boldsymbol{x}\mapsto\text{CATE}[\boldsymbol{x}]$ without any
transport on the left, and $\boldsymbol{x}\mapsto\text{SCATE}[\boldsymbol{x}]$
mutatis mutandis on the right. At the bottom, positive/negative distinction
for the conditional average treatment effect. See Figure 28 in Appendix A.2
for similar graphs when $T$ indicates whether the mother is a smoker or not.

Figure 16: On top, contours of
$\boldsymbol{x}\mapsto\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x},T=0]$ and
$\boldsymbol{x}\mapsto\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x},T=1]$ when
$T$ indicates whether a mother is Afro-American or not, estimated with
logistic GAM models (cubic splines, with more knots and degrees of freedom).
In the middle, contours of the ceteris paribus
$\boldsymbol{x}\mapsto\text{CATE}[\boldsymbol{x}]$ without any transport on
the left, and $\boldsymbol{x}\mapsto\text{SCATE}[\boldsymbol{x}]$ mutatis
mutandis on the right. At the bottom, positive/negative distinction for the
conditional average treatment effect. See Figure 28 in Appendix A.2 for
similar graphs when $T$ indicates whether the mother is a smoker or not.

As briefly discussed earlier, optimal matching or coupling in high dimension
can be computationally intensive, since matrices $n_{0}\times n_{1}$ are
involved. For instance, when $t$ is the sex of the newborn, the cost matrix is
a matrix with 3,000 billion entries. Thus, it is quite natural to consider
sub-sampling techniques (since our dataset is quite large). For convenience,
we can use optimal matching on groups of size $n$, and study the robustness of
estimated, as a function of $n$. Some simulations are mentioned in the
Appendix.

##  Appendix A Appendix

###  A.1 Estimation of CATE in a Gaussian framework

With a correlation $r$ (in the simulations, we considered $r=0.4$), consider
the following SEM,

$\displaystyle{\begin{cases}T=\boldsymbol{1}(\epsilon_{t}<0),\\\
{X}^{m}_{1}={\epsilon}_{1m}+(T=1)\cdot\big{(}2+0.2{\epsilon}_{1m})\\\
{X}^{m}_{2}=r{\epsilon}_{1m}+\sqrt{1-r^{2}}\epsilon_{2m}-0.2(T=1)\cdot\big{(}r{\epsilon}_{1m}+\sqrt{1-r^{2}}\epsilon_{2m}\big{)},\\\
{X}^{c}={\epsilon}_{c}\\\
Y=2+T+{X}^{m}_{1}-{X}^{m}_{2}+{X}^{c}+\epsilon_{y}\\\ \end{cases}}$

where $\epsilon$’s are independent $\mathcal{N}(0,1)$ variables. The two
interventions yield

$do(T=0)$ | $do(T=1)$
---|---
$\displaystyle{\begin{cases}T\leftarrow 0\\\ {X}^{m}_{1}=\epsilon_{1m}\\\ {X}^{m}_{2}=r{\epsilon}_{1m}+\sqrt{1-r^{2}}\epsilon_{2m}\\\ {X}^{c}={\epsilon}_{c}\\\ Y_{T\leftarrow 0}=2+{X}^{m}_{1}-{X}^{m}_{2}+{X}^{c}+\epsilon_{y}\\\ \end{cases}}$ | $\displaystyle{\begin{cases}T\leftarrow 1\\\ {X}^{m}_{1}=2+1.2\epsilon^{\prime}_{1m}\\\ {X}^{m}_{2}=0.8r\epsilon^{\prime}_{1m}+\sqrt{1-r^{2}}\epsilon^{\prime}_{2m}\\\ {X}^{c}={\epsilon}^{\prime}_{c}\\\ Y_{T\leftarrow 1}=3+{X}^{m}_{1}-{X}^{m}_{2}+{X}^{c}+\epsilon^{\prime}_{y}\\\ \end{cases}}$

and if we consider $do(T=0)$, when $X_{1}^{m}=x_{1}$, we have

$\displaystyle{\begin{cases}T\leftarrow 0\\\ {X}^{m}_{1}=x_{1}^{m}\\\
{X}^{m}_{2}=rx_{1}^{m}+\sqrt{1-r^{2}}{\epsilon}_{2m}\\\
{X}^{c}={\epsilon}_{c}\\\ Y_{T\leftarrow
0}=2+{x}^{m}_{1}-(rx_{1}^{m}+\sqrt{1-r^{2}}{\epsilon}_{2m})+{\epsilon}_{c}+\epsilon_{y}\\\
\end{cases}}$

while if we consider $do(T=1)$, when $X_{1}^{m}=x_{1}^{\prime}$

$\displaystyle{\begin{cases}T\leftarrow 1\\\ {X}^{m}_{1}=x_{1}^{\prime}\\\
{X}^{m}_{2}=0.8r(x_{1}^{\prime}-2)/1.2+\sqrt{1-r^{2}}\epsilon^{\prime}_{2m}\\\
{X}^{c}={\epsilon}^{\prime}_{c}\\\ Y_{T\leftarrow
1}=3+x_{1}^{m}-\big{(}0.8r(x_{1}^{\prime}-2)/1.2+\sqrt{1-r^{2}}\epsilon^{\prime}_{2m}\big{)}+{\epsilon}^{\prime}_{c}+\epsilon^{\prime}_{y}\\\
\end{cases}}$

so that

$\displaystyle{\begin{cases}Y_{T\leftarrow
0}=2+(1-r){x}^{\prime}_{1}+\sqrt{1-r^{2}}{\epsilon}_{2m}+{\epsilon}_{c}+\epsilon_{y}\\\
Y_{T\leftarrow
1}=3+1.6r/1.2+(1-0.8r/1.2)x_{1}^{\prime}-\sqrt{1-r^{2}}\epsilon^{\prime}_{2m}+{\epsilon}^{\prime}_{c}+\epsilon^{\prime}_{y}\end{cases}}$

Since ${X}^{m}_{1}=\epsilon^{\prime}_{1m}$ when $T\leftarrow 0$, while
${X}^{m}_{1}=2+1.2\epsilon^{\prime}_{1m}$ when $T\leftarrow 1$, it is
legitimate to assume that if $x_{1,T\leftarrow 0}^{m}=x_{1}$, then
$x_{1,T\leftarrow 1}^{m}=2+1.2x_{1,T\leftarrow 0}^{m}$, in a mutatis mutandis
approach,

$\displaystyle{\begin{cases}Y_{T\leftarrow
0}=2+(1-r){x}_{1}+\sqrt{1-r^{2}}{\epsilon}_{2m}+{\epsilon}_{c}+\epsilon_{y}\\\
Y_{T\leftarrow
1}=3+1.6r/1.2+(1-0.8r/1.2)(2+1.2x_{1})-\sqrt{1-r^{2}}\epsilon^{\prime}_{2m}+{\epsilon}^{\prime}_{c}+\epsilon^{\prime}_{y}\end{cases}}$

Thus,

$ATE=\mathbb{E}[Y_{T\leftarrow 1}-Y_{T\leftarrow 0}]=3$

while the mutatis mutandis CATE is

$CATE(x_{1})=\mathbb{E}[Y_{T\leftarrow 1}|x_{1,T\leftarrow
1}^{m}]-\mathbb{E}[Y_{T\leftarrow 0}|x_{1,T\leftarrow 0}^{m}=x_{1}]$

that is

$CATE(x_{1})=\big{[}3+1.6r/1.2+(1-0.8r/1.2)(2+1.2x_{1})\big{]}-\big{[}2+(1-r){x}_{1}\big{]}$

i.e.,

$CATE(x_{1})=3+0.2(1+r){x}_{1}$

that is linear in $x_{1}^{m}$, with slope $0.2(1+r)$ in this mutatis mutandis
case.

###  A.2 Additional applications (smoker and sex of newborn)

In this section, similar graphs to the one presented earlier are produced.

Figure 17: Distribution of the weight of the newborn infant (in grams) on the
left and distribution of the weight gain of the mother on the right, when the
mother is a smoker or not on top, and depending on the sex of the newborn
infant at the bottom. $t$: mother smoker
---
$x$ (newborn’s weight) | 2000 | 2500 | 3000 | 3500 | 4000 | 4500
$u$ | 2.75% | 7.44% | 24.99% | 64.14% | 91.75% | 98.86%
$\text{CATE}_{0}(x)$ (GAM) | -4.41% | -2.55% | -0.69% | 0.79% | 1.97% | 2.50%
$\widehat{\mathcal{T}}(x)$ | 1775 | 2280 | 2802 | 3317 | 3830 | 4337
$\text{SCATE}(x)$ (GAM) | -0.08% | 1.15% | 1.15% | 0.50% | -1.07% | -4.39%
$\widehat{\mathcal{T}}_{\mathcal{N}}(x)$ | 1786 | 2295 | 2805 | 3314 | 3824 | 4333
$\text{SCATE}_{\mathcal{N}}(x)$ (GAM) | -0.28% | 0.88% | 1.12% | 0.50% | -1.15% | -4.53%
$\text{SCATE}_{\mathcal{N}}(x)$ (kernel) | -0.80% | 0.24% | 1.72% | 0.15% | -1.75% | -2.78%
$t$: sex of the newborn
$x$ (newborn’s weight) | 2000 | 2500 | 3000 | 3500 | 4000 | 4500
$u$ | 2.79% | 7.26% | 23.00% | 60.32% | 90.04% | 98.55%
$\text{CATE}_{0}(x)$ (GAM) | -0.24% | -1.66% | -2.24% | -2.00% | -0.77% | 2.38%
$\widehat{\mathcal{T}}(x)$ | 1960 | 2438 | 2892 | 3374 | 3856 | 4338
$\text{SCATE}(x)$ (GAM) | 0.71% | -0.38% | -0.96% | -2.04% | -3.38% | -5.14%
$\widehat{\mathcal{T}}_{\mathcal{N}}(x)$ | 1947 | 2424 | 2901 | 3377 | 3854 | 4331
$\text{SCATE}_{\mathcal{N}}(x)$ (GAM) | 1.02% | -0.08% | -1.07% | -2.05% | -3.41% | -5.43%
$\text{SCATE}_{\mathcal{N}}(x)$ (kernel) | 2.06% | -0.27% | -1.02% | -2.30% | -3.53% | -5.10%
Table 5: Estimation of the conditional average treatment (CATE), on the
probability to have a non-natural birth ($y$), as a function of the weight of
the baby ($x$, in g.), for different “treatments” ($t$): when the mother is a
smoker, and when the baby is a boy. Several weights $x$ are considered, from
$2$ to $4.5$ kg. $u$ is the probability associated with $x$, in the baseline
population ($t=0$). $\text{CATE}_{0}$ is simply the difference
$\widehat{m}_{1}(x)-\widehat{m}_{0}(x)$, where both $\widehat{m}_{0}$ and
$\widehat{m}_{1}$ are GAMs. $\widehat{\mathcal{T}}(x)$ is the quantile based
transport function
($\widehat{\mathcal{T}}(x)=\widehat{F}_{1}^{-1}\circ\widehat{F}_{0}(x)$),
while $\widehat{\mathcal{T}}_{\mathcal{N}}(x)$ is the Gaussian one. Thus,
$\text{SCATE}(x)$ is the mutatis mutandis CATE
$\text{SCATE}(x)=\widehat{m}_{1}\big{(}\widehat{\mathcal{T}}(x)\big{)}-\widehat{m}_{0}\big{(}x\big{)}$,
while
$\text{SCATE}_{\mathcal{N}}(x)=\widehat{m}_{1}\big{(}\widehat{\mathcal{T}_{\mathcal{N}}}(x)\big{)}-\widehat{m}_{0}\big{(}x\big{)}$,
where both $\widehat{m}_{0}$ and $\widehat{m}_{1}$ are GAMs. Finally, the last
estimate is obtained when $\widehat{m}_{0}$ and $\widehat{m}_{1}$ are simple
local averages, using kernels. $t$: mother smoker
---
$x$ (weight gain of the mother) | 5 | 15 | 25 | 35 | 45 | 55
$u$ | 4.61% | 14.50% | 37.49% | 67.12% | 86.54% | 95.05%
$\text{CATE}_{0}(x)$ (GAM) | 0.06% | 0.84% | 0.69% | -0.08% | -1.26% | -2.59%
$\widehat{\mathcal{T}}(x)$ | 1 | 13 | 25 | 37 | 49 | 60
$\text{CATE}(x)$ (GAM) | 1.60% | 1.21% | 0.69% | 0.14% | -0.38% | -1.08%
$\widehat{\mathcal{T}}_{\mathcal{N}}(x)$ | 1 | 13 | 24 | 36 | 48 | 59
$\text{CATE}_{\mathcal{N}}(x)$ (GAM) | 1.63% | 1.29% | 0.71% | 0.01% | -0.71% | -1.33%
$\text{CATE}_{\mathcal{N}}(x)$ (kernel) | 0.31% | 1.03% | 0.98% | 0.29% | -1.30% | -1.08%
$t$: sex of the newborn
$x$ (weight gain of the mother) | 5 | 15 | 25 | 35 | 45 | 55
$u$ | 4.68% | 14.40% | 36.73% | 65.81% | 85.58% | 94.55%
$\text{CATE}_{0}(x)$ (GAM) | -1.79% | -1.60% | -1.60% | -1.73% | -1.90% | -2.02%
$\widehat{\mathcal{T}}(x)$ | 5 | 15 | 25 | 35 | 45 | 55
$\text{CATE}(x)$ (GAM) | -1.79% | -1.60% | -1.60% | -1.73% | -1.90% | -2.02%
$\widehat{\mathcal{T}}_{\mathcal{N}}(x)$ | 4 | 14 | 24 | 34 | 44 | 54
$\text{CATE}_{\mathcal{N}}(x)$ (GAM) | -1.52% | -1.47% | -1.61% | -1.87% | -2.17% | -2.41%
$\text{CATE}_{\mathcal{N}}(x)$ (kernel) | -1.58% | -1.37% | -1.66% | -1.79% | -2.22% | -2.88%
Table 6: Estimation of the conditional average treatment (CATE), on the
probability to have a non-natural birth ($y$), as a function of the weight
gain of the mother ($x$, in lbs), for different “treatments” ($t$): when the
mother is a smoker, and when the baby is a boy. Several weight gains $x$ are
considered, from $5$ to $55$ lbs. $u$ is the probability associated with $x$,
in the baseline population ($t=0$). $\text{CATE}_{0}$ is simply the difference
$\widehat{m}_{1}(x)-\widehat{m}_{0}(x)$, where both $\widehat{m}_{0}$ and
$\widehat{m}_{1}$ are GAMs. $\widehat{\mathcal{T}}(x)$ is the quantile based
transport function
($\widehat{\mathcal{T}}(x)=\widehat{F}_{1}^{-1}\circ\widehat{F}_{0}(x)$),
while $\widehat{\mathcal{T}}_{\mathcal{N}}(x)$ is the Gaussian one. Thus,
$\text{SCATE}(x)$ is the mutatis mutandis CATE
$\text{SCATE}(x)=\widehat{m}_{1}\big{(}\widehat{\mathcal{T}}(x)\big{)}-\widehat{m}_{0}\big{(}x\big{)}$,
while
$\text{SCATE}_{\mathcal{N}}(x)=\widehat{m}_{1}\big{(}\widehat{\mathcal{T}_{\mathcal{N}}}(x)\big{)}-\widehat{m}_{0}\big{(}x\big{)}$,
where both $\widehat{m}_{0}$ and $\widehat{m}_{1}$ are GAMs. Finally, the last
estimate is obtained when $\widehat{m}_{0}$ and $\widehat{m}_{1}$ are simple
local averages, using kernels.

Figure 18: On the left, evolution of $x\mapsto\mathbb{E}[Y|X_{T\leftarrow
t}=x,T=t]$, estimated using a logistic GAM model, when
$Y=\boldsymbol{1}(\text{non-natural delivery})$, and $X$ is the weight of the
newborn infant, respectively when $T$ indicates whether the mother is a smoker
or not (on top), or whether the the newborn infant is a boy at the bottom. On
the right, evolution of $x\mapsto\text{SCATE}[Y|X=x]$.

Figure 19: On the left, evolution of $x\mapsto\mathbb{E}[Y|X_{T\leftarrow
t}=x,T=t]$, estimated using a logistic GAM model, when
$Y=\boldsymbol{1}(\text{non-natural delivery})$, and $X$ is the weight gain of
the mother, respectively when $T$ indicates whether the mother is a smoker or
not on top, or weather the newborn infant is a girl or not, below. On the
right, evolution of $x\mapsto\text{SCATE}[Y|X=x]$.

Figure 20: On the left, evolution of $x\mapsto\mathbb{E}[Y|X_{T\leftarrow
t}=x,T=t]$, estimated using a logistic GAM model, when
$Y=\boldsymbol{1}(\text{non-natural delivery})$, and $X$ is the weight gain of
the mother, respectively when $T$ indicates whether the mother is a smoker or
not. On the right, evolution of $x\mapsto\text{SCATE}_{\mathcal{N}}[Y|X=x]$.

Figure 21: On the left, evolution of $x\mapsto\mathbb{E}[Y|X_{T\leftarrow
t}=x,T=t]$, estimated using $k$-nearest neighbors, when
$Y=\boldsymbol{1}(\text{non-natural delivery})$, and $X$ is the weight gain of
the mother, when $T$ indicates whether the mother is a smoker or not. On the
right, evolution of $x\mapsto\text{SCATE}_{\mathcal{N}}[Y|X=x]$ with and
without transport.

Figure 22: On the left, evolution of $x\mapsto\mathbb{E}[Y|X_{T\leftarrow
t}=x,T=t]$, estimated using a kernel based local average, when
$Y=\boldsymbol{1}(\text{non-natural delivery})$, and $X$ is the weight of the
newborn infant, respectively when $T$ indicates whether the mother is a smoker
or not (on top), whether the mother is Black or not in the middle, and whether
the newborn infant is a boy at the bottom. On the right, evolution of
$x\mapsto\text{SCATE}_{\mathcal{N}}[Y|X=x]$ with an without transport, based
on a Gaussian transport.

Figure 23: On the left, distribution of $\widehat{\mathcal{T}}(x)$, when $x$
is the weight of a newborn infant, when $n$ goes to $4,000$ (on top) to
$20,000$ (in the middle) and $100,000$ (at the bottom), and when $T$ indicates
whether the mother is a smoker or not. On the right, boxplots of the
estimation of $\text{SCATE}(x)$, with two GAM models, when
$x\in\\{2000,2500,\cdots,4000,4500\\}$.

Figure 24: On the left, distribution of $\widehat{\mathcal{T}}(x)$, when $x$
is the weight of a newborn infant, when $n$ goes to $4,000$ (on top) to
$20,000$ (in the middle) and $100,000$ (at the bottom), and when $T$ indicates
whether the mother is Afro-American or not. On the right, boxplots of the
estimation of $\text{SCATE}(x)$, with two GAM models, when
$x\in\\{2000,2500,\cdots,4000,4500\\}$.

$t$: mother is smoker
---
$t=0$ | $t=1$
non-smoker | smoker
$\boldsymbol{x}=(x_{1},x_{2})$ | $\mathcal{T}_{\mathcal{N}}(\boldsymbol{x})$
weight | gain | weight | gain
2584 | 10.8 | 2353.1 | 7.6
2584 | 46.8 | 2414.8 | 49.5
4152 | 10.8 | 3938.1 | 7.9
4152 | 46.8 | 3999.8 | 49.8

$t$: mother is Afro-American
---
$t=0$ | $t=1$
non-Black | Black
$\boldsymbol{x}=(x_{1},x_{2})$ | $\mathcal{T}_{\mathcal{N}}(\boldsymbol{x})$
weight | gain | weight | gain
2584 | 10.8 | 2392.5 | 7.6
2584 | 46.8 | 2382.0 | 47.6
4152 | 10.8 | 4086.1 | 7.5
4152 | 46.8 | 4075.6 | 47.6

$t$: Sex of the newborn
---
$t=0$ | $t=1$
boy | girl
$\boldsymbol{x}=(x_{1},x_{2})$ | $\mathcal{T}_{\mathcal{N}}(\boldsymbol{x})$
weight | gain | weight | gain
2584 | 10.8 | 2513.4 | 10.2
2584 | 46.8 | 2493.0 | 45.9
4152 | 10.8 | 4012.4 | 10.1
4152 | 46.8 | 3992.0 | 45.8

Table 7: Bivariate optimal transport,
$\boldsymbol{x}\mapsto\mathcal{T}_{\mathcal{N}}(\boldsymbol{x})$, for the
three treatments, for four different individuals $\boldsymbol{x}$ in the
control group.

Figure 25: Joint distributions of $\boldsymbol{X}$ (weight of the newborn
infant and weight gain of the mother), conditional on the treatment $T$, when
$T$ indicates whether the mother is a smoker or not on the left. On the right,
vector field associated with optimal Gaussian transport, in dimension two
(weight of the newborn infant and weight gain of the mother), when the
treatment $T$ indicates whether the mother is a smoker or not. Some numerical
values are given in Table 7. On the right, the origin of the arrow is
$\boldsymbol{x}$ in the control group (non-smoker) and the arrowhead is
$\widehat{\mathcal{T}}_{\mathcal{N}}(\boldsymbol{x}$ in the treated group
(smoker)).

Figure 26: Optimal transport (quantile based) when $X$ is the weight of the
newborn infant on the left, and the weight gain of the mother on the right,
when $T$ indicates whether the mother is a smoker on top, and whether the
newborn infant is a boy at the bottom.

Figure 27: On top, contours of
$\boldsymbol{x}\mapsto\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x},T=0]$ and
$\boldsymbol{x}\mapsto\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x},T=1]$ when
$T$ indicates whether a mother is a smoker or not, estimated with logistic GAM
models (cubic splines). In the middle contours of the ceteris paribus
$\boldsymbol{x}\mapsto\text{CATE}[\boldsymbol{x}]$ without any transport on
the left, and $\boldsymbol{x}\mapsto\text{SCATE}[\boldsymbol{x}]$ mutatis
mutandis on the right. At the bottom, positive/negative distinction for the
conditional average treatment effect.

Figure 28: On top, contours of
$\boldsymbol{x}\mapsto\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x},T=0]$ and
$\boldsymbol{x}\mapsto\mathbb{E}[Y|\boldsymbol{X}=\boldsymbol{x},T=1]$ when
$T$ indicates whether a mother is a smoker or not, estimated with logistic GAM
models (cubic splines, with more knots). In the middle contours of the ceteris
paribus $\boldsymbol{x}\mapsto\text{CATE}[\boldsymbol{x}]$ without any
transport on the left, and $\boldsymbol{x}\mapsto\text{SCATE}[\boldsymbol{x}]$
mutatis mutandis on the right. At the bottom, positive/negative distinction
for the conditional average treatment effect.

In Figure 29, $\text{SCATE}(\boldsymbol{x})$ is estimated for various values
of $\boldsymbol{x}$ ($\boldsymbol{x}=(2500,60)$ on top,
$\boldsymbol{x}=4200,60$ in the middle and $\boldsymbol{x}=(2500,20)$ at the
bottom), depending on sample size $n$. The solid line is the average value,
that is quite stable, but, as expected, the confidence interval is quite large
when $n$ is small. In Figure 30, we can visualize the distribution of
$\widehat{m}_{0}(\boldsymbol{x})$ on top,
$\widehat{m}_{1}(\widehat{\mathcal{T}}(\boldsymbol{x}))$ in the middle, and
$\text{SCATE}(\widehat{\mathcal{T}}(\boldsymbol{x}))=\widehat{m}_{1}(\widehat{\mathcal{T}}(\boldsymbol{x}))-\widehat{m}_{0}(\boldsymbol{x})$
at the bottom, estimated on boostrapped samples of size $n=20,000$. On the
left, $\boldsymbol{x}=(2500,20)$ and on the right $\boldsymbol{x}=(2500,60)$.
The two densities are based on the fact that two GAM models are considered,
with more or less knots. Confidence intervals are obtained by bootstrap.

Figure 29: On the left, estimation of
$\boldsymbol{x}\mapsto\text{SCATE}_{\mathcal{N}}[Y\boldsymbol{x}]$, estimated
using a logistic GAM model with Gaussian transport estimated on $n$
observations (with $n$ increasing from $1,000$ to $150,000$), when
$Y=\boldsymbol{1}(\text{non-natural delivery})$, and $\boldsymbol{X}$ is the
weight of the newborn infant and the weight gain of the mother, respectively
when $T$ indicates whether the mother is a smoker or not (on the left),
whether the mother is Black or not (on the right). On top
$\boldsymbol{x}=(2500,60)$, in the middle $\boldsymbol{x}=(4200,60)$ and at
the bottom, $\boldsymbol{x}=(2500,20)$.

Figure 30: Distribution of $\widehat{m}_{0}(\boldsymbol{x})$ on top,
$\widehat{m}_{1}(\widehat{\mathcal{T}}(\boldsymbol{x}))$ in the middle, and
$\text{SCATE}(\widehat{\mathcal{T}}(\boldsymbol{x}))=\widehat{m}_{1}(\widehat{\mathcal{T}}(\boldsymbol{x}))-\widehat{m}_{0}(\boldsymbol{x})$
at the bottom, estimated on boostrapped samples of size $n=20,000$. On the
left, $\boldsymbol{x}=(2500,20)$ and on the right $\boldsymbol{x}=(2500,60)$.
The two densities are based on the fact that two GAM models are considered.

## Acknowledgments

Arthur Charpentier acknowledges the financial support of the AXA Research Fund
through the joint research initiative use and value of unusual data in
actuarial science, as well as NSERC grant 2019-07077.

Emmanuel Flachaire and Ewen Gallic acknowledge the financial support of the
French National Research Agency Grant ANR-17-EURE-0020, the Excellence
Initiative of Aix Marseille University – A*MIDEX

## References

  * Abrevaya et al. (2015) Abrevaya, J., Hsu, Y.-C. and Lieli, R. P. (2015). Estimating conditional average treatment effects. Journal of Business & Economic Statistics 33: 485–505, doi:10.1080/07350015.2014.975555.
  * Athey et al. (2019) Athey, S., Tibshirani, J. and Wager, S. (2019). Generalized random forests. The Annals of Statistics 47: 1148–1178, doi:10.1214/18-aos1709.
  * Athey and Wager (2019) Athey, S. and Wager, S. (2019). Estimating treatment effects with causal forests: An application. Observational Studies 5: 37–51, doi:10.1353/obs.2019.0001.
  * Brenier (1991) Brenier, Y. (1991). Polar factorization and monotone rearrangement of vector-valued functions. Communications on pure and applied mathematics 44: 375–417, doi:10.1002/cpa.3160440402.
  * Brualdi (2006) Brualdi, R. A. (2006). Combinatorial matrix classes, 13. Cambridge University Press, doi:10.1017/CBO9780511721182.
  * Charpentier (2023) Charpentier, A. (2023). Quantifying fairness and discrimination in predictive models. In Kreinovich, V., SriboonchiNa, S. and Yamaka, W. (eds), Machine Learning for Econometrics and Related Topics. Springer Verlag.
  * Chisholm (1946) Chisholm, R. M. (1946). The contrary-to-fact conditional. Mind 55: 289–307, doi:10.1093/mind/lv.219.289.
  * Cunningham (2021) Cunningham, S. (2021). Causal inference. Yale University Press, doi:10.12987/9780300255881.
  * Davis and Heller (2017) Davis, J. and Heller, S. B. (2017). Using causal forests to predict treatment heterogeneity: An application to summer jobs. American Economic Review 107: 546–50, doi:10.1257/aer.p20171000.
  * Dehejia and Wahba (1999) Dehejia, R. H. and Wahba, S. (1999). Causal effects in nonexperimental studies: Reevaluating the evaluation of training programs. Journal of the American Statistical Association 94: 1053–1062, doi:10.1080/01621459.1999.10473858.
  * Fan et al. (2022) Fan, Q., Hsu, Y.-C., Lieli, R. P. and Zhang, Y. (2022). Estimation of conditional average treatment effects with high-dimensional data. Journal of Business & Economic Statistics 40: 313–327, doi:10.1080/07350015.2020.1811102.
  * Galichon (2016) Galichon, A. (2016). Optimal transport methods in economics. Princeton University Press, doi:10.23943/princeton/9780691172767.001.0001.
  * Goodman (1947) Goodman, N. (1947). The problem of counterfactual conditionals. The Journal of Philosophy 44: 113–128, doi:10.2307/2019988.
  * Hahn (1998) Hahn, J. (1998). On the role of the propensity score in efficient semiparametric estimation of average treatment effects. Econometrica : 315–331doi:10.2307/2998560.
  * Hallin et al. (2021) Hallin, M., Del Barrio, E., Cuesta-Albertos, J. and Matrán, C. (2021). Distribution and quantile functions, ranks and signs in dimension d: A measure transportation approach. The Annals of Statistics 49: 1139–1165, doi:10.1214/20-aos1996.
  * Heckman et al. (1998) Heckman, J. J., Ichimura, H. and Todd, P. (1998). Matching as an econometric evaluation estimator. The Review of Economic Studies 65: 261–294, doi:10.1111/1467-937x.00044.
  * Hernán and Robins (2010) Hernán, M. A. and Robins, J. M. (2010). Causal inference. CRC Press.
  * Hernández-Díaz et al. (2006) Hernández-Díaz, S., Schisterman, E. F. and Hernán, M. A. (2006). The birth weight “paradox” uncovered? American journal of epidemiology 164: 1115–1120, doi:10.1093/aje/kwj275.
  * Higham (2008) Higham, N. J. (2008). Functions of matrices: theory and computation. Society for Industrial and Applied Mathematics, doi:10.1137/1.9780898717778.
  * Hitsch and Misra (2018) Hitsch, G. J. and Misra, S. (2018). Heterogeneous treatment effects and optimal targeting policy evaluation. Available at SSRN 3111957 doi:10.2139/ssrn.3111957.
  * Ho et al. (2007) Ho, D. E., Imai, K., King, G. and Stuart, E. A. (2007). Matching as nonparametric preprocessing for reducing model dependence in parametric causal inference. Political analysis 15: 199–236, doi:10.1093/pan/mpl013.
  * Hsu et al. (2022) Hsu, Y.-C., Lai, T.-C. and Lieli, R. P. (2022). Counterfactual treatment effects: Estimation and inference. Journal of Business & Economic Statistics 40: 240–255, doi:10.1080/07350015.2020.1800479.
  * Imai (2018) Imai, K. (2018). Quantitative social science: an introduction. Princeton University Press.
  * Imbens and Rubin (2015) Imbens, G. W. and Rubin, D. B. (2015). Causal inference in statistics, social, and biomedical sciences. Cambridge University Press, doi:10.1017/cbo9781139025751.
  * Kantorovich (1942) Kantorovich, L. V. (1942). On the translocation of masses. In Doklady Akademii Nauk USSR, 37, 199–201.
  * Künzel et al. (2019) Künzel, S. R., Sekhon, J. S., Bickel, P. J. and Yu, B. (2019). Metalearners for estimating heterogeneous treatment effects using machine learning. Proceedings of the national academy of sciences 116: 4156–4165, doi:10.1073/pnas.1804597116.
  * Li et al. (2018) Li, F., Morgan, K. L. and Zaslavsky, A. M. (2018). Balancing covariates via propensity score weighting. Journal of the American Statistical Association 113: 390–400, doi:10.1080/01621459.2016.1260466.
  * Mach (1893) Mach, E. (1893). The science of mechanics: A critical and historical exposition of its principles. Open court publishing Company.
  * McCann (1999) McCann, R. J. (1999). Exact solutions to the transportation problem on the line. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 455: 1341–1380, doi:10.1098/rspa.1999.0364.
  * Monge (1781) Monge, G. (1781). Mémoire sur la théorie des déblais et des remblais. Histoire de l’Académie Royale des Sciences de Paris .
  * Morgan and Winship (2014) Morgan, S. L. and Winship, C. (2014). Counterfactuals and causal inference. Cambridge University Press, doi:10.1017/cbo9781107587991.
  * Pearl and Mackenzie (2018) Pearl, J. and Mackenzie, D. (2018). The book of why: the new science of cause and effect. Basic books.
  * Powers et al. (2018) Powers, S., Qian, J., Jung, K., Schuler, A., Shah, N. H., Hastie, T. and Tibshirani, R. (2018). Some methods for heterogeneous treatment effect estimation in high dimensions. Statistics in medicine 37: 1767–1787, doi:10.1002/sim.7623.
  * Rosenbaum and Rubin (1983) Rosenbaum, P. R. and Rubin, D. B. (1983). The central role of the propensity score in observational studies for causal effects. Biometrika 70: 41–55, doi:10.1093/biomet/70.1.41.
  * Rubin (1973) Rubin, D. B. (1973). Matching to remove bias in observational studies. Biometrics : 159–183doi:10.2307/2529684.
  * Rubin (1974) Rubin, D. B. (1974). Estimating causal effects of treatments in randomized and nonrandomized studies. Journal of educational Psychology 66: 688, doi:10.1037/h0037350.
  * Santambrogio (2015) Santambrogio, F. (2015). Optimal transport for applied mathematicians. doi:10.1007/978-3-319-20828-2.
  * Stuart (2010) Stuart, E. A. (2010). Matching methods for causal inference: A review and a look forward. Statistical Science 25: 1, doi:10.1214/09-sts313.
  * Villani (2003) Villani, C. (2003). Topics in optimal transportation, 58. American Mathematical Society, doi:10.1090/gsm/058.
  * Villani (2009) Villani, C. (2009). Optimal transport: old and new, 338. Springer Berlin, Heidelberg, doi:10.1007/978-3-540-71050-9.
  * Vowels et al. (2022) Vowels, M. J., Camgoz, N. C. and Bowden, R. (2022). D’ya like DAGs? A survey on structure learning and causal discovery. ACM Computing Surveys (CSUR) 55, doi:10.1145/3527154.
  * Wager and Athey (2018) Wager, S. and Athey, S. (2018). Estimation and inference of heterogeneous treatment effects using random forests. Journal of the American Statistical Association 113: 1228–1242, doi:10.1080/01621459.2017.1319839.
  * Wilcox (1993) Wilcox, A. J. (1993). Birth weight and perinatal mortality: the effect of maternal smoking. American journal of epidemiology 137: 1098–1104, doi:10.1093/oxfordjournals.aje.a116613.
  * Wilcox (2001) Wilcox, A. J. (2001). On the importance—and the unimportance—of birthweight. International journal of epidemiology 30: 1233–1241, doi:10.1093/ije/30.6.1233.