The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
44,126 characters
Non-linear Triple Changes Estimator for Targeted Policies
\maketitle
\begin{abstract}
The renowned difference-in-differences (DiD) estimator relies on the assumption of `parallel trends,' which does not hold in many practical applications.
To address this issue, the econometrics literature has turned to the triple difference estimator.
Both DiD and triple difference are limited to assessing average effects exclusively.
An alternative avenue is offered by the changes-in-changes (CiC) estimator, which provides an estimate of the entire counterfactual distribution at the cost of relying on (stronger) distributional assumptions.
In this work, we extend the triple difference estimator to accommodate the CiC framework, presenting the `triple changes estimator' and its identification assumptions, thereby expanding the scope of the CiC paradigm.
Subsequently, we empirically evaluate the proposed framework and apply it to a study examining the impact of Medicaid expansion on children's preventive care.
\end{abstract}
\section{Introduction}
In the domains of econometrics and quantitative social sciences, the \emph{difference-in-differences} (DiD) estimator and its extensions have emerged as indispensable tools for estimating causal effects in observational studies.
With its roots traced back to the work of \citet{snow1855mode} and further popularized by seminal works such as \citet{ashenfelter1984using} and \citet{card1993minimum},
DiD has provided a reliable method for evaluating the effect of policy interventions over time.
In particular, DiD stands out in studies where inference based on controlling for confounders or using instrumental variables is deemed unsuitable, and where pre-treatment information is available.
Researchers have actively pursued various extensions of the DiD estimator, aimed at enhancing its applicability and robustness \cite{athey2006identification, sofer2016negative, callaway2021difference, roth2023parallel}.
The DiD estimator hinges on the assumption of parallel trends, which states that the (average) outcome of both the control and treatment groups share exactly the same evolution trend in the absence of the treatment.
Mathematically speaking, this assumption translates to
\begin{equation}\label{eq:common-trends}
\mathbb{E}[Y^0(t_1)-Y^0(t_0)\,\vert\, D=1] = \mathbb{E}[Y^0(t_1)-Y^0(t_0)\,\vert\, D=0],
\end{equation}
where $Y^0(t)$ denotes the potential outcome associated with the absence of the treatment at time $t$, and $D$ denotes the assigned treatment.
In particular, $D=1$ and $D=0$ represent the treatment and control groups, respectively.
The validity of this assumption can be challenged, and parallel trends might get violated due to unobserved time-varying factors or dynamic changes in the study context.
This violation may result in biased estimates.
As a response to potential deviations from parallel trends, researchers have turned to more flexible models, such as the \emph{triple difference} framework \cite{gruber1994incidence}.
The triple difference estimator can be formulated as the difference between two DiD estimators.
Intuitively, the difference of two DiD estimators is unbiased, provided both estimators have the same bias.
Indeed, sole purpose of of subtracting the second DiD estimator is to debias the first one.
Despite the prevalent use of the triple difference estimator, especially over the last two decades \cite{raifman2018association, sakurai2020relationship, han2016effect, chen2020triple, tai2001racial}, only recently a formal presentation of the framework and its identification assumptions was provided \cite{olden2022triple}.
\import{./figures/}{fig0.tex}
The triple difference framework guarantees only the identification of the \emph{average treatment effect} (on the treated).
Another issue that arises is that the triple difference estimator is biased if the outcomes in the pre- and post-treatment or among control and treatment groups are measured on a different scale.
The \emph{changes-in-changes} (CiC) framework \citep{athey2006identification}, on the other hand, is scale-invariant and yields the identification of the counterfactual outcome probability distribution.
To do so, CiC framework requires additional assumptions beyond Eq.~\eqref{eq:common-trends}.
It postulates that there exists a unique monotone mapping $T_d$ which maps the probability measure over $Y^0(t_0)$ to that over $Y^0(t_1)$ within group $D=d$ for $d\in\{0,1\}$ (see Figure \ref{fig:one}.)
Moreover, these two mappings are assumed to be identical.
That is, $T_1(y) = T_0(y)$ for every $y$, or since the mappings are bijective,
\begin{equation}\label{eq:tid}
T_1\circ T_0^{-1} = \textrm{Id},
\end{equation}
where $\textrm{Id}$ represents the identity map.
Eq.~\eqref{eq:tid} states that there is \emph{no drift} in the evolution trend across groups.
~Note that Eq.~\eqref{eq:tid}, was not explicitly stated in \citet{athey2006identification} as an assumption but can be derived as a consequence of the following four assumptions:
\begin{assumption}[Model assumption]\label{as:model}
The potential outcomes $Y^0(t)$ can be modelled using a production function $h$ of a latent variable $U$, which models the individual characteristics:
\begin{equation*}
\forall t\in\{0,1\}:\quad Y^0(t) = h(U; t).
\end{equation*}
\end{assumption}
In particular, \ref{as:model} posits that $h$ does not depend on the group assignment ($D$).
\begin{assumption}[Strict monotonicity]\label{as:monotone}
The function $h(\,\cdot\,; t)$ is strictly increasing in $U$ for every $t\in\{t_0,t_1\}$.
\end{assumption}
\begin{assumption}[Time invariance]\label{as:invariance}
Within every subgroup, the distribution of the latent variable $U$ does not change over time.
\end{assumption}
\begin{assumption}[Latent support overlap]\label{as:sup}
The support of the latent variable $U$ in the treatment group is a subset of its support in the control group.
\end{assumption}
In the case of a one-dimensional outcome, mappings $T_d$ can be expressed as
\begin{equation}\label{eq:td}
T_d= F^{-1}_{Y^0(t_1)\,\vert\, D=d}\circ F_{Y^0(t_0)\,\vert\, D=d}
\end{equation}
where $F_{Y^0(t)\,\vert\, D=d}(\cdot)$ represents the cumulative density function of $Y^0(t)$ in group $D=d$.
Under regularity conditions, $T_d$
is the unique monotone map that pushes forward the probability measure over $Y^0(t_0)$ to that over $Y^0(t_1)$ in group $D=d$ \citep{villani2009optimal, santambrogio2015optimal}.
Combining Equations \eqref{eq:tid} and \eqref{eq:td}, it is straightforward to identify the counterfactual distribution of $Y^0(t_1)$ in the treatment group, $F_{Y^0(t_1)\,\vert\, D=1}(y)$.
Specifically,
\[
\begin{split}
F_{Y^0(t_1)}&_{\,\vert\, D=1}(y)=
F_{Y^0(t_0)\,\vert\, D=1}\circ
F^{-1}_{Y^0(t_0)\,\vert\, D=0}\circ
F_{Y^0(t_1)\,\vert\, D=0}(y)
,
\end{split}
\]
which matches Eq.~(9) in the original work of \citet{athey2006identification}.
The no-drift assumption specified in Eq.~\eqref{eq:tid} can be challenged in practice, especially in scenarios where the treatment or exposure is directed toward a specific sub-population.
This situation arises in studies on the impact of targeted interventions, such as a criminal justice initiative ($D$), on recidivism rates,
where this intervention is exclusively administered to individuals with specific criminal histories or risk profiles. Naturally, it is expected that the time evolution of counterfactual recidivism rates will exhibit significant disparities between the control and treated groups.
Similar challenges arise when the eligibility criteria is narrow within the context of social programs such as welfare or housing assistance, which target specific demographic groups.
Recognizing the strengths and weaknesses of both the triple difference and CiC frameworks, we propose a novel estimator that combines the best of both worlds.
Our proposed `triple changes' estimator aims to overcome the limitations of DiD and the triple difference estimator by leveraging the scale invariance and strong identification results of CiC, while introducing a more flexible mapping assumption that allows us to relax \eqref{eq:tid}.
We briefly present the contributions of this work.
\begin{itemize}[leftmargin=*]
\item We formally present and analyze the triple changes estimator as an extension to CiC, and provide the necessary assumptions for its point identification first in the scalar case.
We then discuss how to generalize our results to high-dimensional outcomes by harnessing theory of optimal transport.
\item We provide several partial identification results under relaxed versions of our proposed point identifiability assumptions.
Further, we show the validity of analogous results for the classic CiC framework as a special case of our derivations.
\item We introduce a finite-sample estimator for the average treatment effect on the treated within our framework and analyze its asymptotic behaviour.
\item We conduct an empirical evaluation of our estimator on both synthetic and real datasets.
\end{itemize}
This paper is organized as follows.
Section \ref{sec:model} reviews the necessary background and the setup of the study.
In Section \ref{sec:one}, the identification of our estimand of interest is studied for a scalar outcome.
Section \ref{sec:inf} provides an estimator for the latter and studies its asymptotic properties.
In Section \ref{sec:ot}, we extend our work to high-dimensional outcomes using theory of optimal transport.
Numerical evaluations are presented in Section \ref{sec:exp}.
\subsection{Causal model}\label{sec:model}
We consider a study where we have access to data from two sources, e.g., two states of the united states, or two cities, or any two separate populations.
These two sources of data will be denoted by $S=s_0$ and $S=s_1$ throughout.
We assume that a treatment (e.g., a health-care policy) is administered in one state, without loss of generality in $S=s_1$, and not in the other.
In both states, the individuals are partitioned into two cohorts, namely, $D=d_1$ and $D=d_0$, signifying the individuals that are eligible and not eligible for receiving the treatment, respectively\footnote{In a classic controlled trial, these would correspond to the treatment and control groups, respectively.}.
The outcome is measured in two time points, namely $t_0<t_1$, where the eligible individuals in state $s_1$ receive the treatment at an infinitesimal amount of time after $t_0$\footnote{Note that we do not limit our setting to panel data.
In particular, the individuals for which the outcome is measured may differ across time points.}.
We denote by $Y^{D=d_0}(t)$ and $Y^{D=d_1}(t)$ the potential outcome variables associated with the outcome at time $t$ in the absence, and in the presence of treatment, respectively.
To improve readability, we will often use the short-hands $Y^0(t)$ and $Y^1(t)$ for $Y^{D=d_0}(t)$ and $Y^{D=d_1}(t)$, respectively.
We denote the observed outcome at time $t$ by $Y(t)$.
Throughout, we make the following standard consistency assumption \citep{rubin1980randomization}.
\begin{assumption}[Consistency]\label{as:cons}
At each time $t$, the realized outcome $Y(t)$ is determined as
\[Y(t)=\sum_d\mathbbm{1}\{D=d\}\cdot Y^{D=d}(t),\]
where $\mathbbm{1}\{\cdot\}$ denotes the indicator function.
\end{assumption}
The estimand of interest is the effect of treatment on the treated, i.e., the group corresponding to $S=s_1,D=d_1$.
As such, we target learning the probability measure over the counterfactual outcome $Y^0$ within this subgroup:
\[F(Y^0\mid S=s_1,\:D=d_1),\]
where $F$ denotes the cumulative density function.
When clear from context, we use the shorthand $F_{Y^0\mid s_1,d_1}$ instead.
Throughout, we assume that random variables are defined over a compact domain, and that densities are absolutely continuous with respect to the Lebesgue measure.
\section{One-dimensional Estimator}\label{sec:one}
\import{./figures/}{figure1}
We commence our analysis by considering cases where the outcome of interest, $Y$, is one-dimensional, i.e., a scalar.
In contrast to \citet{athey2006identification}, we posit that the outcome of an individual can be determined by a combination of the latent variable $U$ (with a common support across groups), and the group to which the individual belongs.
This adjustment relaxes the model assumption of the CiC framework, as outlined below.
\begin{customass}{1}[Model assumption]\label{as:model2}
At each state $S=s$ and treatment group $D=d$,
the potential outcomes $Y^0(t)$ can be determined through a production function $h_{s,d}(\cdot\:;t)$ of a latent variable $U$, which models the individual characteristics:
\begin{equation*}
\forall t\in\{0,1\}:\quad Y^0(t) = h_{s,d}(U; t).
\end{equation*}
\end{customass}
A key distinction between our setup and the classic CiC framework lies in relaxing \ref{as:model} to \ref{as:model2}, allowing the production functions $h_{s,d}(\cdot)$ to be \emph{group-specific}.
In applications involving targeted treatment assignments, \ref{as:model2} emerges as a more sensible assumption.
Assumptions \ref{as:monotone} and \ref{as:invariance} are adapted analogously:
\begin{customass}{2}[Strict monotonicity]\label{as:monotone2}
The production functions $h_{s,d}(\,\cdot\,; t)$ are strictly increasing in $U$ for every $t\in\{t_0,t_1\}$, and every $s,d$.
\end{customass}
\begin{customass}{3}[Time invariance]\label{as:invariance2}
Within every subgroup, the distribution of the latent variable $U$ does not change over time.
That is, $\forall u, \forall s,d$,
\[F_{U\vert S=s, D=d, T=t_1}(u) = F_{U\vert S=s, D=d, T=t_0}(u).\]
\end{customass}
Additionally, we prefer to formulate the overlap assumption in relation to the potential outcomes rather than the latent variable $U$.
This preference arises from the broader accessibility and interpretability of the support of the outcome, as opposed to that of latent characteristics.
\begin{customass}{4}[Outcome support overlap]\label{as:sup2}
The potential outcomes $Y^0(t_0)$ and $Y^0(t_1)$ are defined over domains $\mathbb{Y}_0$ and $\mathbb{Y}_1$, respectively, which are common across subgroups $S\in\{s_0,s_1\}, D\in\{d_0,d_1\}$.
Moreover, $\mathbb{Y}_0\subseteq\mathbb{Y}_1$.
\end{customass}
Assumptions \ref{as:model2} through \ref{as:sup2} establish the existence of four distinct monotone maps that push forward the density of $Y^0(t_0)$ to that of $Y^0(t_1)$ in each subgroup.
In particular, let $T_{s,d}$ denote the monotone map that pushes forward the density of $Y^0(t_0)$ in the group corresponding to $S=s, D=d$ to the density of $Y^0(t_1)$ in the same group (see Figure \ref{fig:four}.)
Specifically, $T_{s,d}$ can be expressed in terms of the $h(\cdot)$ functions as\footnote{For the purposes of this section, Eq.~\eqref{eq:td3} can be considered as the definition of maps $T_{s,d}$.}
\begin{equation}\label{eq:td3}
T_{s,d}(y) = h_{s,d}\big( h_{s,d}^{-1}(y; t_0) ; t_1\big),
\end{equation}
or equivalently, in terms of the cumulative density functions,
\begin{equation}\label{eq:td2}
T_{s,d} = F^{-1}_{Y^0(t_1)\,\vert\, s,d}\circ F_{Y^0(t_0)\,\vert\, s,d}.
\end{equation}
It is noteworthy that under \ref{as:cons}, $F_{Y^0(t)\mid s,t}=F_{Y(t)\mid s,t}$ for every subgroup except $\{s_1,d_1\}$.
Consider the map
$T^*_s=T_{s,d_1}\circ T^{-1}_{s,d_0}$.
$T^*_s$ can be interpreted as the non-linear drift between the maps from $Y^0(t_0)$ to $Y^0(t_1)$ among control and treatment groups in any state $s$ (see Figure \ref{fig:four}.)
Identification in CiC is achieved by assuming that there is no drift, i.e., $T^*_s$ is the identity map (see Eq.~\ref{eq:tid}).
We shall proceed by relaxing this assumption as follows.
\begin{assumption}[State-independent drifts]\label{as:stindep}
The drift between the mappings of potential outcomes $Y^0(t_0)$ to $Y^0(t_1)$ in the control and treatment groups is independent of the state.
Formally,
\[T^*\coloneqq T_{s_0,d_1}\circ T^{-1}_{s_0,d_0} \equiv T_{s_1,d_1}\circ T^{-1}_{s_1,d_0}.\]
\end{assumption}
Specifically, rather than assuming \emph{no drift}, we have relaxed the assumption to \emph{equal drift} across the two states.
\ref{as:stindep} can also be expressed in terms of the production functions, albeit at the cost of interpretability:
\begin{multline}\label{eq:stindep}
h_{s_0,d_1}\Big\{
h_{s_0,d_1}^{-1}\big\{
h_{s_0,d_0}\big(
h_{s_0,d_0}^{-1}(\cdot;t_1)
;t_0\big)
;t_0\big\}
;t_1\Big\}
=\\
h_{s_1,d_1}\Big\{
h_{s_1,d_1}^{-1}\big\{
h_{s_1,d_0}\big(
h_{s_1,d_0}^{-1}(\cdot;t_1)
;t_0\big)
;t_0\big\}
;t_1\Big\}.
\end{multline}
We are now ready to state our identification result.
The complete set of proofs for our results can be found in Appendix \ref{apx:proofs}.
\begin{restatable}{theorem}{thmid}\label{thm:id}
Under assumptions \ref{as:model2} - \ref{as:sup2} and \ref{as:cons} - \ref{as:stindep}, the cumulative density function of the missing counterfactual $Y^0(t_1)$ in the group $S=s_1,D=d_1$ is identified as:
\begin{equation}\label{eq:thm1}\begin{split}
F_{Y^0(t_1)}&_{\mid s_1,d_1}(y)
=\\&
F_{Y(t_0)\mid s_1,d_1}
\circ
T^{-1}_{s_0,d_1}
\circ
T_{s_0,d_0}
\circ
T^{-1}_{s_1,d_0}
(
y)
,
\end{split}
\end{equation}
where $T_{s,d}$ is given by Eq.~\eqref{eq:td2}.
\end{restatable}
\begin{remark}
To avoid unnecessarily heavy notation, we did not discuss the observed covariates.
However, an identical analysis can be done after adjusting for the observed covariates, $X$.
In particular, production functions may depend on the observed covariates, as long as their monotonicity in $U$ is maintained for every $x$ in the domain of $X$.
\ref{as:invariance2}, \ref{as:sup2} and \ref{as:stindep} need to be valid conditioned on $X$ in this scenario, and Eq.~\eqref{eq:thm1} must hold when every term is conditioned on $X$.
\end{remark}
\begin{remark}
We articulated our identifiability assumptions in accordance with the original work of \citet{athey2006identification}.
An alternative way of presenting the assumptions would be to do it akin to the \emph{quantile-quantile equi-confounding bias} assumption proposed by \citet{ghassami2022combining}.
In our context, this would translate to directly assuming the existence of monotone maps $T_{s,d}$ based on Eq.~\eqref{eq:td2} instead of drawing conclusions from \ref{as:model2}-\ref{as:sup2} to establish it.
\end{remark}
\begin{remark}
We formulated production functions to model the potential outcomes under no treatment ($Y^0$), whereas we left the other potential outcome, $Y^1$, unrestricted.
Due to symmetry, one could model $Y^1$ using production functions and leave $Y^0$ unrestricted.
This scenario might arise for instance, in an study where 3 out of 4 cohorts receive treatment.
One should exercise greater caution in such cases however, since in practice, the treatment may have effects that significantly alter the composition of the population under study.
Under these circumstances, assuming monotone production functions for $Y^1$ could be a more drastic assumption.
\end{remark}
\subsection{Relaxing monotonicity}
As mentioned earlier, \ref{as:monotone2} (and its counterpart, \ref{as:monotone} in \citet{athey2006identification}) is an untestable assumption, and might be drastic to impose in certain applications.
The monotonicity of functions $h_{s,d}(\cdot;t)$ in $U$ has two implications:
(i) these functions are bijective, establishing a well-defined inverse for them;
(ii) the property that $\mathbbm{P}(h_{s,d}(U;t)\leq y)=\mathbbm{P}(U\leq h^{-1}_{s,d}(y;t))$, which is repeatedly utilized in proving the point identification result in Theorem \ref{thm:id} (see Appendix \ref{apx:proofs}).
While the bijectivity of $h_{s,d}(\cdot;t)$ appears to be essential for our framework to work, we can relax (ii).
~Let us first rephrase the monotonicity assumption in the equivalent form:
\begin{equation}\label{eq:rephmono}
\langle u_0- u_1, h_{s,d}(u_0;t)- h_{s,d}(u_1;t)\rangle\geq0,\quad \forall u_0,u_1.
\end{equation}
This assumption can be relaxed as follows.
\begin{assumption}[$\epsilon$-monotonicity]\label{as:asm}
Functions $h^{-1}_{s,d}(\cdot;t)$ are well-defined. Additionally, for any $u$ in the support of $U$,
\begin{equation}\label{eq:asmonotone}\mathbbm{P}\big(\langle U- u, h_{s,d}(U;t)- h_{s,d}(u;t)\rangle<0\mid s,d\big)\leq\frac{\epsilon}{2}.\end{equation}
\end{assumption}
For example, the monthly income of an individual in terms of her age after adjusting for the other covariates can be a $\epsilon$-monotone function.
In general, monthly income increases due to promotions and inflation.
However, temporary unemployment and retirement can affect this trend.
See Figure \ref{fig:eps} for a visualization.
\begin{figure}
\centering
\includegraphics[width=0.6\textwidth]{figures/agerevenue.pdf}
\caption{Gross monthly income versus age after adjusting for other covariates. }
\label{fig:eps}
\end{figure}
We expect that under this relaxation, point identification cannot be achieved.
However, the following partial identification result holds.
\begin{restatable}{proposition}{prppartial}\label{prp:partial}
Under assumptions \ref{as:model2}, \ref{as:invariance2}, \ref{as:sup2} and \ref{as:cons} - \ref{as:asm}, for any $y\in\mathbb{Y}_1$,
\begin{equation*}
\begin{split}
F_{Y(t_0)\mid s_1,d_1}&\circ
\mathrm{\underline{\Phi}}^{-1}_{s_0,d_1}
\circ
\mathrm{\underline{\Phi}}_{s_0,d_0}
\circ
\mathrm{\underline{\Phi}}^{-1}_{s_1,d_0}
(
y)
-\epsilon
\\&\hspace{3cm}\leq
F_{Y^0(t_1)\mid s_1,d_1}(
y)
\leq\\
&\hspace{4.5cm}F_{Y(t_0)\mid s_1,d_1}\circ
\mathrm{\overline{\Phi}}^{-1}_{s_0,d_1}
\circ
\mathrm{\overline{\Phi}}_{s_0,d_0}
\circ
\mathrm{\overline{\Phi}}^{-1}_{s_1,d_0}
(
y)
+\epsilon,
\end{split}
\end{equation*}
where
\[
\begin{split}
\mathrm{\underline{\Phi}}_{s,d}(y)=F^{-1}_{Y(t_1)\vert s,d}\big(F_{Y(t_0)\vert s,d}(y)-\epsilon\big),\\
\mathrm{\underline{\Phi}}^{-1}_{s,d}(y)=F^{-1}_{Y(t_0)\vert s,d}\big(F_{Y(t_1)\vert s,d}(y)-\epsilon\big),\\
\mathrm{\overline{\Phi}}_{s,d}(y)=F^{-1}_{Y(t_1)\vert s,d}\big(F_{Y(t_0)\vert s,d}(y)+\epsilon\big),\\
\mathrm{\overline{\Phi}}^{-1}_{s,d}(y)=F^{-1}_{Y(t_0)\vert s,d}\big(F_{Y(t_1)\vert s,d}(y)+\epsilon\big).
\end{split}
\]
\end{restatable}
\begin{remark}
Note that strict monotonicity is a special case of \ref{as:asm}, which corresponds to $\epsilon=0$.
Accordingly, Proposition \ref{prp:partial} reduces to Theorem \ref{thm:id} with the choice of $\epsilon=0$.
\end{remark}
\begin{remark}
An analogous result applies to the original framework of CiC when \ref{as:monotone} is relaxed to $\epsilon$-monotonicity.
See Appendix \ref{apx:CiC} for details.
\end{remark}
\subsection{Relaxing time invariance}
In certain applications, \ref{as:invariance2} may also be violated.
In particular, in repeated cross-sections, it might be challenging to maintain the same distribution among the cases under study.
Even if possible, this may both reduce the number of available samples and result in selection bias, making inference more complicated.
As such, we consider a relaxation of \ref{as:invariance2} and derive a partial identification result in this setting.
We assume that the distribution of the latent variable may change over time, but this change is bounded in Kolmogorov (aka KS) distance \cite{kolmogorov1933sulla}.
\begin{restatable}[$\delta$-invariance]{assumption}{asdelta}\label{as:delta}
Within every subgroup, the Kolmogorov distance of the distribution of the latent variable across $T=t_0$ and $T=t_1$ is bounded by $\delta$.
That is,
\[\sup_u\big\vert F_{U\vert S=s,D=d,T=t_1}(u)-F_{U\vert S=s,D=d,T=t_0}(u)\big\vert\leq\delta.\]
\end{restatable}
Note again that when $\delta=0$, \ref{as:delta} reduces to \ref{as:invariance2}.
\begin{restatable}{proposition}{prpdelta}\label{prp:delta}
Under assumptions \ref{as:model2}, \ref{as:monotone2}, \ref{as:delta}, \ref{as:sup2}, and \ref{as:cons}-\ref{as:stindep}, for any $y\in\mathbb{Y}_1$,
\begin{equation*}\small
\begin{split}
F_{Y(t_0)\mid s_1,d_1}&\circ
\mathrm{\underline{\Psi}}^{-1}_{s_0,d_1}
\circ
\mathrm{\underline{\Psi}}_{s_0,d_0}
\circ
\mathrm{\underline{\Psi}}^{-1}_{s_1,d_0}
(
y)
-\delta
\\&\leq
F_{Y^0(t_1)\mid s_1,d_1}(
y)
\leq\\
&F_{Y(t_0)\mid s_1,d_1}\circ
\mathrm{\overline{\Psi}}^{-1}_{s_0,d_1}
\circ
\mathrm{\overline{\Psi}}_{s_0,d_0}
\circ
\mathrm{\overline{\Psi}}^{-1}_{s_1,d_0}
(
y)
+\delta,
\end{split}
\end{equation*}
where
\[
\begin{split}
\mathrm{\underline{\Psi}}_{s,d}(y)=F^{-1}_{Y(t_1)\vert s,d}\big(F_{Y(t_0)\vert s,d}(y)-\delta\big),\\
\mathrm{\underline{\Psi}}^{-1}_{s,d}(y)=F^{-1}_{Y(t_0)\vert s,d}\big(F_{Y(t_1)\vert s,d}(y)-\delta\big),\\
\mathrm{\overline{\Psi}}_{s,d}(y)=F^{-1}_{Y(t_1)\vert s,d}\big(F_{Y(t_0)\vert s,d}(y)+\delta\big),\\
\mathrm{\overline{\Psi}}^{-1}_{s,d}(y)=F^{-1}_{Y(t_0)\vert s,d}\big(F_{Y(t_1)\vert s,d}(y)+\delta\big).
\end{split}
\]
\end{restatable}
\begin{remark}
Even more generally, we can simultaneously relax \ref{as:monotone2} to \ref{as:asm} and \ref{as:invariance2} to \ref{as:delta}.
See Proposition \ref{prp:partialgen} in Appendix \ref{apx:partial} for the partial identification result for this case.
\end{remark}
\subsection{Marginal contrasts and joint counterfactuals}\label{sec:func}
Theorem \ref{thm:id} guarantees the identification of the probability density of the missing counterfactual, $Y^0(t_1)$ in the group corresponding to $S=s_1,D=d_1$.
Having access to this density, we can compute any \emph{marginal contrast estimand} \cite{franks2019flexible}.
A marginal contrast estimand is an estimand that can be expressed as a functional of the marginal distribution of the counterfactual outcomes.
This includes a vast majority of commonly used estimands, such as average treatment effects, conditional average treatment effects, quantile treatment effects, risk ratios, etc.
However, in certain applications, more information is desired.
Examples of estimands that are not marginal contrasts include the \emph{distribution} of the treatment effect, individual-level treatment effects, or the quantiles of the treatment effect.
As a concrete example, consider $Z \coloneqq Y^1(t_1) - Y^0(t_1)$.
The density of $Z$ is not identifiable from merely the marginal densities of $Y^1(t_1)$ and $Y^0(t_0)$.
In particular, the treatment may have non-zero effects on a fraction of the population even if the two potential outcomes have identical distributions.
Estimands that are not marginal contrasts require stronger identifiability assumptions in general.
In this section, we discuss a stronger version of \ref{as:invariance2} that can lead to the identification of such estimands.
\begin{customassu}{3}[Strong time invariance]\label{as:invariance3}
The latent variable $U$ does not change over time, and there is no loss to follow-up.
\end{customassu}
\ref{as:invariance3} is stronger than \ref{as:invariance2}, in the sense that it completely rules out the possibility of any changes in the latent variable itself, or the cohort of the study.
In contrast, \ref{as:invariance2} would allow for loss to follow-up, as long as similar individuals were recruited for the study.
\ref{as:invariance2} even accommodates a repeated cross-sections study.
\ref{as:invariance2} also allows for the evolution of the latent variable, as long as its distribution among the study population remains the same.
On the other hand, the stronger assumption \ref{as:invariance3} allows for the identification of a wider range of causal estimands.
In particular, under \ref{as:invariance3}, the joint density of the counterfactuals $\big(Y^0(t_1), Y^{1}(t_1)\big)$ is identified.
\begin{restatable}{proposition}{prpjoint}\label{prp:joint}
Under assumptions \ref{as:model2}, \ref{as:monotone2}, \ref{as:invariance3}, \ref{as:sup2}, and \ref{as:cons} - \ref{as:stindep}, the joint density of $Y^0(t_1)$ and $Y^1(t_1)$ in the group corresponding to $S=s_1, D=d_1$ is identified as
\begin{equation*}
\begin{split}
&F_{Y^0(t_1), Y^1(t_1)\mid s_1,d_1}(y^0,y^1)
=\\&
F_{Y(t_0), Y(t_1)\mid s_1,d_1}\big(F^{-1}_{Y(t_0)\mid s_1,d_1} \circ F_{Y^0(t_1)\mid s_1,d_1}(y^0),y^1\big),
\end{split}\end{equation*}
where $F_{Y^0(t_1)\mid s_1,d_1}(\cdot)$ is given by Equation \eqref{eq:thm1}.
\end{restatable}
\section{Inference}\label{sec:inf}
In this section, we discuss the estimation aspect of our framework with a focus on the average effect of treatment on the treated.
More formally, we consider the estimation of
\[\tau\coloneqq\mathbb{E}[Y^1(t_1)-Y^0(t_1)\vert S=s_1,D=d_1].\]
Note that under the identifiability assumptions of Section \ref{sec:one}, $\tau$ is identified as
\begin{multline}\label{eq:tau}
\tau = \mathbb{E}[Y(t_1)\vert s_1,d_1] -
\mathbb{E}[F^{-1}_{Y(t_1)\mid s_0,d_1}
\circ\\
F_{Y(t_0)\mid s_0,d_1}
\circ
F^{-1}_{Y(t_0)\mid s_0,d_0}
\circ
F_{Y(t_1)\mid s_0,d_0}
\circ\\
F^{-1}_{Y(t_1)\mid s_1,d_0}
\circ
F_{Y(t_0)\mid s_1,d_0}\big(
Y(t_0)
\big)\vert s_1,d_1].
\end{multline}
We make the following assumption on the data generating mechanism to render estimation feasible.
\begin{assumption}\label{as:estimation}
Conditioned on $S=s,D=d,T=t$, variables $Y(t)$ are continuous random variables defined on a shared bounded domain $[\underline{y},\overline{y}]$, with continuously differentiable density functions $f_{sdt}$, where $f_{sdt}$ is bounded from above and away from $0$, and $\partial f_{sdt}/\partial y$ is bounded.
For all $s,d,t$, $p_{sdt}=\mathbb{P}(S\!=s,D\!=d,T\!=t)>0$, and given $s,d,t$, the samples $Y_i(t)$ are independent draws from $f_{sdt}$.
\end{assumption}
Let $N$ be the total number of observed outcome samples.
Akin to \citet{athey2006identification}, we build a finite-sample estimator for $\tau$ based on empirical estimators of cumulative density functions.
In particular, let $\{Y_{sd,i}(t)\}_{i=1}^{N_{sdt}}$ denote the independent samples of the outcome at time $t$ in group $S=s, D=d$.
We define
\begin{equation}\label{eq:hatf}
\hat{F}_{Y(t)\vert s,d}(y) \coloneqq N_{sdt}^{-1}\sum_{i=1}^{N_{sdt}}\mathbbm{1}\{Y_{sd,i}(t)\leq y\},\quad\text{and,}
\end{equation}
\begin{equation}\label{eq:hatfinv}
\hat{F}^{-1}_{Y(t)\vert s,d}(u) \coloneqq \inf\{y\in\mathbb{Y}_t, \hat{F}_{Y(t)\vert s,d}(y)\geq u\}.
\end{equation}
Finally, the estimator for $\tau$ is built as:
\begin{multline}\label{eq:hattau}
\hat{\tau}\coloneqq
N^{-1}_{s_1d_1t_1}\sum_{i=1}^{N_{s_1d_1t_1}}Y_{s_1d_1,i}(t_1) - \\N^{-1}_{s_1d_1t_0}\sum_{i=1}^{N_{s_1d_1t_0}}
\hat{F}^{-1}_{Y(t_1)\mid s_0,d_1}
\circ
\hat{F}_{Y(t_0)\mid s_0,d_1}
\circ\\
\hat{F}^{-1}_{Y(t_0)\mid s_0,d_0}
\circ
\hat{F}_{Y(t_1)\mid s_0,d_0}
\circ
\hat{F}^{-1}_{Y(t_1)\mid s_1,d_0}
\circ\\
\hat{F}_{Y(t_0)\mid s_1,d_0}\big(
Y_{s_1d_1,i}(t_0)
\big).
\end{multline}
The following theorem establishes the consistency and asymptotic normality of $\hat{\tau}$.
\begin{restatable}{theorem}{thmconsistency}\label{thm:consistency}
Under \ref{as:estimation}, $\hat{\tau}-\tau = \mathcal{O}_p(N^{-\frac{1}{2}})$, and \begin{multline}
\sqrt{N}(\hat{\tau}-\tau)\overset{D}{\rightarrow}\mathcal{N}(0,
\dfrac{V_0}{p_{s_1d_1t_1}}+
\dfrac{V_1}{p_{s_0d_1t_1}}+
\dfrac{V_2}{p_{s_0d_1t_0}}+\\
\dfrac{V_3}{p_{s_0d_0t_0}}+
\dfrac{V_4}{p_{s_0d_0t_1}}+
\dfrac{V_5}{p_{s_1d_0t_1}}+
\dfrac{V_6}{p_{s_1d_0t_0}}+
\dfrac{V_7}{p_{s_1d_1t_0}}
),
\end{multline} where $\{V_i\}_{i=0}^7$ are given by Eq.~\eqref{eq:varianceterms}.
\end{restatable}
The expression for the variance and the discussion on its estimation are postponed to Appendix \ref{apx:asymptotic}.
\section{Optimal Transport Representation and High-dimensional Extension} \label{sec:ot}
So far, we focused on the special case of scalar outcome for ease of presentation.
In this section, we discuss the generalization of our results to cover multi-dimensional outcomes.
To this end, we first review the one-dimensional case from an optimal transport point of view.
Let $\eta_{s,d}$ and $\mu_{s,d}$ denote the probability measures over $Y^0(t_0)$ and $Y^0(t_1)$ in the group corresponding to $S=s,D=d$, respectively.
Brenier's theorem \cite{brenier1991polar} implies that there exists a unique mapping $T$ such that $T_\#\eta_{s,d}=\mu_{s,d}$, i.e., $T$ pushes forward $\eta_{s,d}$ to $\mu_{s,d}$, and $T$ is the gradient of a convex function.
Moreover, $T$ is the optimal transport map with quadratic cost.
More precisely, let $\Gamma$ denote the space of joint distributions over $\big(Y^0(t_0),Y^0(t_1)\big)$ conditioned on $S=s, D=d$, that agree with the marginal densities $\eta_{s,d}$ and $\mu_{s,d}$.
The optimization problem
\begin{equation}
\inf_{\gamma\in\Gamma}\int\vert\vert y_0-y_1\vert\vert^2d\gamma(y_0,y_1)
\end{equation}
has a unique solution $\gamma^*$, where
$\big(Y^0(t_0),Y^0(t_1)\big)\sim \gamma^*$ if and only if $Y^0(t_0)\sim \eta_{s,d}$ and $Y^0(t_1) = T\big(Y^0(t_0)\big)$, $\eta_{s,d}-a.s$.\footnote{Note that we have omitted the implicit conditioning on $S=s,D=d$ in our notation to improve readability.}
It is straightforward to verify that a one-dimensional function is the gradient of a convex function if and only if it is monotone.
Therefore, the monotonicity of $h_{s,d}(\cdot)$ and the monotonicity of $T_{s,d}(\cdot)$ (as a consequence of the latter) imply $T_{s,d}\equiv T$.
In other words, the mapping $T_{s,d}$ is identified as the optimal transport map with quadratic cost that pushes forward $\eta_{s,d}$ to $\mu_{s,d}$.
Indeed, monotonicity can be slightly relaxed.
The identifiability of the map $T_{s,d}$ in one dimension is guaranteed under \emph{co-monotonicity} of $h_{s,d}(\cdot;t_0)$ and $h_{s,d}(\cdot;t_1)$:
\begin{equation}\label{eq:comonotone}
\langle u_0-u_1, h_{s,d}(u_0;t_0)-h_{s,d}(u_1;t_1)\rangle\geq 0, \:\: \forall u_0,u_1.
\end{equation}
In order to achieve identifiability results in higher dimensions based on Brenier's theorem, we need to make sure that $T_{s,d}$ is the gradient of a convex function.
It is known that a function is the gradient of a convex function if and only if it is \emph{cyclically monotone} \cite{rockafellar1970convex}.
A function $h$ is said to be cyclically monotone if for any sequence $x_0,\dots,x_n$ in its domain,
\begin{equation}\label{eq:cyclic}
\sum_{i=0}^n\langle h(x_i), x_i-x_{i+1}\rangle\geq0,\end{equation}
where $x_{n+1}=x_0$.
For $n=1$, Eq.~\eqref{eq:cyclic} reduces to Eq.~\eqref{eq:rephmono}.
With this preliminary discussion in place, the identification result in higher dimensions can be stated as follows.
\begin{customassu}{2}[Co-cyclic monotonocity]\label{as:cocyclic}
Functions $h_{s,d}(\cdot;t_0)$ and $h_{s,d}(\cdot;t_1)$ are co-cyclically monotone.
That is, for any sequence $u_0,\dots, u_n$ in their common domain,
\[\sum_{i=0}^n\langle h_{s,d}(u_i;t_0), h_{s,d}(u_i;t_1)-h_{s,d}(u_{i+1};t_1)\rangle\geq0,\]
where $u_{n+1}=u_0$.
\end{customassu}
\begin{restatable}{theorem}{thmhighd}\label{thm:highd}
Under assumptions \ref{as:model2}, \ref{as:cocyclic}, \ref{as:invariance2}, \ref{as:sup2}, and \ref{as:cons} - \ref{as:stindep}, the probability measure over the missing counterfactual, i.e., $\mu_{s_1,d_1}$, is identified as
\begin{equation}
\mu_{s_1,d_1} = (T^*\circ T_{s_1,d_0})_\#\eta_{s_1,d_1},
\end{equation}
where $T_{s,d}$ is the Brenier map that pushes forward $\eta_{s,d}$ to $\mu_{s, d}$ for every $s\in\{s_0,s_1\}, d\in\{d_0,d_1\}$, and $T^*$ is the Brenier map that pushes forward $(T_{s_0,d_0\#}\eta_{s_0,d_1})$ to $\mu_{s_0,d_1}$.
\end{restatable}
Theorem \ref{thm:highd} is in essence a generalization of Theorem \ref{thm:id}.
To see this, note that the Brenier map that pushes forward $\eta_{s,d}$ to $\mu_{s,d}$ in one dimension is precisely $F^{-1}_{Y^0(t_1)\vert s,d}\circ F_{Y^0(t_0)\vert s,d}$.
However, the co-cyclic monotonicity assumption \ref{as:cocyclic} is not as easy to interpret as \ref{as:monotone2}.
To address this challenge, we follow the proposition of \citet{torous2021optimal} based on the following result.
\begin{proposition}[\citealp{saks2005weak}]\label{prp:saks}
Let $K$ and $F$ be a convex and a finite subset of $\mathbb{R}^d$ respectively. A function $T:K\to F$ is cyclically monotone if and only if it is monotone.
\end{proposition}
Proposition \ref{prp:saks} implies that in our setting, if the densities $\eta_{s,d}$ are supported on convex sets and $Y^0(t_1)$ is finite-valued in every subgroup, then \ref{as:cocyclic} reduces to co-monotonicity \eqref{eq:comonotone}.
In other words, restricting $\eta_{s,d}$ and $\mu_{s,d}$ densities to be defined over convex and finite sets, respectively, Assumption \ref{as:cocyclic} of Theorem \ref{thm:highd} can be replaced by the easily interpretable assumption of Eq.~\eqref{eq:comonotone}.
\begin{figure*}[t]
\centering
\begin{subfigure}[b]{0.375\linewidth}
\includegraphics[width=\textwidth]{figures/linear.pdf}
\caption{Linear production functions.}
\label{fig:syn-linear}
\end{subfigure}\hspace{.4cm}
\begin{subfigure}[b]{0.405\linewidth}
\includegraphics[width=\textwidth]{figures/nonlinear.pdf}
\caption{Nonlinear production functions.}
\label{fig:syn-nonlinear}
\end{subfigure}
\caption{Relative bias of the estimators evaluated on (\subref{fig:syn-linear}) a linear, and (\subref{fig:syn-nonlinear}) a non-linear model.}
\label{fig:syn}
\end{figure*}
\section{Simulation Studies}\label{sec:exp}
Our empirical analysis is structured into two main sections\footnote{
The code to reproduce the results of this paper are accessible at https://github.com/SinaAkbarii/Triple-Changes.}.
In the first part, we compare the estimation error of the triple changes estimator and DiD, triple difference, and CiC estimators using synthetically generated datasets, where we know the ground truth: i.e., the treatment effect on the treated.
To evaluate the performance of each estimator, we focus on the relative bias metric, defined as
$\varepsilon=\big\vert 1-\frac{\widehat{\tau}}{\tau}\big\vert,$
where
$\tau$ is the average treatment effect on the treated, given by \eqref{eq:tau}, and $\widehat{\tau}$ represents the estimate provided by the respective estimator.
In the second part, we apply our estimator to data from \emph{National Survey on Children's Health} (NSCH)\footnote{https://mchb.hrsa.gov/national-survey-childrens-health-questionnaires-datasets-supporting-documents} to assess the effect of Medicaid expansion under Affordable Care Act (ACA) in the united states on preventive care for children.
\subsection{Synthetic data}
We begin with a simple linear model where the latent variable $U$ given $s,d$ is sampled from a Gaussian distribution with mean $\nu_{s,d}$ and variance $1$ (see Appendix \ref{apx:exp} for a complete table of parameters $\nu_{s,d}$, as well as further details of our experiment setup,) and the production functions are defined as
\begin{equation}\label{eq:prod}
h_{s,d}(u;t)\coloneqq 2u+\big(\frac{1+s}{4}+\frac{d-0.5}{2}\big) t.
\end{equation}
The actual outcome at time $t_1$ in the treated group $(s_1,d_1)$ is also sampled from a Gaussian distribution with mean $2.75$ and variance $1$.
CiC and triple changes estimators use estimates of the cumulative density functions and their inverses.
We implemented two versions of each of these estimators, namely a non-parametric estimator using the empirical estimator of density functions (see Eq.~\ref{eq:hatf} and \ref{eq:hatfinv},) and a model-based version using the maximum likelihood estimator given a class of distributions.
The class of distributions in this section are specified correctly, i.e., the true model lies in the class.
See appendix \ref{apx:exp} for a setting where this class is misspecified.
In Figure \ref{fig:syn}, the triple difference and triple changes estimators are denoted by DDD and CCC, respectively.
As depicted by Figure \ref{fig:syn-linear}, DiD and CiC estimators exhibit persistent bias, not converging to zero even with increasing sample size.
In contrast, the biases of triple difference and triple changes estimators approach zero with as the sample size grows.
Furthermore, since the triple difference estimator only relies on empirical averages, it shows slightly lower bias in small sample sizes.
To add non-linearities to the previous model, we modified the data generating mechanism, including two of the production functions.
Specifically, for $(s_0,d_1,t_1)$, and $(s_1,d_1,t_1)$, the production functions were modified to
\[h_{s,d}(u;t) = 0.1\exp\Big(2u+\big(\frac{1+s}{4}+\frac{d-0.5}{2}\big) t\Big),\]
whereas the other six groups were generated according to Eq.~\eqref{eq:prod}.
The outcomes for this model are illustrated in Figure \ref{fig:syn-nonlinear}. Notably, under the nonlinear model, the triple-difference estimator yields biased estimates, while the triple-changes estimators remain (asymptotically) unbiased. Furthermore, owing to the complexity of the model, empirical estimators of density functions exhibit more bias compared to their MLE-based counterparts, particularly in low-sample regimes.
\subsection{Application to NSCH data}
The Medicaid expansion under Affordable Care Act was designed to extend Medicaid coverage to more low-income citizens in the US.
This expansion raised the income threshold for Medicaid eligibility to 138\% of the federal poverty level (FPL).
While Louisiana ($s_1$) adopted this expansion in 2016, the neighboring states of Mississippi and Texas ($s_0$) are yet to do so.
We utilized publicly available anonymous NSCH data for the years 2016 ($t_0$) and 2017 ($t_1$) to assess the impact of Medicaid expansion on children's access to preventive healthcare, specifically analyzing the change in the frequency of doctor visits. Individuals with FPL $\leq100\%$ were classified as eligible for the expansion ($d_1$), while those with FPL $\geq140\%$ were considered not eligible ($d_0$). Data within the uncertainty margin between these thresholds was excluded.
Applying the triple difference and triple changes estimators on the frequency of doctor visits for preventive purposes with 1000 bootstraps resulted in means of $0.170$ and $0.145$, respectively, with $90\%$ confidence intervals of $[-0.004, 0.344]$ and $[-0.01, 0.331]$, respectively.
Our findings suggest that the adoption of Medicaid expansion in Louisiana increased the likelihood of children undergoing an annual preventive care visit, compared to those in non-expansion states. This aligns with the conclusion of \citet{roy2020impact}.
\section{Concluding Remarks}
We formally presented the triple changes estimator for assessing the treatment effect on the treated in observational studies.
We derived a set of necessary assumptions for point identification and discussed partial identification results under relaxed versions of these assumptions.
Future research avenues could include exploring partial identification under alternative assumptions, investigating extensions to time series settings, and conducting statistical analyses to enhance the robustness and applicability of the estimator.
\bibliographystyle{plainnat}
\bibliography{biblio}
\clearpage