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.
51,017 characters · 25 sections · 44 citation commands
Optimal Transport Weights for Causal Inference
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 \fi
\if10 {
} \fi
{\it Keywords:} balancing weights, semiparametric efficiency, IPW, Sinkhorn divergence
\doublespacing
Inverse propensity score weighting (IPW) uses the predicted probabilities of treatment to estimate causal effects. Under the correct model, these weights will lead to distributional balance between treatment groups and thereby to unbiased effect estimates. Unfortunately, adequate distributional balance is a challenge to determine, so researchers often assess performance by measuring the balance of observed covariate functions li_balancing_2018.
Rather than using weights that indirectly balance selected covariate functions, other authors construct weights that achieve such balance by design Hainmueller2012, Imai2014, Zubizarreta2015; however, the particular aspect of the covariates that will achieve distributional balance, especially in finite samples, is usually unknown. Weights built using reproducing kernel Hilbert spaces (RKHS) may get around this difficulty because they theoretically balance all covariate functions Li2021. Nonetheless, RKHS weights necessitate the tuning of unknown hyperparameters without an obvious metric to assess performance. Some authors tune the RKHS weights such that they balance predictions from an outcome model Kallus2016,Kallus2018, meaning the performance of the weights in turn depends on the performance of these outcome models. However, all of these balancing methods fail to capture something fundamental about the covariates: the joint distribution.
Ultimately, if the distributions are the same between treatment groups, then all functions of the covariates will be balanced, leading to unbiased treatment effect estimates. But researchers may have certain covariate functions that they believe are important a priori. As such, researchers may not want to trade off balance on these quantities for better overall distributional balance. Fortunately, we do not have to choose between balancing distributions or covariate functions---we can do both.
Causal Optimal Transport (COT) is our proposed method that balances the joint distribution and any other selected covariate functions of interest in one framework. As we document, our method is semiparametrically efficient and performs well in simulation studies compared to competing methods---especially when the propensity score and outcome models are misspecified. We also demonstrate that several methods in the literature are actually special cases of COT, meaning our framework can be seen as an interpolation between several seemingly unconnected methods. Finally, we show how optimal transport methods can nonparametrically impute the missing potential outcomes by estimating the conditional mean outcome function at a $n^{-2/(2d' + 9)}$ rate, where $d' = \lceil 5d/4\rceil$ and $d$ is the dimension of the covariate space. \if00{Code to implement the methods discussed in this work is available in the new R package causalOT, found at \url{http://www.github.com/ericdunipace/causalOT}.}\fi
The paper proceeds as follows: in Section (ref), we describe the setting and assumptions necessary for identification. Then we introduce optimal transport and COT in Section (ref). In Section (ref), we demonstrate how COT unifies several existing methods, and in Section (ref), we provide simulation results demonstrating the utility of COT. Section (ref) presents a case study utilizing our method in a study of post-partum hemorrhage. Finally, we offer our concluding remarks in Section (ref).
We adopt the potential outcomes framework of Neyman and Rubin Splawa-Neyman1923, Rubin1974. Assume that we have an independent, identically distributed (iid) sample of $n \in \Nat$ units from some population. Let $Z$ be a binary variable that denotes receiving either a treatment ($Z = 1$) or control ($Z = 0$) condition. $Y(0)$ and $Y(1)$ are the potential outcomes, and $Y = Z \cdot Y(1) + (1-Z) \cdot Y(0)$ is the observed outcome defined on a space $\mcY \subseteq \R$. The confounders are $X \in \mcX \subseteq \R^d$. We will assume we have $n_0 \in \Nat$ control units and $n_1 \in \Nat $ treated units giving $n = n_0 + n_1$ total observations from this sample. Denote $\bolda_z = \sum_{i: Z_i = z}^n \delta_{x_i} a_i$ as the empirical distribution for treatment $Z = z$ and $\bolda = \sum_i^n \delta_{x_i} a_i$ as the empirical distribution for the full sample. Let $\alpha_z$ and $\alpha$ be the corresponding population distributions.
Finally, we assume the space $\mcX$ has a distance metric between observations, $d_{\mcX}(x_i, x_j) \in \R_+$. We will define a generic cost function as $c(x_i, x_j) = d_{\mcX}(x_i, x_j) ^p$ with $p \geq 1$. As an example, if $d_\mcX$ is the Euclidean distance and $p = 2$, then $c$ is the squared-Euclidean distance. From this function, we then construct a pairwise cost matrix $\boldC \in \R_+^{n \times m}$ between each unit $i$ and $j$: $\boldC_{ij} = c(x_i, x_j).$
There are several potential causal contrasts of interest over these populations but we focus on the sample average treatment effect (ATE):
Unfortunately, we cannot estimate Eq. (ref) since one of the potential outcomes is missing for each individual.
Instead, we need to use the information in the source population to get valid treatment effect estimates. That is, we desire $\E\left\{Y(z) \indicator(Z = z) \cdot w\right\} = \E\left\{Y(z)\right\},$ for some function $w$. A common way to do this is to use an importance sampling weight $w = \frac{d\alpha}{d\alpha_z}$, which is also known as the Radon-Nikodym derivative of $\alpha$ with respect to $\alpha_z$.
With a known $w$, we estimate Eq. (ref) as
and such that the weights sum to one in the treatment and control groups: $\sum_i w_i Z_i = \sum_i w_i (1 -Z_i) = 1$.
To identify these estimators, we need several assumptions, which we formalize below.
These standard assumptions allow us to use the observed data to estimate the desired treatment effects in the target sample. Assumption (ref) allows us to use the observed outcomes and not consider interference between units, while Assumption (ref) gives common support between treatment populations. With these conditions, we now turn to optimal transport and COT.
The popularity of optimal transport methods have exploded in recent years thanks to several recent theoretical and methodological advances Cuturi2013,Peyre2019, but the field dates back centuries. We frame our discussion in terms of empirical samples $\bolda_z$ and $\bolda$ to align with the rest of the paper but these quantities can be arbitrary samples for general optimal transport problems.
The original optimal transport problem formulated by Monge1781 involves finding optimal maps between distributions. Define such a map as a function $T: \mcX \mapsto \mcX$ and such that $\int_{\mcX } g(x) d\alpha_z = \int_{\mcX} g(T(x)) d\alpha$ for all measurable functions $g$. We denote the corresponding push-forward operator from $\alpha$ to $\alpha_z$ as $T_\# \alpha = \alpha_z$. The Monge1781 formulation of the optimal transport problem is then
where $T_\# \alpha = \alpha_z$. Unfortunately, this problem can be intractable to solve in practice since the mapping must be injective.
To alleviate this issue, the Kantorovich1942 formulation instead considers probabilistic assignments between distributions $\bolda_z$ and $\bolda$ using a transport matrix $\boldP$:
where $\boldU(\bolda_z, \bolda)$ is the set of joint distributions with margins $\bolda_z$ and $\bolda$. This metric is a proper distance that obeys the triangle inequality and metrizes the convergence in distribution, i.e. $\ot{\bolda_z}{\bolda}{} = 0 \iff \bolda_z = \bolda$ (Proposition 2.3, Peyre2019). When $c(x, x') = d_{\mcX}(x, x')^p$, as is the case for our setting, then Eq. (ref) is also known as the $p\,$-Wasserstein distance. Unfortunately, this problem is known to have a decaying convergence with increasing dimension Weed2017 and also to suffer from a large computational complexity Cuturi2013.
Conveniently, regularized optimal transport offers improved rates of asymptotic convergence Genevay2018, Mena2019 and computational speed Altschuler2017 by adding a convex penalty to the objective function:
Common penalties for $H_\lambda$ include an entropy penalty, $\lambda \boldP_{ij} \log \boldP_{ij}$ Cuturi2013, or an $L_2$ penalty, $\frac{\lambda}{2}\boldP_{ij}^2$ Blondel2018. The solutions to this problem converge to the solutions from Eq. (ref) as $\lambda \to 0$, while as $\lambda \to \infty$, the solutions put equal weight on every entry in $\boldP$.
To adjust for the fact that $\ot{\bolda_z}{\bolda}{\lambda} \neq 0$, genevay_learning_2018 introduced the Sinkhorn divergence for entropy penalized optimal transport:
This has the advantage that $S_\lambda (\bolda_z, \bolda) = 0 \iff \bolda_z = \bolda$ Feydy2018, while retaining the computational and theoretical advantages of regularized optimal transport.
Finally, we can still use Eqs. (ref) or (ref) to construct a map as in Eq. (ref). In finite samples, this function can be estimated from the Kantorovich formulation as
This mapping is alternatively known as the barycentric projection Peyre2019. For the squared-Euclidean cost, this map equals $\frac{1}{a_j}\sum_{i: Z_i = z} \boldP_{ij} X_i$, or the weighted mean of the observations in the sample who received treatment $Z=z$. For an $L_1$ cost, $T_{\bolda \mapsto \bolda_z}$ is the weighted median of the corresponding $X_i$. Under an $L_2$ cost, this map will also converge to the optimal Monge map provided one of underlying measures is continuous Ambrosio2005. With these general properties established, we now turn to our proposed method.
We define the COT problem as
where $\boldw_z$ is the empirical measure $\sum_{i: Z_i = z} \delta_{x_i} w_i $, $\Delta_n$ is the simplex with $n$ vertices, and $S_\lambda$ is defined in (ref). The COT weights will then be the weights that minimize $\cotobj{\bolda}{\lambda}$. In a slight abuse of notation, we have the following marginal distribution:
In addition to seeking distributional balance, a researcher may also know a set of functions that he or she thinks are important to balance a priori for valid causal estimates. These functions may include a hypothesized outcome model or the moments of the covariates. Define $B_k(\cdot): \mcX \mapsto \R$ for $k \in \{1,...,K\}$ as these $K$ functions of interest. We can then add an additional constraint to the problem in Eq. (ref) to approximately balance these important functions between samples:
Of course, other formulations of the problem are possible and we detail some of them in Appendix (ref) of the Supplementary Materials. However, we find that in practice the formulation in Eq. (ref) has the best performance in terms of bias and variance.
We now discuss the convergence of our weights to the distribution of interest. First, we define the importance sampling weights as $\breve{w}_i^\star = \frac{d\alpha(X_i)}{d\alpha_z(X_i)}$ and define the self-normalized importance sampling weights as $w_i^\star = \frac{1}{n}\breve{w}_i^\star/\sum_{i} \frac{1}{n}\breve{w}_i^\star.$ In our setting, $\breve{w}_i^\star =1/$ $\p(Z_i = z \given X_i) $. Further, let $\delta_{n}$ be the smallest value of the balancing function constraints at which the importance sampling weights satisfy the condition in Eq. (ref) for sample size $n$. We also rely on some additional assumptions to prove the convergence of the COT weights.
These assumptions enforce some regularity on the constituent parts of COT. First, the cost function must exist and be continuously differentiable. Second, the measures are either subgaussian or defined on subsets of the real numbers. Third, the penalty term $\lambda$ cannot go to zero too quickly if at all. And fourth, if using balancing constraints, there needs to be a value at which the importance sampling weights satisfy the constraints. With these assumptions, we have our first theorem.
This theorem says that the COT weights converge to the distribution of the target sample at a $\sqrt{n}$-rate, which also has implications for the efficiency of estimators based on COT, as we will see in the next section. A proof of this theorem is provided in Appendix (ref) of the Supplementary Materials.
Finally, Theorem (ref) also gives the following corollary.
The corollary follows as a consequence of the Radon-Nikodym Theorem and the fact that the Radon-Nikodym derivatives are almost surely unique.
For statistical inference, we turn our attention to the asymptotic distribution of Eq. (ref) and its variance. We assume the following conditions also hold.
This assumption has several important parts. We assume that the second moment of the outcome is finite and that the residual is $L_1$-integrable, which are not strong assumptions for real data. Then we require one of two additional conditions to hold. The first potential condition is that the convergence of $\boldw_{\text{COT}}$ to $\bolda$ occurs at a faster than $\sqrt{n}$-rate. We note that this is not actually that strong of an assumption in practice since the convergence to $\alpha$ happens at a $\sqrt{n}$-rate and COT is directly targeting the empirical distribution $\bolda$. Thus, we expect the convergence to the empirical distribution to be faster than $1/\sqrt{n}$, which is what we observe in practice. The second potential condition requires that the basis functions $B$ encompass the true conditional mean and that the empirical means of $B$ converge faster than $1/\sqrt{n}$. We observe in practice that the convergence of the basis functions is actually possible with relatively small sample sizes, making this assumption very plausible; however, outcome models are typically not known---though this can be ameliorated by using nonparametric models. With these conditions, we have our next theorem.
This result follows from the fact that the expansion of the bias $\hat{\tau} - \tau$ has the form of the semiparametrically efficient score function. We defer a proof to Appendix (ref) of the Supplementary Materials.
Theorem (ref) also gives us the following corollary.
In practice, one can check if Corollary (ref) holds by examining both the optimal transport distance between distributions $\boldw_{\text{COT}}$ and $\bolda$ as well as the balance of the hypothesized outcome models between samples of the covariate functions that determine the assumed outcome models . If the hypothesized outcome models are well-balanced, then there is little benefit to model augmentation. A manifestation of this phenomenon can be seen in the simulations in Section (ref) where adding in model augmentation does not change the estimates from using COT even for sample sizes as low as $500$.
One of the advantages of COT is that it provides a method to impute the missing potential outcomes, if so desired. Moreover, methods based on this estimator will also converge to the correct treatment effect, even when using weights from other methods.
We can construct a transportation matrix, $\boldP$, a posteriori for Eq. (ref) by solving $\ot{\boldw_z}{\bolda}{\lambda'},$ for any $\lambda' > 0$ using the appropriate weights from treatment group $Z=z$. Then the missing potential outcomes can be estimated by the barycentric projection in Eq. (ref): $\hat{Y}_j(z) = \argmin_\nu \sum_i \indicator(Z_i = z) c(Y_i, \nu) \boldP_{ij}$. In practice, we do not have to use the same cost function used to estimate $\boldP$ but using a squared-$L_2$ cost gives us the following theorem
There are several things to observe about this theorem and its assumptions. First, the assumptions on the cost and distribution allow us to connect the optimal transport solutions to the Monge maps of Eq. (ref), but also allow us to give a rate for our theorem. Second, using penalized optimal transport ensures that we average out the errors $Y - \mu_z$. Third, in a related manner, the Lipschitz continuity of the outcome means that this averaging out of the errors will still achieve good estimates of $\mu_z$. Fourth and finally, this theorem also suggests that $\mu_z$ is like a Monge map, $T$, between potential outcomes.
Unfortunately, these imputations are not necessarily useful by themselves.
Clearly, care must be taken when using these estimators. Proofs are located in Appendix (ref) and (ref) of the Supplementary Materials.
In this section, we turn to the practical considerations of optimizing the COT weights. Namely, we discuss the tuning of the hyperparameters and estimation of the weights.
\paragraph{Hyperparameter tuning.} Our goal is to select the hyperparameters $\lambda$ and $\delta$ so that we achieve the best distributional balance without over-fitting the current data. To do so, we propose a bootstrap based tuning procedure detailed in Algorithm (ref).
We justify this procedure in two ways. First, practitioners probably do not have an ideal weight penalty in mind based on subject matter knowledge. Second, because the COT weights target the Radon-Nikodym derivative, this tuning procedure will select the hyperparameters that lead to weights robust to sampling variation and better approximation of these population level quantities. We present an empirical examination of this tuning algorithm in Appendix (ref) of the Supplementary Materials that demonstrates its effectiveness at selecting the optimal $\lambda$.
\paragraph{Weight estimation.} Given the known complexity of estimating optimal transport distances, huling_energy_2020 raise the concern that methodology like COT will not be feasible. Fortunately, these concerns are addressed by using regularized optimal transport.
In our simulations, we find that estimating the COT weights only takes a few seconds for a 1000 observations. Eq. (ref) can be solved by alternating Sinkhorn divergence calculations in GeomLoss Feydy2018 and optimization steps on the weights. With balancing constraints, we use the Frank-Wolfe algorithm to optimize the weights Frank1956; without balancing constraints, we can use an LBFGS algorithm.
COT is also well-suited to the case where the target estimand is for a separate sample entirely. This is because the weights can be calibrated to target any arbitrary set of samples. The only additional assumption for convergence is that there is common support between distributions.
Finally, COT is easily extended to more than two treatments as long as the treatment values are discrete. This is because COT estimates weights separately for each treatment group.
The COT framework is actually related to several other methods in the literature and can be seen as an interpolation between all of them, as we detail below.
\paragraph{Synthetic control method} Abadie2003 first proposed the synthetic control method (SCM) as a way of performing counterfactual inference for a single treated unit, $j$. The objective function is \[\min_{w: w \trans \ones_n = 1, w_i \geq 0} \left \| \sum_i X_i w_i - X_j \right \|_2^2,\] which is the same objective as the Monge map in Eq. (ref) when $c$ is the squared-Euclidean distance and $T$ has the corresponding form of the barycentric projection in Eq. (ref). This means that SCM is actually estimating a version of the COT problem with $\lambda = 0$.
For a proof, see Appendix (ref) in the supplementary materials. One potential drawback of using SCM versus the formulation used for COT can be seen in the following simple example in Figure (ref). COT favors the nearest point while the SCM method utilizes the points further away, which could be a problem if the response surface looks like Figure (ref). To avoid this, SCM could incorporate a modified objective that directly models both the barycentric projection and distance between units as in Perrot2016.
\paragraph{Nearest neighbor matching} When done with replacement, nearest neighbor matching (NNM) is also a reformulation of the COT problem with $\lambda=0$: $\min_{\mathbf{P} \geq 0} \; \sum_{i: Z_i = z, j} \boldC_{ij}\boldP_{ij}$ subject to $\sum_{i,j} \boldP_{ij} \indicator(Z_i = z) = 1$ and $\boldP \trans \ones_{n} = \bolda.$ This will seek to find the unit $i$ that is closest in terms of $\boldC$ for each unit $j$ since this minimizes the total cost, which is the definition of NNM with replacement. Each observation will simply be weighted by the number of times it is matched, divided by the total sample size $n$.
NNM has some advantages and drawbacks relative to more general COT. Positive aspects of the method are that quick to estimate and also corresponds to an easily understood quantity of physical matching familiar to many researchers. However, the method can have poor convergence properties if $c \neq d_\mcX(\cdot, \cdot)^p$ where $p > d/2$ Fournier2015. This means that additional assumptions are necessary to ensure adequate convergence. Further, the weights will be given as rational numbers and as such, we would expect estimators based on them to have higher variance than weights without such constraints.
\paragraph{Optimal matching and MIP matching} Optimal Matching Rosenbaum1989 corresponds to COT weights with integer solutions and $\lambda = 0$: $\min_{\mathbf{P} \in \{0,1\}} \sum_{i: Z_i = z, \,j} \boldC_{ij}\boldP_{ij}$ subject to $L \leq \sum_{i} \boldP_{ij} \indicator(Z_i = z) \leq U, \; \forall j$ and $1 \leq L \leq U \leq n.$ We can turn this into Mixed Integer Program (MIP) matching by adding constraints: \[\left | \sum_{ij} \frac{B_k(X_i) \boldP_{ij}\indicator(Z_i = z)}{n M_j } - \frac{1}{n} \sum_{j} B_k(X_j) \right| \leq \delta_k, \, \, \forall k \in \{1,...,K\},\] where $M_j = \sum_{i} \boldP_{ij}$ is the number of matches for unit $j$ Zubizarreta2012. These methods also imply a re-weighting of treatment group $Z=z$ since the weights on element $i$ will be $n \inv \sum_j \boldP_{ij} /M_j$.
These methods have a similar flavor to NNM but with additional linear constraints. As such, they would likely share some of its benefits and drawbacks. As an advantage, these methods again yields matches which have an easy interpretation; however, as a disadvantage, the additional linear constraints will slow down the problem estimation. Moreover, the problem has additional tuning parameters $L$, $U$, and $\delta_k$ not present in NNM. Similar to NNM, we would again expect weights based on rational numbers to have higher variance.
\paragraph{Energy distance.} The Energy Distance (ED) is defined as $\mcE(\bolda_z, \bolda) = \frac{2}{n^2} \sum_{i:Z_i = z, \,j} d_\mcX(x_i, x_j)^p$ $- \frac{1}{n^2} d_\mcX(x_i, x_i)^p - \frac{1}{m^2} d_\mcX(x_j, x_j)^p,$ where $p \geq 1$. Then $S_\lambda(\bolda_z, \bolda)\to \frac{1}{2} \mcE(\bolda_z, \bolda)$ as $\lambda \to \infty$ Feydy2018. Thus, Energy Balancing Weights huling_energy_2020 are a special case of COT.
One advantage of this method is that there is no tuning parameter necessary to estimate the weights. However, we might assume that there would be an advantage to interpolating between all of these various methods. Indeed, in our experiments we find that the larger values of $\lambda$ do not approximate the true inverse propensity score as well as intermediate values. See Appendix (ref) in the Supplementary Materials for an empirical evaluation.
\paragraph{Mean Maximum Discrepancy.} Optimal transport is also related to the mean maximum discrepancy (MMD) through the ED. The MMD is equal to $\mcM = 0.5 \int_{\mcX \times \mcX} k(x, x') d \phi(x),$ for $\phi = \alpha - \beta$. For some reproducing kernel Hilbert space $k$ and for a distance $d_\mcX$ defined as $d_\mcX(x, x') = \frac{1}{2} k(x, x') + \frac{1}{2} k(x, x') - k(x, x')$, MMD is equivalent to the ED Feydy2018 and, therefore, to COT. To our knowledge, there has not been a proposed weighting method based on the MMD but we would expect it to have performance similar to that of the ED and COT with large values of $\lambda$.
To evaluate the finite sample performance of the proposed weighting methodology, we use the simulation study originally presented in Hainmueller2012. For each setting, we run 1000 experiments with a sample size of $n = 512$. The estimand of interest is the ATE.
For additional experiments examining the convergence of COT, its confidence interval coverage, and the efficacy of Algorithm (ref), see Appendix (ref) of the Supplementary Materials.
\paragraph{Study design.} We generate six covariates $X_1, ..., X_6$ from the following distributions
In this study, the last three covariates are mutually independent of each other and also of the first three covariates.
The treatment indicator is generated as $Z = \indicator(X_1 + 2 X_2 - 2 X_3 - X_4 - 0.5 X_5 + X_6 + \nu > 0),$ where $\nu$ is drawn from one of three distributions leading to different degrees of overlap: $\nu \sim \N(0,100)$ for high overlap, $\nu \sim g(\chi_5^2)$ for medium overlap, and $\N(0,30)$ for low overlap. The function $g$ in the medium-overlap setting gives the $\chi_5^2$ draws expectation 0.5 and variance 67.6. We expect the scenarios that will lead to the highest bias to be in the low-overlap setting where there is a strong separation between treatment groups and in the medium-overlap setting where the errors are leptokurtic.
Given the treatment indicator $Z $ and the covariates $X_1,...,X_6$, we draw the outcome $Y$ from $Y(0) = Y(1) = (X_1 + X_2 +X_5 )^2 + \eta,$ with $\eta \sim \N(0,1)$. There are two things to note about this outcome model. The first is that there is no effect of the treatment at the unit level and hence the ATE is 0. The second is that a linear outcome model should be biased.
\paragraph{Methods under examination.} We compare our methodology to several other weighting methods commonly used in the literature. The first method we consider is a logistic regression (GLM) using only first order terms. We also consider balancing methods such as the covariate balancing propensity score (CBPS) of Imai2014 and the stable balancing weights (SBW) of Zubizarreta2015 both targeting mean balance. Finally, we also utilize SCM and NNM.
For the COT weights, we include two variations both using an $L_2$ metric with standardized covariates. The first only balances the joint distribution (no constraints or “none”) and the second demonstrates basis function balancing by targeting the joint distribution as well as mean balance (“means”).
In our simulations, we consider three estimators to target the ATE. The first is known as the H\'{a}jek estimator Hajek1988 and is simply a weighted mean with sum to one weights as in Eq. (ref). The second is an augmented or doubly robust estimator of Robins1994, and the third is a weighted least squares estimator---both only including linear terms of the covariates. We do not include a barycentric projection estimator like in Eq. (ref) since under an $L_2$ metric it gives the same result as the H\'ajek estimator (see Proposition (ref)).
We now turn our attention to the results. Due to its nonparametric to semiparametric formulation, we expect COT to do better than other methods when the true propensity score model diverges from a logistic regression---e.g., in the medium-overlap scenario.
Indeed, the COT methods have the lowest RMSE across all overlap scenarios, as we can see in Table (ref). Further, across the medium and low overlap settings, COT has the lowest bias as well; in the high overlap scenario, there is a negligible difference between COT, SBW, and GLM.
We should also note that COT gives estimates that do not vary between estimators. This is because the COT already balances the basis functions used in the augmented estimator and weighted least squares. Therefore, there is no difference between running an outcome regression utilizing linear terms of the covariates and the H\'{a}jek estimator. This is a function of the fact that by balancing distributions of the covariates, COT will also balance functions of the covariates. Similar phenomenon can also be observed with SBW.
In this section, we apply our methodology to a real data set. There is growing interest in the literature to utilize libraries of randomized control trials (RCTs) to evaluate new interventions, the idea being that running new RCTs to evaluate every new intervention is expensive and time consuming Schmidli2020. These studies, alternatively called externally controlled trials or synthetic control group trials, compare a set of study subjects receiving a treatment to a group of individuals external to the trial at hand who did not receive the intervention of interest. The participants used for the control group can be taken from a variety of sources such as an observational study, electronic medical records, or from historical clinical trial data Davi2020. To demonstrate this in practice, we present an analysis utilizing data originally from a multi-site RCT discussed by Blum2010.
The original study was a double-blind, non-inferiority trial that exposed 31,055 women to prophylactic oxytocin during labor at five hospitals across Burkina Faso, Egypt, Turkey, and Vietnam. The 807 women in this group with uncontrolled blood-loss after delivery---a condition known as post-partum hemorrhage or PPH---were then randomized to receive either 800 milligrams misoprostol (treatment condition) or 40 international units oxytocin (control condition). There were 407 and 402 women in each treatment group, respectively. The primary outcome for the study was whether blood loss was controlled within 20 minutes after PPH diagnosis. The authors measure several important confounders like maternal age, blood loss at treatment, whether cord traction was maintained, maternal hemoglobin, whether the mother is currently married, whether the cord was clamped early, fetal gestational age, whether labor was augmented, whether labor was induced, maternal education, number of previous live births, whether the placenta was delivered prior to hemorrhage, and whether a uterine massage was given.
We modify the study in a couple of ways to make it more similar to an externally controlled trial. For each site, we separate the paired treatment groups and attempt to estimate effects using the units from other sites---e.g., for the misoprostol group in Egypt we remove the oxytocin group from Egypt and attempt to estimate a causal effect using the oxytocin groups from the other sites. In this manner, we generate effect estimates and confidence intervals for each treatment group at each site.
To estimate the treatment effects, we use COT with hyperparameter tuning as in Algorithm (ref) and the squared-Euclidean distance as the cost function. We also compare COT to GLM, CBPS, SBW, SCM, and NNM. For GLM, CBPS, and SBW we utilize all first and second covariate moments. Finally, our estimator is that of Eq. (ref).
For our estimates to be valid, we require that there be no unmeasured confounding but also that the estimates are “transportable,” i.e. we are able to take estimates from the other hospitals in the external group and “transport” them to site of interest. This requires that conditional on the observed covariates there are no other variables that can effect the outcome and treatment indicator ($d$-separation holds, Pearl2013).
Amazingly, the three optimal transport flavored methods are able to achieve estimates close to the original treatment effects on average. In Figure (ref), we can see that COT, NNM, and SCM do the best job in terms of average bias; moreover, these are the only methods that have estimates inside the original confidence interval. Of these, we note that COT has the least overall bias across all treatment groups and sites.
In terms of inference, COT also performs best across all sites and treatment groups. Table (ref) displays how well the calculated confidence intervals cover the original treatment effects and also if the calculated estimates are inside the original confidence interval. In both cases, COT has the highest percentage of confidence intervals covering the original treatment effect and estimates inside the original confidence interval at $60\%$ each, respectively. And while NNM had a good overall bias, only 20% of its confidence intervals covered the true treatment effect and its estimates were inside of only 20% of the original confidence intervals.
We have described a new tool for the estimation of causal effects in observational studies: Causal Optimal Transport. This method allows for checks of distributional overlap and model-free weight estimation that is semiparametrically efficient. We also showed how several other methods for causal inference are closely related to COT and that COT can be seen as an interpolation between these methods. In our simulation study, we demonstrated that this methodology performs well even when both the outcome and propensity score models are misspecified. Compared to other common weighting methods, COT generally has lower bias and lower root mean-squared error under model misspecification.
There are several areas for future research. First, the sensitivity of COT to the choice of cost function remains to be elucidated. Second, selecting covariates through typical model selection frameworks such as an $L_1$ penalized regression is not obvious given that COT does not generate clear predictive models, though this may not matter given the nonparametric nature of the weights. Third, more work needs to be done to extend this framework to time series data, but further connections to SCM may offer a way forward.
\if00 {
The author would like to thank Claire Chaumont, Gang Liu, Aar\'on Sonabend, Lorenzo Trippa, and Jos\'e Zubizarreta for helpful comments and feedback on an earlier version of this manuscript. This research was funded by generous support from NIH grant 5T32CA009337-40, the Department of Biostatistics at the Harvard T.H. Chan School of Public Health, and the David Geffen Scholarship from the David Geffen School of Medicine at UCLA. } \fi
The Supplementary Materials contain proofs of the theorems in the paper, further empirical studies, an additional case study, and proofs for other formulations of COT. The last two s Proofs are found in Appendix (ref), and empirical studies of the convergence, confidence interval coverage, and efficacy of the tuning algorithm are located in Appendix (ref). Appendix (ref) is the additional case study and Appendix (ref) are the proofs for other formulations of COT.
\printbibliography