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.
68,057 characters · 17 sections · 72 citation commands
An optimal transport approach to estimating causal effects via nonlinear difference-in-differences
Difference-in-differences estimators are among the most widely-used approaches to estimate causal effects of discrete treatment interventions in observational data abadie2005semiparametric, lechner2011did, imbens2015causal, roth2022s. The classical setting we considered in this article is one in which the researcher observes the outcomes of interest in a group of treated units and a corresponding control group of untreated units over time, i.e. a binary treatment setting.\footnote{The proposed method can straightforwardly be extended to finitely many treatments. We do not address the question of continuous treatments in this article.} At one fixed point in time the treatment group undergoes an intervention and stays treated thereafter. The fundamental problem consists of isolating the causal effect of this intervention from the existing natural trend the treatment group undergoes irrespectively. The method of difference-in-differences achieves this under the assumption that the units in the control group have the same fundamental trend as the treated units would have without intervention. Under this assumption, the difference of the average outcome of the treated group between pre- and post-intervention periods net the difference of the control groups average outcomes pre-and post intervention manages to isolate the true average causal effect of the intervention.
The arguably two most fundamental methods for estimating causal effects in this setting are the classical method of Difference-in-Differences (DiD) abadie2005semiparametric, lechner2011did, wing2018designing as described above and its generalization beyond average and aggregate effects, the Changes-in-Changes (CiC) estimator athey2006identification.
DiD is a linear approach designed to estimate average (or aggregate) causal effects abadie2005semiparametric, heckman1990varieties over all units within a group. This method can be proved to correctly identify counterfactuals under the “parallel trends” assumption abadie2005semiparametric, roth2020parallel, ryan2019now under which the average natural trend is assumed to be the same across both the control and treatment group. This idea has been used in many areas of science where capturing the average treatment effect is sufficient to formulate informative causal conclusions. Recent applications include quantifying the effect of public health measures in response to COVID-19 bacon2020using and estimating how irrigation farmers adapt their watering to fixed usage limits drysdale2018adaptation.
The classical linear DiD estimator only provides estimates of average and aggregate treatment effects. This can be limiting in many settings where treatment heterogeneity is important, i.e., where different units can react differently to the same level of intervention. The CiC estimator athey2006identification addresses this issue for univariate outcomes. It extends the fundamental idea of the DiD estimator by estimating changes in the entire probability distribution of the respective units within a group. This permits the estimation of the entire counterfactual law of the treated units had they not received the treatment, hence allowing for a general form of heterogeneity in their response to treatment. The cost is a slightly stricter assumption on the evolution of the unobservable distribution as we specify below.
In its current form the CiC estimator is only applicable in settings with univariate outcomes, as it relies heavily on the definition of quantile functions. This can be restrictive for modern applications where the outcomes of interest are often multivariate. Examples range from A/B testing in digital marketing campaigns where outcomes are a combination of features such as click-through rate, time-per-page, etc., to measuring intervention effects such as gene knock down on a population of cells measured in high-dimensional gene space. In principle, the CiC estimator may be extended to higher dimensions through tensorization, which estimates treatment effects independently for each coordinate, but this solution fails to capture correlations that are often important in causal discovery. We demonstrate a simple linear setting with synthetic data where this is the case in Section (ref) below. \\
{\bf Our contribution.} We recast the the CiC estimator using tools from the theory of optimal transportation. This perspective allows us to extend this methodology to handle multivariate observables after introducing novel structural assumptions on the natural trend underlying both treatment and control groups. In particular, we note that the causal model proposed in athey2006identification readily implies that both populations, treatment and control, evolve between pre- and post-treatment periods via optimal transport maps (Theorem (ref)); our methodology relies on a natural higher-dimensional extension of this via optimal transport theory. In fact, we show that if the corresponding functions mapping the unobservables to the potential outcomes are cyclically comonotone then the counterfactual distributions and causal effects can be estimated from data. In particular, this allows for the unobservable variables to be of the same dimension as the dimension of the outcome variables and not just univariate.
The counterfactual distributions can be estimated consistently from data using recent results from statistical optimal transport and implemented efficiently using recent advances in computational optimal transport peyre2019computational. We demonstrate the benefits of this extension by comparing it to the simple multivariate extension of the CiC estimator both on artificial and real data. As illustrated in Section (ref), a tensorized version of CiC can be inconsistent and even estimate opposite correlation structures in the multivariate distribution of counterfactual outcomes. We also revisit the classical dataset of card1994minimum that has sparked an intense debate about the effects of raising the minimum wage on employment. Our ability to jointly handle the number of part-time and full-time workers as a bivariate outcome reveals an interesting substitution effect from the perspective of labor economics: when the minimum wage is raised, restaurants increase full-time labor and decrease part-time labor suggesting a substitution effect.
We recall in this section the general causal model on which the classical difference-in-differences and changes-in-changes estimators are built, and which our model extends. See Figure (ref) for illustration.
We first formalize a stochastic model for an experimental design in which two groups are measured before and after some intervention. Sampled from a larger population, each unit $i$ in a realized experiment is characterized by four vector potential outcomes, their level of treatment received, and a set of time points at which potential outcomes were observed. This model includes randomized control trials with non-compliance and general observational studies. Adapting the notation of athey2006identification, we define an indicator random variable $G_i$ for a unit's adopted treatment arm as well as random vectors $Y_{i;C,0}$ and $Y_{i;C,1}$ to model unit $i$'s potential control group pre- and post-intervention outcomes observable when $G_i = 0$; likewise $Y_{i;T,0}$ and $Y_{i;T,1}$ are the unit's observable potential outcomes when $G_i = 1.$ These potential outcomes are assumed to lie in a subspace of $\R^d,$ and each unit has indicator random variables $T_{i,0}$ and $T_{i,1}$ denoting whether an outcome was observed in each study period.
We assume that each potential outcome vector is generated by a deterministic function of a latent random vector also in $\R^d.$ While any given unit's latent vector may change over time, the population distribution of continuous latent variables is time-invariant. Furthermore, these latent distributions can differ arbitrarily between control and treatment arms. We denote this latent distribution $\nu$ for controls and $\nu^\star$ for treated. Therefore we write for control units the latent random vectors $U_{i,t} \stackrel{d}{\sim} \nu, t=0,1$ and treated $U_{i,t} \stackrel{d}{\sim} \nu^\star$. We can think of $U_{i,t}$ as capturing all of the intrinsic and unobservable characteristics of unit $i$, such as an employee's skills and motivation, which influence the outcome of interest at time $t$, for instance their salary. This assumption may be less appropriate in settings where external shocks could change the latent distributions over time. For instance in a study measuring individuals' opinion of the economy, a latent distribution of intrinsic perception reasonably shifts over time. Our model can accommodate exogenous noise through the assumption that latent $U_{i,t}$ is the additive perturbation of a unit's time-invariant baseline latent state $\mu_i$ by a stationary noise process $\varepsilon_{i,t}$. Correlation between these two terms is allowed as long as the distribution of realized latent variables is time-invariant within each treatment arm. These latent variable assumptions are general enough to include many forms of mechanistic structure researchers may know about their population.
As Figure 1 illustrates, three global production functions $h_0$, $h_1$, and $h_1^\star$ characterize the evolution of both treatment arms over time. At $t=0$ the treatment arms are assumed to have identical factors influencing the translation of a latent variable to observed potential outcome and thus share $h_0$. For controls $h_1$ describes any changes due the natural trend over evolving time. For treated units, $h_1^\star$ captures both drift over time and the causal effect of interest. If $h_1 \circ (h_0^{-1})$ exists, it describes exactly how any given control group unit evolves over time with a fixed latent state. We can therefore obtain a counterfactual outcome distribution with this function under an assumption that the treatment arm units would have evolved in the absence of treatment just as controls do. It is notable that this proposed model does not rely on additional information from observed covariates. As such, it does not require classical unconfoundedness assumptions. However, the question about the assignment mechanism in these difference-in-differences setups is a delicate one, see for instance ghanem2022selection, marx2022parallel.
As units in the same treatment arm share a latent variable distribution, we define random vectors with the population distribution such that $Y_{i;C,0} \stackrel{d}{\sim} Y_{C,0}$ for all units such that $G_i = 0$. The random vectors $Y_{C,1}$, $Y_{T,0}$, and $Y_{T,1}$ are defined similarly for observations adopting the appropriate treatment arms. Each population-distributed random variable of potential outcomes induces a measure or law, which we will utilize in application of results from optimal transport theory. We adopt the notation $\Prob_{Y_{G,T}}$ for this measure. The action of each production function can be described as a pushforward of measure. For example we write $\Prob_{Y_{C,0}} = h_{0\#} \nu$ and $\Prob_{Y_{T,1}} = h_{1\#}^\star \nu^\star$. We make an assumption that both $\Prob_{Y_{C,0}}$ and $\Prob_{Y_{T,0}}$ are absolutely continuous with respect to Lebesgue measure. In Section (ref) we discuss possible changes to our assumptions which would allow our methodology to handle discrete outcomes at the cost of point identification of the counterfactual distribution. We note that the discrete case in general does not allow for point-identification, which is already evident in the univariate setting, where athey2006identification provide bounds on discrete counterfactuals. Our optimal transport framework illuminates this issue even further.
Additionally assuming that $h_0$ is invertible, we have $\nu=(h_0^{-1})_{\#} Y_{C,0}$ so that
In the sequel we will propose a class of production functions such that the multivariate counterfactual distribution for the treated, distributed as a random variable ${Y}_{T,1}^\dagger$, can be identified as
A consistent estimator for $\sd$ and hence this counterfactual pushforward is presented utilizing optimal transport. Figure (ref) shows that the production function $h_1^\star$, which generates post-intervention outcomes for the treated, does not need to be specified to recover the map $\sd$ and its counterfactual pushforward $Y_{T,1}^{\dagger}$. As discussed in Section (ref), further assumptions on $h_1^\star$ enable the treatment effect map $\sT$ to be consistently estimated and allow for the estimation of non-linear treatment effects.
Given two measures $\Prob_{Y_{C,0}}, \Prob_{Y_{C,1}}$ which can be estimated from control group data pre- and post-intervention, there exist infinitely many maps $\sd$ such that $\Prob_{Y_{C,1}} = \sd_{\#} \Prob_{Y_{C,0}}$ (ref) holds. To identify a unique such map $\sd$, one must make additional assumptions on the causal model.
Recall that the source of randomness in the proposed causal model are latent random vectors $U_{i,t}$ which can be thought of as summarizing a unit's intrinsic characteristics. If we not only assume that the distribution $\nu$ of control group latent vectors is time-invariant but further that the latent variables are constant over time (i.e., $U_{i,0}=U_{i,1}$ for all control units $i$), then we have $Y_{i;C,1}=\sd(Y_{i;C,0})$; this map can be estimated via multivariate regression from independent copies of the pair $(Y_{C,0}, Y_{C,1})$. In many practical applications, pairs of potential outcomes $(Y_{i;C,0}, Y_{i;C,1})$ are observed, but in general the assumption that $U_{i,0}=U_{i,1}$ for all control units is more tenuous. Moreover, this setup leaves out an important class of problems where the units in the control group can differ between $t=0$ and $t=1$. A canonical example arises in genomics where in single-cell RNA-Seq data each unit corresponds to a cell. Measuring the gene expression levels of a cell is a destructive process, and as a result a given cell may be only measured once. Therefore, coupled potential outcomes of the form $(Y_{i;C,0}, Y_{i;T,1})$ are not available. Instead, one observes independent copies of $Y_{C,0}$ as well as independent copies of $Y_{C,1}$. While difficult, this problem can be solved at the cost of additional assumptions on the natural trend $\sd$. This problem is sometimes referred to as uncoupled regression and it arises in various applications, far beyond the genomics example described above RigWee19, balabdaoui2021unlinked, slawski2022permuted.
In a setting with uncoupled scalar data, the map $\sd$ is unique when it is monotone with known monotonicity (either increasing or decreasing). The CiC estimator enforces this by assuming both production functions, and therefore $h_1 \circ h_0^{-1}$, are monotone increasing. In fact, the unique monotone increasing function $\sd$ such that $\Prob_{Y_{C,1}}=\sd_\# \Prob_{Y_{C,0}}$ is given by $F_1^{-1}\circ F_0$ where $F_t$ denotes the cumulative distribution function (CDF) of $Y_{C,t}$ for $t=0,1$. Estimation of $\sd$ from data and plugging that estimator into Equation (ref) to yield an estimated counterfactual measure is the key idea of the CiC estimator. As previously noted, this estimator can be extended to multivariate outcomes through tensorization, which applies it independently to each coordinate. Joint structure between the coordinates is incompatible with the scalar distribution and quantile functions used in the estimator, as well as with the assumption of scalar production functions.
Recently in the univariate setting, other approaches for modeling the control group time trend have been proposed which generalize the “parallel trends” assumption to allow for general heterogeneity. These models in general make stronger or less interpretable assumptions than the monotone production functions of athey2006identification. callaway2019quantile directly assume that the copulas between $\Prob_{Y_{C,0}}$ and $\Prob_{Y_{T,0}}$ as well as $\Prob_{Y_{T,0}}$ and $\Prob_{{Y}_{T,1}^\dagger}$ are the same; roth2020parallel assume that the pointwise differences between the corresponding cumulative distribution functions are equal, i.e. that $F_1(x)-F_0(x)={F}_1^\dagger(x)-F_0^\star(x)$ for all $x\in\mathbb{R}$; bonhomme2011recovering restrict the heterogeneity of the model to be additively separable and assume that the pointwise differences between the corresponding logarithms of the characteristic functions are equal. In the next section we introduce cyclical comonotonicity, a natural extension of the monotonicity requirement in athey2006identification to multivariate settings, and hence the multivariate notion of parallel trends closest to the framework of that work.
In the following section, we show that under the CiC model the univariate natural trend map $\sd$ is unique and therefore identifiable as an optimal transport map between $\Prob_{Y_{C,0}}$ and $\Prob_{Y_{C,1}}$. We leverage this insight to identify unique natural trend maps in higher dimensions under the same model. In particular, we propose an assumption of cyclically comonotonicity, a multivariate notion of monotonicity between the production functions, which extends the CiC estimator's assumption of monotone increasing production functions.
Identifiability of the causal model presented in the previous section hinges on the uniqueness of the monotone increasing map $\sd$ such that $\sd_{\#} \Prob_{Y_{C,0}} = \Prob_{Y_{C,1}}$. This uniqueness follows from a fundamental result in the theory of optimal transportation known as Brenier's theorem.
For $d=1$, a map $T$ is the gradient of a convex function if and only if it is non-decreasing. The quantile map $F_{Y_{C,1}}^{-1}\circ F_{Y_{C,0}}$ has the appropriate domain and range measures and is non-decreasing by the properties of CDFs; therefore it is the unique optimal transport solution.
Theorem (ref) demonstrates that the structural properties of the univariate CiC model can be extended to higher dimensions by leveraging the theory of optimal transportation. In particular, the Brenier map is unique and hence identifiable between pairs of sufficiently regular measures in $\mathbb{R}^d$. We therefore introduce regularity assumptions on the production functions $h_{\cdot,\cdot}$ in the model which imply that the optimal drift is the Brenier map between $\Prob_{Y_{C,0}}$ and $\Prob_{Y_{C,1}}$. The counterfactual distribution $\Prob_{Y_{T,1}^\dagger}$ in Figure (ref) can be identified by first estimating the unique natural trend $\sd$ between the observable measures $\Prob_{Y_{C,0}}$ and $\Prob_{Y_{C,1}}$ in the control group and then applying that transformation to the treatment group pre-intervention potential outcome measure $\Prob_{Y_{T,0}}.$
For our main result, we assume that all four the potential outcome measures $\Prob_{Y_{\cdot,\cdot}}$ and the counterfactual measure $\Prob_{Y_{T,1}^\dagger}$ are absolutely continuous with respect to the Lebesgue measure. This assumption for pre-intervention measures is needed to apply Theorem (ref), and for the post-intervention measures it helps enable identification. To satisfy technical assumptions on the regularity of the optimal transport problem that often arise in classical results caffarelli2000contraction, we enforce through Assumption (ref) that measures are supported on convex sets. Furthermore, if the support $K_0^\star$ of $\Prob_{Y_{T,0}}$ is not contained in the support $K_0$ of $\Prob_{Y_{C,0}}$, one cannot infer the counterfactual distribution $\Prob_{Y_{T,1}^\dagger}$ at the parts of $K_0^\star$ that fall outside $K_0$ without extrapolation via structural assumptions. This restriction on the supports is a standard assumption that prevents the need for extrapolation athey2006identification. We summarize our structural assumptions as follows:
When $\sd$ is an optimal transport map, Theorem (ref) implies $\sd$ is the gradient of a convex function. Recall that in higher dimensions, gradients of convex functions from $\mathbb{R}^d\to\mathbb{R}$ are cyclically monotone rockafellar1997convex. A map $T: \mathbb{R}^d\to\mathbb{R}^d$ is cyclically monotone if for any positive integer $m$ and any cycle $u_1, \ldots, u_m, u_{m+1} \equiv u_1$ in its domain, it holds
Therefore, the natural trend map is identifiable when it is cyclically monotone. For $m=2$, Equation (ref) reduces to the simple notion of multivariate monotonicity defined as
Furthermore, in the univariate case these two definitions are equivalent; in the CiC model this implies that the monotone map $h_1 \circ h_0^{-1}$ is also cyclically monotone, which explains why the assumption of monotonicity of the functions $h$ in the unobservable is sufficient.
In the context of our latent model, the monotonicity assumption on the production functions for the CiC estimator implies that larger values of the unobservable latent variable correspond to strictly larger potential outcomes. This is a common assumption in economic models matzkin2003nonparametric, chernozhukov2005iv, imbens2007nonadditive, but can be restrictive in some settings. For instance, consider a classroom experiment where the treatment variable is class size and the outcomes of interest are standardized test scores before the intervention and 5 years after. Monotone production functions imply that individuals with higher unobserved ability $U$ will receive strictly higher test scores. Specifically, the assumption of being able to rank individuals in terms of one univariate unobservable $U$ is a strong restriction.
A careful inspection of the data generating process reveals that monotone assumptions on the production functions allow the latent variable to be entirely abstracted away. This actually follows from the weaker notion of comonotonicity such that $h_0,h_1:\R^d \to \R^d$ are comonotone if $\langle h_0(x)-h_0(y),h_1(x)-h_1(y)\rangle\geq 0$ for all $x,y \in \R^d$. Note that with an identity production function comonotonicity reduces to classical monotonicity. To get a sense of the strength of this relaxation, consider univariate differentiable functions. In this case comonotonicity implies locally that their derivatives $h_0'\cdot h_1'\ge 0$ have the same sign and furthermore imposes additional global constraints. As a concrete example, if $h_0$ is a polynomial production function and $h_1$ its pointwise scaling by some $\gamma > 0$, then this pair of functions is comonotone because they have the same sign between all zeros and hence the same signed difference between any pair of points. This example emphasizes that $h_0$ and $h_1$ need not be individually monotone themselves.
We have now seen that that the monotone production functions assumption in athey2006identification implies that the natural trend is cyclically monotone, which enables identification of the counterfactual. Furthermore, the production functions are comonotone, which allows the model to include general latent space distributions. However, these two functional assumptions only collapse to monotonicity in the scalar case. To extend their model to higher dimensions, we propose an assumption of cyclically comonotone production functions which achieves both these desired properties.
Just as cyclical monotonicity collapses to monotonicity when $m=2$, cyclical comonotonicity collapses to comonotonicity in that case. Whenever $h_0$ has an inverse $h_0^{-1}$, condition (ref) implies the map $\sd=h_1 \circ h_0^{-1}$ is also cyclically monotone. This result is shown in the proof of Theorem (ref), which we will state shortly. By Theorem (ref), $\sd$ is the unique Brenier map such that $\Prob_{Y_{C,1}}=\sd_{\#}\Prob_{Y_{C,0}}$.
Cyclical comonotonicity is one of several concepts to extend comonotonicity to higher dimensions. One set of extensions imposes a total ordering on the outcome space; however, these imply that the copula of outcomes pre- and post-intervention are identical, a strong assumption on the causal model. Note that the tensorized CiC assumes this shared copula structure. More recently, optimal transportation theory has been used to propose other multivariate notions of comonotonicity. A particular contribution is $\mu$-comonotonicity from galichon2010comonotonic. Similar to our model, the authors assume a latent variable and production functions mapping it to potential outcomes. Both production functions $h_0$ and $h_1$ are assumed to have optimal transport structure, as they lie in the graph of a gradient of a convex function. In general, the composition of optimal transport maps, such as $h_1 \circ h_0^{-1}$, is not an optimal transport map, except in special circumstances boissard2015distribution. puccetti2010multi provide an overview of these ideas.
Cyclical monotonicity is a natural extension of monotonicity to multivariate settings with many uses in economic theory and mechanism design ashlagi2010monotonicity, rochet1987necessary, econometrics shi2018estimating, chernozhukov2021identification, and statistics ruschendorf1996c. Cyclical comonotonicity is even weaker than this, as it does not imply that each production function itself is cyclically monotone, which makes it a natural and weak extension of the monotonicity assumption in the univariate setting. In particular, it is weaker than the assumption that individuals can be ranked in terms of a univariate unobservable variable (ability in the classroom example) and that this ranking persists in the observed outcome distribution. We argue that this makes it a more natural assumption than the monotonicity assumption of the CiC estimator in practical settings. Examples of cyclic monotonicity in causal inference have been highlighted in chernozhukov2021identification, shi2018estimating, and gunsilius2022condition, among many others. More fundamentally, rochet1987necessary shows that the maximizer of quasi-linear utility functions is cyclically monotone. For a concise survey, we refer to galichon2017survey.
From a practical perspective, cyclical comonotonicity is an important weakening of cyclic monotonicity in that it allows for unobservables of general dimension in practice. Allowing for multivariate unobserved heterogeneity in general causal models has been a long-standing challenge: in most models, the unobservable $U_{i,t}$ is required to be univariate athey2006identification, matzkin2003nonparametric, imbens2007nonadditive. By relying on the concept of cyclical comonotonicity, we can relax this assumption by only requiring that the unobservable is of the same dimension as the observable variable. Note that even this can be relaxed by using recent advances in optimal transport theory mccann2020optimal. The cost of this relaxation is a more complicated model and definitions from geometric measure theory like the area and coarea formula, which is why we focus on the simpler setting of unobservables of the same dimension in this article.
\paragraph*{Discrete outcome distributions} Mirroring the univariate setting athey2006identification, our methodology does not guarantee a point identified-counterfactual distribution when outcomes are discrete-valued. Theorem (ref) guarantees a unique control group natural trend with optimal transport structure when the measure of pre-intervention outcomes is absolutely continuous with respect to the Lebesgue measure. Indeed, unique optimal transport maps between discrete measures need not exist. For illustration, consider transporting a measure distributed on vertices of the unit square with equiprobability at each main diagonal vertex to a one with equiprobability at the anti-diagonal vertices. With respect to squared distance, it is equivalent to transport the mass between vertical or horizontal pairs of vertices. Furthermore, it is admissible to split mass and transport $0.25$ probability mass from each main diagonal vertex to each anti-diagonal vertex. This solution with fractional structure is an example of a non-deterministic optimal transport plan (cf. deterministic optimal transport maps), which can be interpreted as a probabilistic assignment rule. Discrete-valued optimal transport problems often have a minimum cost achieved by multiple probabilistic plans but no deterministic maps. This reality implies that point-identification in the discrete case need not be possible in general without restrictive further assumptions.
This optimal transport view hence also explains the lack of point-identification in the univariate setting, matching the conclusions of athey2006identification. The extension of CiC to discrete outcomes in that setting provides bounds on the counterfactual outcome distribution due to the same challenges described above. Their estimator for that case is heavily dependent on notions of quantile functions without a direct multivariate extension. Identification only becomes possible with additional assumptions about conditional independence and exogenous covariates, which the authors admit are restrictive. In the discrete multivariate setting, similar partial identification results are possible through carefully application of recent developments in optimal transport theory. For example auricchio2022structure bound the supremum cost of moving any single point mass in $\ell^p$-cost transport problems. With one possible counterfactual distribution identified by a solution to the minimization, uniform bounds can then be found by applying the aforementioned bound on the possible pointwise deviation. Point-identification can also be reestablished in a similar manner to the univariate setting by introducing additional information from covariates under strong assumptions. The focus of this article is to introduce an optimal transport framework for point-identification, which is why we focus on the case of absolutely continuous measures throughout.
Based on the extensions of the causal model introduced in the previous sections, we can formally state the identification result for the counterfactual distribution of the treated units had they not received treatment.
The proof of Theorem (ref) is relegated to Appendix (ref). Another potential quantity of interest is the actual counterfactual random variable $Y_{T,1}^\dagger$. Identifying this quantity requires more structure in the model. Indeed, in the setup of the causal model represented in Figure (ref), we only have that $\Prob_{Y_{C,1}}=\sd_\# \Prob_{Y_{C,0}}$ but not necessarily that $Y_{C,1}=\sd (Y_{C,0})$ because we do not assume\footnote{We note in passing that if we were to assume that $U_{i,0}=U_{i,1}$ for all units, we could, in fact, test the assumption that the natural trend $\sd$, and hence the production functions $h_0,h_1$, are monotone increasing from independent copies of paired data of the form $(X_0,Y_1)$. We would do so by checking if there exists an increasing function $\sd$ that indeed interpolates these points. In practice, such a function seldom exists and allowing a time-variable $U_{i,t}$ accounts for deviations from this ideal situation.} that $U_{i,0}=U_{i,1}$ for all units. However, assuming that both $h_1^\dagger$ and $h_1$ are cyclically comonotone allows to identify $Y_{T,1}^\dagger$, which follows from the fact that $T:K_1^\star\to K_1^\dagger$ will be identified for the fixed unobservable $U_{i,1}$ for treated units.
The assumption that $h_1$ and $h_1^\star$ are cyclically comonotone can be hard to satisfy in practice, which is why the main parameter of interest is the counterfactual law $\Prob_{Y_{T,1}^\dagger}$ and not the random variable $Y_{T,1}^\dagger$. It is worth noting that having cyclically comonotone pairs $h_1^\star$ and $h_1$ as well as $h_0$ and $h_1$ does not necessarily imply $h_1^\star$ and $h_0$ are cyclically comonotone.
In general, the assumptions of our model match those in athey2006identification when outcomes are scalar and extend their structure in higher dimensional cases. Our production function setup illustrated in Figure (ref) exactly matches that of CiC. In the scalar case, our assumption of comonotone production functions reduces to monotone production functions; although this is weaker than CiC's requirement that $h_0$ and $h_1$ are strictly monotone, the requirement of Theorem (ref) that $h_0$ has a well-defined inverse restricts that pre-intervention production function to be strictly monotone as well. In one dimension, a strictly monotone and a weakly monotone production function pair remains cyclically comonotone, so our model does provides a relaxation of the control group outcome generation process. Here it is the assumption that $h_0$ has an inverse which forces the production functions to be themselves monotone, as contrasted by our previous example of two positively proportional polynomials being comonotone. In higher dimensions, however, the invertibility of $h_0$ does not imply both $h_0$ and $h_1$ are themselves cyclically monotone. If this were the case, then both production functions would be optimal transport maps, which can be seen by considering cyclical monotonicity as the special case of cyclical comonotonicity when $h_0$ is the identity. This is precisely the notion $\mu$-comonotonicity introduced by galichon2010comonotonic.
For the latent distributions, our extension retains the time-invariance within groups assumption of CiC. While we do not make the same assumption of nested latent supports $\nu^\star \subset \nu$ as athey2006identification, our assumption that $K_0^\star \subset K_0$ is equivalent. To see why note that both treatment arm share the same pre-intervention production function $h_0$. Our assumption that the support sets are convex is an additional requirement not needed for CiC but necessary to apply theoretical results about multivariate optimal transport. athey2006identification allow for an additional case with discrete latent distributions, but our assumption that $h_0$ is invertible with a continuous co-domain implies that the latent variable measures must be continuous.
\paragraph*{Examples of functionals that can be analyzed with the proposed method} One of the advantages of the method is to deal with multivariate outcome distributions. Many functionals of interest in different areas of the (social) sciences are based on the information of the joint distribution of several outcomes of interest, from social inequality decancq2013weights, duclos2006robust, bourguignon2003measurement to risk management ekeland2012comonotonic, embrechts2002correlation to decision theory hallin2022measure, to name a few.
More generally, the proposed estimator provides consistent estimates of functionals that fundamentally rely on the copula structure of the outcome distributions. Examples of this are multivariate stochastic dominance garcia2019review, multivariate risk embrechts2006bounds, and functionals that encompass interrelation between random variables such as mutual information or correlation. We say that a distribution $F$ dominates another distribution $G$ in the $\alpha$-th order, written as $F\succsim G$ if $D_\alpha^F(x)\leq D_\alpha^G(x)$, where \[D_1^F(x)\equiv F(x) \qquad\text{and}\qquad D_\alpha =\int_{S(x)} D_{\alpha-1}(t)d t,\] with $S(x)\coloneqq \left\{y\in \mathbb{R}^d_+: y\leq x\right\}$, see garcia2019review for details. An example of a multivariate risk measure is the multivariate lower-orthant value-at-risk at level $\alpha\in[0,1]$ of a multivariate increasing function $G$, which is defined as embrechts2006bounds \[\underline{\text{VaR}}_\alpha(G) = \partial\left\{x\in\mathbb{R}^d: G(x)\geq\alpha\right\},\] where $\partial A$ denotes the boundary of the set $A$. An analogous definition holds for the upper-orthant counterpart. The mutual information $I(X,Y)$ of two random variables $X$ and $Y$ is defined as shannon1948mathematical \[I(X,Y) = \text{KL}\left(P_{X,Y} || P_X\otimes P_Y\right),\] where $\text{KL}$ is the Kullback-Leibler divergence and $P_X\otimes P_Y$ is the independence coupling.
The above examples are only a short list of potentially interesting functionals one can analyze using the proposed method. All of these functionals rely on the information of the corresponding copula structure, which cannot be identified with a tensorized version of the changes-in-changes estimator unless in the limited setting of completely independent coordinates described previously.
In practice we only observe observations from each of the four distributions $\Prob_{Y_{C,0}}, \Prob_{Y_{C,1}}, \Prob_{Y_{T,0}}$, and $\Prob_{Y_{T,1}}$. In order to distinguish between population measures and empirical measures, we define the corresponding empirical measures of $\Prob_{Y_{C,0}}$ as $\hat{\mu}_{0,n}$, of $\Prob_{Y_{C,1}}$ as $\hat{\mu}_{1,l}$, and of $\Prob_{Y_{T,0}}$ as $\hat{\mu}_{0,m}^\star$ with corresponding sample sizes $n,l,m$. We denote the estimator of the counterfactual $\Prob_{Y_{T,1}^\dagger}$ by $\hat{\mu}_{1,m}^\dagger$. Estimating Brenier maps from data is a subject of intense activity both from theoretical and practical angles forrow2019statistical, gunsilius2021convergence, hr2021minimax, delara2021consistent, ManBalNil21,manole2021sharp,pooladian2021entropic,VacMuzRud21,muzellec2021near,deb2021rates. Further details about estimating Brenier maps from data can be found in Appendix (ref). In this section we state propositions about the consistency of the estimated counterfactual distribution and its rate of convergence in expected Wasserstein distance to the population distribution. The $2$-Wasserstein distance between measures $\mu$ and $\nu$, which we denote as $W_2\left(\mu,\nu\right)$, is the square root of the infimum transportation cost (ref) between the two measures. Proofs of these propositions are deferred to Appendix $\ref{sec:proofs}$.
The observations form empirical measures $\hat{\mu}_{0,n}$, $\hat{\mu}_{0,m}^\star$, and $\hat{\mu}_{1,l}$ with $n$, $m$, and $l$ denoting the respective sample sizes. These samples are assumed to be mutually independent. An estimator for the natural trend, $\sdhat$ is a minimizer of $W_2(\hat{\mu}_{0,n}, \hat{\mu}_{1,l})$. A plug-in estimator for the counterfactual distribution is given by $\sdhat_{\# \hat{\mu}_{0,m}^\star}$. However, $\sdhat$ is non-trivially defined only on the support of $\hat{\mu}_{0,n}$. Therefore in practice we project $\hat{\mu}_{0,m}^\star$ onto the support of $\hat{\mu}_{0,n}$ through one-nearest neighbor Euclidean matching, itself a 2-Wasserstein optimal transport problem. The “rounded” measure of treatment outcomes pre-intervention is denoted $\tilde{\mu}_{0,m}^\star$. This leads to our proposed estimator for the counterfactual distribution: ${\sdhat_{\# \tilde{\mu}_{0,m}^\star} =: \hat{\mu}_{1,m}^\dagger}$.
Under the assumptions on population measures stated in Assumption (ref), the Wasserstein distance between the estimated and population counterfactual distribution converges uniformly to zero. This enables a consistency statement about that distance converging to probability to zero, as well as a stronger statement about the Wasserstein-distance of the estimated and population natural trend maps.
Under mild further assumptions, we show that the expected 2-Wasserstein distance between the population and estimated counterfactual measures converges to $0$ at a $O(\text{min}(m,n)^{\nicefrac{1}{d}})$ rate.
Combining these assumptions with those previously stated for consistency under Assumption (ref) yields the following claim.
Obtaining the large sample distribution of optimal transport maps is largely an open problem. However, rates of convergence are known deb2021rates, hr2021minimax, which enables inference on these functionals through subsampling methods. Under mild assumptions politis2001asymptotic, romano2012uniform, shao2013fixed, subsampling theory provides confidence regions and hypothesis tests with strong asymptotic guarantees. Even if the rate of convergence for the functional is unknown, that rate may still be estimated.
In this section, we demonstrate the performance of the optimal transport-based estimator with two numerical experiments and apply it to a classical causal inference dataset.
We first present a simple simulation experiment which demonstrates that improper modeling of dependency between coordinates of multivariate outcomes can lead to highly biased estimates. In particular, our proposed method manages to estimate the correct bivariate counterfactual distributions, while the multivariate extension of the CiC estimator through tensorization does not. Consider the following set of linear production functions in which map latent vectors in $\mathbb{R}^2$ to outcomes without intervention in $\mathbb{R}^2$: \[h_0(u) =
u\qquadand\qquad h_1(u) =
u.\]
These production functions are a cyclically comonotone but not element-wise monotone as required for the CiC estimator. Choosing Beta distributions as the laws of the latent variables, we generate $n=3000$ independent outcomes from the control population before and after the intervention, as well as from the treatment group before the intervention and its unobserved counterfactual. Further details about the implementation can be found in Appendix (ref).
The recovered marginals in Figure (ref) and the kernel density estimator (KDE) of the counterfactual joint distribution in Figure (ref) demonstrate that our methodology estimates the true counterfactual distribution almost exactly while the CiC estimator cannot. Notably, the CiC estimator recovers a joint distribution with a mirrored dependence structure.
As previously noted, the tensorized CiC estimates counterfactuals under the strong assumption that the copula of outcomes pre- and post-intervention are identical. Each coordinate of the potential outcomes vectors is assumed to be a monotone increasing function of the latent random vector's corresponding coordinate. When there are interactions between latent coordinates during potential outcome generation the multivariate CiC will not identify counterfactual distributions. We believe that proper accounting for correlation structures contributes to our seemingly novel substitution result when reanalyzing the classical minimum wage data from card1994minimum in Section (ref).
In the following numerical experiment, we consider the performance of our estimator against the multivariate tensorized CiC estimator as the dimensionality of the observations grows. Optimal transport's strong recovery of higher-dimensional joint distributions, illustrated by Figure (ref), distinguishes it.
Recall from Theorem $\ref{thm:brenier}$ that optimal transport maps can be characterized as the gradients of convex functions. Consider the function $f(x,y) = x^2 y^{-1}$ which is convex when $y > 0$. As the sum of convex functions is convex, summing $f$ evaluated at multiple pairs of coordinates in a random vector constructs a convex function from $\mathbb{R}^d$ to $\mathbb{R}$. Formally we consider a class convex functions $g(x; \theta) = \sum_{i=1}^p x_{(i_1)}^2 x_{(i_2)}^{-1}$ where $x \in \mathbb{R}_{++}^d$, $x_{(i)}$ represents the $i^{\text{th}}$ coordinate of $x$, and $i_j \in [[d]]$ for $j \in \{1, 2\}$; the parameter $\theta$ contains the $p$ pairs of indices to be compared. The gradient of $g$ is straightforward to compute. In this simulation, we specify $h_0$ to be the identity and $h_1$ to be $\nabla g$. There are $d$ coordinate pairs randomly selected each run. Both treatment arms have $n=3,000$ units drawn from Beta latent distributions. Further implementation details can be found in Appendix (ref).
Our experiment is designed quantify how well the two methods estimate a joint counterfactual distribution. For computational reasons, we focus on randomly selected pair-wise coordinate interactions. To this end, we construct the empirical CDF (eCDF) of the true counterfactual, its optimal transport estimate, and its CiC estimate over a uniform mesh of $10,000$ points. We compute the mean absolute difference between the true and estimated eCDFs over each point in the mesh and have $20$ runs per trial. As demonstrated in Figure (ref), our proposed method better estimates the quantiles of a plane in the outcome space. Table (ref) shows that optimal transport has a mean absolute difference consistently smaller than changes-in-changes with a smaller standard deviation as well. Our method's monotone increasing loss likely occurs due to the curse of dimensionality and the difficulty of matching vectors in very high dimensions with only $3000$ samples. Regardless, optimal transport still has substantially smaller MAD and standard deviation than CiC when $d=100$.
On April 1, 1992, New Jersey raised its minimum wage from the Federal level of $\$ 4.25$ per hour to $ \$ 5.05$ per hour. card1994minimum (CK henceforth) investigated the effect this increase had on employment in fast food restaurants with a case-control experimental design, taking the bordering region of eastern Pennsylvania, where the minimum wage remained at $\$ 4.25$, as a control group. The original study employs the standard DiD estimator to conclude that the higher minimum wage led to increased employment in New Jersey restaurants. This result spawned much debate within the economics community about the effect of minimum wage policies on employment (e.g., neumark2000minimum, dube2010minimum, meer2016effects, neumark2014revisiting, and references therein).
The treatment effect of interest in CK and many subsequent reevaluations is the change in full-time equivalent employees (FTE), which is defined via ${\text{FTE} = \text{FT} + .5 \text{PT}}$, where FT (PT) is the number of full-time (part-time) employees neumark2000minimum,lu2004optimal. We use the proposed multivariate extension of the changes-in-changes estimator to estimate a bivariate treatment effect in terms of full- and part-time employees, hence fully dissecting the causal effect of raising the minimum wage on these subgroups.
We analyze the number of full- and part-time employees as drawn from an underlying continuous distribution. In the original Card and Krueger data, “there are 28 records that report fractional full-time employees (the fraction is always one-half), 29 records that report fractional numbers of part-time employees, and one record that reports a fractional number of managers” ehrenberg1995myth. These partial observations could represent ambiguity in the classification of a worker as full- or part-time ehrenberg1995myth. As the restaurant employee responding to the phone interview was not necessarily the same at both measurements, neumark2000minimum argue, ”there is no reason to believe that the responses in the first and second waves are based on the same `definition' of employment.” Therefore these fractional observations can also be considered a degree of belief statement. As discussed above, deterministic optimal transport maps are the unique solution in the continuous case while the discrete case also admits probabilistic optimal transport plans as solutions. In practice we use a linear program implemented in the Python package POT which finds a deterministic plan, matching the structure of our theoretical results.
The original CK study finds a positive average treatment effect of $2.76$ FTE jobs added in the treatment group compared to the control group, a result reproduced by neumark2000minimum using different methodology. neumark2000minimum also compute an average treatment effect in full- and part-time jobs with separate regressions and find treatment effects of $3.16$ gained full-time jobs and $0.60$ lost part-time jobs. These classical contributions only focus on aggregate outcomes over all restaurants, irrespective of the size of the respective restaurant. An exception is ropponen2011reconciling who reanalyzes the data with the CiC estimator, taking into account the heterogeneity in the size of the fast-food restaurants. He finds the average treatment effect bounded in $[0.90,1.70]$. Ropponen also notes in New Jersey that large restaurants pre-intervention (measured in FTE employees) have a negative treatment effect while small restaurants pre-intervention have a positive treatment effect. Our analysis recovers this trend for both FT and PT employees in New Jersey. As we also generate counterfactual controls, we find an opposite trend for Pennsylvania restaurants. This opposite trend is intriguing and worth exploring further as it does not seem to be explained by substitution effects of employees moving across the state line to find jobs.
The previously discussed subsampling approach to inference was applied to generate symmetric 95 % confidence intervals around each point estimate. The only interval to contain zero is for the DiD treatment effect on part-time employment. Furthermore, we run subsampled hypothesis tests to check whether the point estimates are significantly different across estimators. Of the four tests between optimal transport to an alternative, only the difference in treatment effect on part-time employment estimated with CiC was not significantly different from zero at the $p = 0.05$ level. For both the confidence intervals and hypothesis tests, we use a block size of 300 and 10,000 replications.
The average treatment effects estimated using the OT and CiC counterfactual distributions depend only on the marginals, not the full joint structure. However as illustrated in Figure (ref), the strong copula assumptions between pre- and post-intervention periods for CiC can lead to markedly different estimated marginals. That structure implies restaurants with the largest pre-intervention full/part-time employment will also have the largest counterfactual full/part-time employment. The wider tails of the CiC treatment effect marginals compared to OT suggest that additional joint structure in this setting leads to counterfactuals with less variance from the observed potential outcome.
The proposed multivariate extension of the changes-in-changes estimator allows us to jointly estimate the effects on full- and part-time employment while accounting for the heterogeneity in restaurant size. For the results in this section, restaurants are represented in $\mathbb{R}^2$ by our outcomes of interest, the number of full- and part-time employees. In the supplementary materials we provide an additional experiment with restaurants represented as higher dimensional vectors including additional measurements about wages, prices, and operations. We discard $19$ units with missing outcomes, leaving $76$ control and $315$ treatment restaurants. The transport plans for both groups are solved with a linear program, and we apply nearest neighbor matching between treatment and control pre-intervention samples to calculate counterfactual outcomes.
Our optimal transport analysis suggests that fast food restaurants responded to an increased minimum wage by substituting full-time employees for part-time ones with a net gain of full-time equivalent employees. Our results seem to indicate that the negative effect on part-time employees is more pronounced than previously estimated. The quantile plot in Figure (ref) suggests this effect applies throughout the distribution of restaurants, not just the mean: fixing the number of full-time employees, counterfactual quantiles tend to have more part-time employees than treatment group quantiles at the same level; likewise fixing the number of part-time employees, treatment group quantiles tend to have more full-time employees. In their original paper CK suggest that restaurants may respond to higher labor costs with more full-time employees because they tend to be older, more skilled, and more productive, thus a better investment of capital. Furthermore, the univariate methods may underestimate the number of part-time employees lost, a conclusion supported by the higher dimensional analysis in Table (ref) from the supplementary materials.
We have introduced a general method based on optimal transport theory for causal inference in observational treatment and control study designs. It combines two desirable properties: it is designed for multivariate settings while at the same time capturing the heterogeneity in treatment response of individuals. In particular, it complements both the classical difference-in-differences estimator, which is applicable in multivariate settings but only captures average effects, and the changes-in-changes estimator athey2006identification, which allows for general treatment heterogeneity but is only applicable in univariate settings.
Showcasing the utility of the proposed method, we revisit the classical Card & Krueger dataset by decomposing the treatment effects for full-and part-time employees. We find that fast-food restaurants responded to an increased minimum wage by substituting full-time employees for part-time employees, even when accounting for restaurant size. This provides a novel insight by dissecting the relationship between full- and part-time employees while accounting for the heterogeneity in restaurant size. More generally, since our proposed method is able to consider multivariate outcomes, it can help uncover interesting causal effects and nonlinear relations in a wide array of applications.
William Torous was supported as an undergraduate through MIT's UROP program with funding from the Lord Foundation and the Paul E. Gray UROP Fund. Florian Gunsilius was supported by a MITRE faculty research award. Philippe Rigollet was supported by NSF grants IIS-1838071, DMS-2022448, and CCF-2106377.
The authors state no conflict of interest.
The datasets analyzed during the current study are available at David Card's public website, \url{https://davidcard.berkeley.edu/data_sets.html}.
Implementations of our methodology in $\texttt{R}$, $\texttt{Python}$, and $\texttt{STATA}$ are under development.