Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
99,220 characters · 22 sections · 51 citation commands
Optimal Transport for Counterfactual Estimation: A Method for Causal Inference
\
{\bf Keywords} Causality; Conditional Average Treatment Effects (CATE); Counterfactual; Mutatis Mutandis; Optimal Transport; Quantiles
In pearl2018book, a “ladder of causation” is introduced, to describe the three levels of causal reasoning. The first level, named “{\em association}”, discusses associations (not to use the word “{\em correlation}”) between variables. Questions such as “{\em 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 “{\em intervention}”. Reasoning on this level answers questions of the form “{\em if I make the intervention $T$, how will this affect the level of the outcome $Y$?}” For example, the question “{\em 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 “{\em 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 “{\em would that person have had surgery if she had been Afro-American?}”
Consider, as in rubin1974estimating or hernan2010causal, 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 imbens2015causal or imai2018quantitative, or $y^1$ and $y^0$ in morgan2014counterfactuals or cunningham2021causal, even $y_{t=1}$ and $y_{t=0}$ in pearl2018book), 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 (ref).
We will use the term “treatment” (and letter $t$) even if interventions are not possible, so it is no {\em 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 charpentier2023fairness, such as the race of an individual, e.g., “{\em 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\leftarrow0}\neq y^\star_{T\leftarrow1}$. More specifically, the causal effect for individual $i$ is ${\tau_i=y^\star_{i,T\leftarrow1} - y^\star_{i,T\leftarrow0}}$. The {average treatment effect} (ATE) can the be defined as follows: $$ \tau = \text{ATE} = \mathbb{E}\big[ Y^\star_{i,T\leftarrow1} - Y^\star_{i,T\leftarrow0}\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\leftarrow1} - y^\star_{i,T\leftarrow0}. $$ 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\leftarrow1} - Y^\star_{T\leftarrow0}\big\vert \boldsymbol{X}=\boldsymbol{x}\big] $$ that can be written $$ \tau(\boldsymbol{x}) = \text{CATE}(\boldsymbol{x}) = \mathbb{E}\big[ Y^\star_{T\leftarrow1} \big\vert \boldsymbol{X}=\boldsymbol{x}\big]- \mathbb{E}\big[ Y^\star_{T\leftarrow0} \big\vert \boldsymbol{X}=\boldsymbol{x}\big], $$ as introduced in hahn1998role and heckman1998matching. More recently, hitsch2018heterogeneous used that measure to quantify heterogeneous treatment effects to evaluate optimal targeting policies, as well as powers2018some and fan2022estimation. wager2018estimation, athey2019estimating and athey2019generalized suggested to use random forests to estimate this quantity, inspired by davis2017using. See also kunzel2019metalearners or hsu2022counterfactual 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})$. rosenbaum1983central suggested a strong “ignorable treatment assignment” assumption defined as a conditional independence between $(Y^\star_{T\leftarrow0},Y^\star_{T\leftarrow1})$ and $T$, conditional on the covariates $\boldsymbol{X}$.
In Section (ref), and more specifically in Section (ref), 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 (ref), we suggest an extension from the standard {\em 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 {\em 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 (ref), 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 “{\em 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 (ref), while in Section (ref), 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 (ref)).
In Section (ref), we will focus on the case where only one covariate $x$ is considered. We will start with classical matching techniques in Section (ref), 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 (ref), 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 (ref), we will discuss the case where the two groups have different sizes, that will be called optimal “coupling”. In Section (ref), we will define an estimator, the {\em 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 (ref), we will illustrate this on probability to have a non-natural baby delivery, on our dataset.
In Section (ref), 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 (ref), we will define the optimal transport problem for any number of dimensions and then, in Section (ref), 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 (ref). In Section (ref), 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 (ref), 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 (ref), those techniques will be further discussed in the context of the application to non-natural birth\footnote{See \href{https://github.com/3wen/counterfactual-estimation-optimal-transport}{ https://github.com/3wen/counterfactual-estimation-optimal-transport} for more details}..
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}$.
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 mach1893science), to answer a question such as “{\em what if $t$ had taken another value?}”, without being able to make an experiment. chisholm1946contrary introduced the idea of “{\em contrary-to-fact conditional}”, coined as “{\em counterfactual}” in goodman1947problem. 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 (ref)).
The Markov assumption, on causal networks, states that each variable is conditionally independent of its non-descendants, given its parents. In Figure (ref), 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., vowels2021d for a survey on techniques used to discover causal structures).
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 (ref). 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 pearl2018book):
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 {\em 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 {\em 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:
More specifically, when we ask the question “{\em 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.
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_iy_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 abrevaya2015estimating, 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_iy_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_iy_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 li2018balancing.
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):
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 (ref), 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).
The two interventions yield
more precisely, in that model with three covariates, $\boldsymbol{X}^m=(X_1^m,X_2^m)$, and since $$ \boldsymbol{\Sigma}_t=
and \boldsymbol{\Sigma}_t^{1/2}=
$$ we can write
and therefore $${
}.$$ Hence, $$ ATE = \mathbb{E}[Y_{T\leftarrow 1}-Y_{T\leftarrow 0}]=\gamma. $$ For conditional average treatment effects, $${
}.$$ {\em Ceteris paribus}, we suppose that $x_1'=x_1$, then $$ 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] =ATE +\delta x_1 +\kappa, $$ where $$
. $$ {\em 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'=x_{1:T\leftarrow 1}=\mu_{11} + \sigma_{11}(\sigma_{01}^{-1}[x_1-\mu_{01}])$. Therefore, {\em mutatis mutandis},$$ CATE_{mm}(x_1) =ATE+ \delta' x_1 +\kappa', $$ where $$
, $$ so that we can also write $$ CATE_{mm}(x_1) = CATE_{cp}(x_1)+\big(dx_1+k\big). $$
In Figure (ref), the horizontal orange line is the true average treatment effect (ATE). The green line is the true {\em ceteris paribus} CATE, while the blue line is the true {\em 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.
Let us now consider the dataset of all deliveries in the U.S. in 2013.\footnote{\href{https://www.cdc.gov/nchs/data_access/Vitalstatsonline.htm}{https://www.cdc.gov/nchs/data_access/Vitalstatsonline.htm}} Those data have been intensively used to discuss the “low birth weight paradox”. As explained in wilcox1993birth,wilcox2001importance, 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$. hernandez2006birth discussed the birth weight paradox based on causal directed acyclic graphs as a conceptual framework. Multiple causal models have been considered. Figure (ref) illustrates four situations, using directed acyclic graphs. In the first case (Figure (ref)a), 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 (ref)b). To increase the plausibility of this scenario, some known common causes of lower birth weight and mortality, denoted $z$, can be added (Figure (ref)c). In this third case, hernandez2006birth 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 (ref)d).
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 (ref), 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}=
$ given $y$ can be visualized in Figure~\ref{fig:densite-1x3-conditional-k}. To illustrate various techniques based on optimal transport, we will consider $CATE(\boldsymbol{x})$, $$\tau(\boldsymbol{x}) = \text{CATE}(\boldsymbol{x}) = \mathbb{P}\big[ Y^\star_{T\leftarrow1} =1\big\vert \boldsymbol{X}=\boldsymbol{x}\big]- \mathbb{P}\big[ Y^\star_{T\leftarrow0}=1 \big\vert \boldsymbol{X}=\boldsymbol{x}\big], $$ for several possible ``treatment'' $t$, that can be visualized in Figure~\ref{Fig:NN:1ex}, 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~\ref{fig:densite-2x3-conditional-t}. In a nutshell, we want to address the following questions ``{\em 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 “{\em if the mother had been smoking?}”
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 (ref), 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 (ref), 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 (ref), 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 (ref), we will discuss the case where the two groups have different sizes, that will be called optimal “coupling”. In Section (ref), we will define an estimator, the {\em 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 (ref), we will illustrate this on the probability that a non-natural baby delivery occurs.
To estimate the average treatment effect $\displaystyle{\tau = \mathbb{E}\big[ Y^\star_{T\leftarrow1} - Y^\star_{T\leftarrow0}\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 (ref).
This algorithm, introduced by rubin1973matching, is described in stuart2010matching under the name “1:1 nearest neighbor matching”, and properties are discussed in ho2007matching or dehejia1999causal 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 (ref) $$ \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 (ref)) or based on a kernel approach ($\omega_i(x) =K(|{x}-{x}_i|)$ for some kernel $K$).
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.
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 brualdi2006combinatorial. 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 kantorovich1942translocation. 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 santambrogio2015optimal. 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$.
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\lbrace 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\rbrace. $$ This set of matrices is a convex polytope (see brualdi2006combinatorial). The optimal coupling is matrix $P^\star$ solution of $$ \min_{P\in U(\boldsymbol{a}_0,\boldsymbol{a}_1)} \left\lbrace \langle C,P\rangle \right\rbrace, $$ 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_0n_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 $$
\boldsymbol{p} = A\boldsymbol{p} =(\boldsymbol{a}_0,\boldsymbol{a}_1)^\top=
, $$ where $A$ is some $(n_0+n_1)\times(n_0n_1)$ matrix. The optimal matching problem is then simply $$ \min\left\lbrace \boldsymbol{c}^\top\boldsymbol{p} \right\rbrace 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)$
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.
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$.
Thus, $\text{CATE}(x) =\text{QCATE}(F_0(x))$.
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 $$
An algorithm to compute that estimator is Algorithm (ref) (in higher dimension).
In Figure (ref), 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.
In Figure (ref), 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 (ref), 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.
In Figure (ref), 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 (ref) and (ref), 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 (ref) and when $x$ is the weight gain of the mother in Figure (ref). 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 (ref) when $x$ is the weight of the newborn, and Table (ref) 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.
In Figure (ref), 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 (ref) when $x$ is the weight of the newborn, and Table (ref) when $x$ is the weight gain of the pregnant mother.
In Figure (ref), 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 (ref) when $x$ is the weight of the newborn, and Table (ref) when $x$ is the weight gain of the pregnant mother.
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 (ref), we will define properly the optimal transport problem (in any dimension). Then, in Section (ref), 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 (ref). In Section (ref), we will present the extension when the two groups have different sizes. In Section (ref), 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 (ref).
In the mathematical formulation of monge1781memoire'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 “{\em push-forward}” measure, $$ \mathbb{P}_1(A)= \mathcal{T}_{\#}\mathbb{P}_0(A)= \mathbb{P}_0\big(\mathcal{T}^{-1}(A)\big),~\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}'=\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 brenier1991polar, see villani2009optimal or galichon2016optimal) 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\lbrace\int_{\mathbb{R}^k} \|\boldsymbol{x}-\mathcal{T}(\boldsymbol{x})\|^2d\mathbb{P}_0(\boldsymbol{x})\right\rbrace, $$ 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\lbrace\int_{\mathbb{R}^k} \gamma(\boldsymbol{x},\mathcal{T}(\boldsymbol{x}))d\mathbb{P}_0(\boldsymbol{x})\right\rbrace, $$ 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 mccann1999exact.
In higher dimension, for a quadratic cost, one can prove (see villani2003optimal,villani2009optimal or galichon2016optimal) that $\mathcal{T}^\star=\nabla \psi$ where $\psi$ is a convex function.
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 kantorovich1942translocation. 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.
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 (ref) some $n_0\times n_1$ matrices will be involved, and similar problems are considered
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 {\em mutatis mutandis} CATE. This algorithm is given by Algorithm (ref), used to compute the Average Treatment Effect.
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 hallin2021distribution, to define multivariate quantiles (and therefore to extend concepts defined in Section (ref)). 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 villani2003optimal 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 higham2008functions.
The algorithm to compute that estimate is Algorithm (ref).
We should probably stress here that, in the very general case, we should transport {\em only} endogenous variables $\boldsymbol{x}^m$ (or mediators) and not exogenous ones $\boldsymbol{x}^c$ (or coliders), as discussed in Section (ref) (and Figure (ref)).
The left-hand side of Figure (ref) 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 (ref) in Appendix (ref) 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 (ref), 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 (ref) in Appendix (ref). 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 (ref) 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 {\em ceteris paribus} $\boldsymbol{x}\mapsto\text{CATE}[\boldsymbol{x}]$ without any transport (on the left), and $\boldsymbol{x}\mapsto\text{SCATE}[\boldsymbol{x}]$ {\em 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 (ref) 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 {\em 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).
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.