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.
89,576 characters · 11 sections · 37 citation commands
Distributional Treatment Effect with Latent Rank Invariance
Keywords: distributional treatment effect, proximal inference, finite mixture, \\ nonnegative matrix factorization, Neyman orthogonality.
JEL classification codes: C13
The fundamental limitation that we cannot simultaneously observe the two potential outcomes\textemdash treated potential outcome and untreated potential outcome\textemdash for a given unit makes the task of identifying the distribution of treatment effect particularly complicated. Thus, instead of estimating the entire distribution of treatment effect, researchers often estimate some summary measures of the treatment effect distribution, such as the average treatment effect (ATE) or the quantile treatment effect (QTE). These summary measures provide insights into the treatment effect distribution and thus help researchers with policy recommendations. However, there still remain a lot of questions that can only be answered with the distribution of the treatment effect: e.g., is the treatment Pareto improving?; how heterogeneous is the treatment effect at the unit level? Also, the distribution is important in empirical contexts where participation cannot be mandated. To anticipate the participation rate, we need to identify the share of people who are better off under the treatment regime and thus would select into treatment.
Consider a potential outcome setup with a binary treatment: $$ Y= D \cdot Y(1) + (1-D) \cdot Y(0). $$ $Y(1)$ is the treated potential outcome, $Y(0)$ is the untreated potential outcome, and $D \in \{0,1\}$ is the binary treatment variable. The questions above correspond to testing $H_0: F_{Y(1)-Y(0)}(0)=0$ and estimating $\mathrm{Var} \big( Y(1)-Y(0) \big)$. Note that these quantities, $F_{Y(1)-Y(0)}(0)$ and $\mathrm{Var} \big( Y(1) - Y(0) \big)$, all come from the distribution of individual-level treatment effect $Y(1) - Y(0)$. To answer questions that relate to the distributional concerns in policy recommendation more broadly, I focus on the following two parameters of interest:
The first parameter is the joint distribution of the two potential outcomes and the second parameter is the marginal distribution of the treatment effect. For the rest of the paper, I refer to these quantities as the distributional treatment effect (DTE) parameters.\footnote{Some previous works in the literature use the terminology `distributional effect' to discuss parameters that are a functional of the marginal distributions of the potential outcomes; e.g., FP2016. To avoid confusion, I will reserve the expression `distributional' to only when the object involves the joint distribution of the two potential outcomes.}
When we believe that there is no dependence between the two potential outcomes, meaning that a realized value of the treated potential outcome has no information on the individual-level heterogeneity and thus has no predictive power for the untreated potential outcome and vice versa, identification of the joint distribution of the two potential outcomes becomes trivial. Once we identify the marginal distributions of the two potential outcomes, the joint distribution becomes their product. However, this assumption is extremely restrictive. Thus, I instead assume conditional independence, by assuming a scalar latent variable that captures the individual-level heterogeneity in terms of the dependence between the two potential outcomes. For illustration, consider a simple additive model: the two potential outcomes are constructed with a unit-level latent variable $U \in \mathcal{U}\subset \mathbb{R}$ and two treatment-status-specific random shocks $\varepsilon(1)$ and $\varepsilon(0)$:
When
we can characterize the joint distribution of the two potential outcome as follows: $$ \Pr \left\lbrace Y(1) \leq y_1, Y(0) \leq y_0 \right\rbrace = \mathbf{E} \left[ \Pr \left\lbrace Y(1) \leq y_1 | U \right\rbrace \cdot \Pr \left\lbrace Y(0) \leq y_0 |U \right\rbrace \right]. $$ Thus, the task of identifying the joint distribution of the two potential outcomes becomes that of identifying the condtional distribution of $\varepsilon(1)$ given $U$, the conditional distribution of $\varepsilon(0)$ given $U$, and the marginal distribution of $U$.
To identify the conditional distribution of $\varepsilon(d)$ given $U$ and the marginal distribution of $U$, I assume that there are two additional proxy variables $X,Z$ that are conditionally independent of each other and the potential outcomes, given $U$. This identififcation strategy is drawn from the nonclassical measurement error literature and the proximal inference literature: see HS2008,MGT,D2023,K2023,N2022 and more. In the simple example (ref)-(ref), the proxy variables $X,Z$ will shift $\mu(U)$ independently of $\big( \varepsilon(1), \varepsilon(0) \big)$, allowing us to decompose the variation of $Y(d)$ into the variation of $U$ and the variation of $\varepsilon(d)$.
Additionally, since I do not adopt the `measurement error' interpretation on the proxy variables as in the nonclassical measurement error literature, I assume that there exists a functional of the conditional distribution of the potential outcomes given the latent variable, which strictly increases in the latent variable $U$. An example of such a functional is conditional expectation. Suppose that the two conditional expectations $\mathbf{E} \left[Y(1)|U=u\right]$ and $\mathbf{E}\left[Y(0)|U=u\right]$ are strictly increasing in $u$. In this example, the latent variable $U$ can be thought of as the rank of the conditional expectations $\mathbf{E}[Y(1)|U]$ and $\mathbf{E}[Y(0)|U]$; hence `latent rank invariance.' The conditional independence assumption and the latent rank invariance assumption are the key assumptions in identification.
In developing estimators for the distributional treatment effect parameters, I additionally assume a finite support on $U$. The finite support assumptions has several merits. Firstly, it motivates a simple estimation method based on a nonnegative matrix factorization algorithm. Secondly, the conditional independence assumption can be interpreted as finite mixture whose properties are well-studied in the literature. Lastly, under the finite support assumption, the identification of the DTE parameters reduces down to a GMM model with quadratic moments, giving us some insights on how the DTE parameters are identified. Though I assume that $U$ is finitely discrete in the estimation, the identification result does not require such a restriction and I develop an alternative estimation method based on sieve maximum likelihood for a setup with continuous $U$, in the Appendix subsection (ref).
The estimation procedure is two-step. In the first step, I estimate the conditional probability $\Pr \{U=u|Z=z\}$, using the nonnegative matrix factorization. In the second step, I identify a DTE parameter with a moment condition involving probabilities $\Pr \{Y \leq y | D=d, Z=z\}$ and $\Pr \{U=u|Z=z\}$. The former probability is directly observed from the dataset and the latter is estimated in the first step. Thus, the estimator can be thought of as a plug-in GMM estimator, where nuisance parameters are estimated in the first-step nonnegative matrix factorization. Asymptotic normality of the distributional treatment effect parameters is established. In deriving asymptotic normality, I construct a moment condition that satisfies Neyman orthogonality to be robust to the first-step estimation error from the nonnegative matrix factorzation.
This paper makes contribution to the distributional treatment effect literature by proposing a framework where the joint distribution of the potential outcomes and thus the marginal distribution of treatment effect are point identified, without imposing any functional form assumptions. This is in contrast to the partial identification results in the literature: FP2010,FSS2014,FR2019,FL2021 and more. There exist several notable point identification results: HSC,CHH. These point identification results either assume independence on potential outcomes conditioning on observables only, or assume structural assumptions on treatment and/or potential outcomes: e.g. a Roy model for treatment and a factor structure for potential outcomes. In terms of estimation, WP2006,N2023 also develop DTE estimators; unlike this paper, they both build on the point identification result without latent conditioning variable and develop a deconvolution-based estimator. BBDM provides an insightful overview on recent developments on both identification and estimation in program evaluation literature regarding distributional concerns.
This paper also makes contribution to the nonclassical measurement error/proximal inference literature and the finite mixture literature: HS2008,HKS2014,MGT,D2023,K2023,N2022 and more. In terms of identification, the latent rank invariance assumption provides an alternative assumption in labeling the latent variable that uses the information from the outcome variable $Y$, as opposed to the “measure of location” assumption suggested in HS2008 that uses the information from the proxy variable $X$. Alternatively, this paper can be thought of as adding an additional identifying assumption\textemdash conditional independence between $Y(d)$ and $X$\textemdash to narrow down the identified set of HKS2014 to a singleton. In terms of the asymptotic theory on the estimator, this paper is in a similar setup as H2008, assuming a finite support. Unlike existing estimators based on the principal component anaylsis, the estimation strategy based on the nonnegative matrix factorization as proposed in this paper has guarantee that the estimated conditional distributions are indeed nonnegative and sum-to-one. The $\sqrt{n}$-consistency is proven in this paper so that future works may build upon the nonnegative matrix factorization estimator.
The rest of the paper is organized as follows. Section (ref) discusses the identification result for the joint distribution of the two potential outcomes. Section (ref) explains the estimation method for the two DTE parameters and develops asymptotic theory for the estimators. Section (ref) contains Monte Carlo simulation restuls and Section (ref) applies the estimation procedure to an empirical dataset from JMR.
An econometrican observes a dataset $\left\lbrace Y_{i}, D_i, X_i, Z_i \right\rbrace_{i=1}^n$ where $Y_{i}, X_i, Z_i \in \mathbb{R}$ and $D_i \in \{0,1\}$. $Y_i$ is an outcome variable, $D_i$ is a binary treatment variable and $X_i, Z_i$ are two proxy variables. The outcome $Y_{i}$ is constructed with two potential outcomes.
In addition to $\left(Y_i(1), Y_i(0), D_i, X_i, Z_i\right)$, there is a latent variable $U_i \in \mathcal{U} \subset \mathbb{R}$. $U_i$ plays a key role in putting restrictions on the joint distribution of $Y_i(1)$ and $Y_i(0)$ and overcoming the fundamental limitation that we observe only one potential outcome for a given unit. The dataset comes from random sampling: $\left( Y_i(1), Y_i(0), D_i, X_i, Z_i, U_i \right) \overset{iid}{\sim} \mathcal{F}$.
Firstly, I assume conditional random assignment on the treatment $D_i$ and exclusion restriction on the proxy variable $Z_i$.
Assumption (ref) assumes that the treatment is as good as random with regard to the potential outcomes and $X_i$ after conditioning on the latent variable $U_i$. In this sense, Assumption (ref) is a restriction on treatment endogeneity. In addition, Assumption (ref) assumes that the proxy variable $Z_i$ does not have any additional information on the potential outcomes after conditioning on the latent variable $U_i$, satisfying exclusion restriction. Note that Assumption (ref) does not impose any restriction on the dependence between $Z_i$ and $D_i$. The proxy variable $Z_i$ may still depend on treatment. This assumption imposes nontrivial restriction on treatment endogeneity in a non-experimental context. Thus, for the rest of the paper, I focus on a randomly assigned treatment, limiting my attention to randomized controlled trials. The following condition is a sufficient condition for Assumption (ref).
Remark. A sufficient condition for Assumption (ref) is
Thus, in a randomized controlled trial setup, Assumption (ref) is satisfied when there are one proxy variable that is independent of the treatment $D_i$ and another proxy variable that only depends the latent heterogeneity $U_i$ and the treatment $D_i$.
When $U_i$ is observed, Assumption (ref) identifies numerous treatment effect parameters such as average treatment effect (ATE), quantile treatment effect (QTE) and more. However, even when $U_i$ is observed, we still cannot identify the distribution of treatment effect from Assumption (ref) since Assumption (ref) does not tell us anything about the dependence between $Y_i(1)$ and $Y_i(0)$.
To impose restrictions on the joint distribution of $Y_i(1)$ and $Y_i(0)$ and have more identifying power, I assume that the latent variable $U_i$ captures all of the dependence between the two potential outcomes and the proxy variable $X_i$.
Assumption (ref) assumes that the two potential outcomes $Y_i(1)$ and $Y_i(0)$ and the proxy variable $X_i$ are mutually independent of each other conditioning on $U_i$. Given $U_i$, the proxy variable $X_i$ does not give us additional information on the distribution of the potential outcomes. Note that the latent variable $U_i$ lies in $\mathbb{R}$ as do $Y_i(1)$ and $Y_i(0)$. This excludes a non-binding case where $U_i = \big( Y_i(1), Y_i(0) \big)$.
When $U_i$ is observed, Assumptions (ref)-(ref) identify the joint distribution of the two potential outcomes and various distributional treatment effect parameters. Examples include the variance of the treatment effect $\mathrm{Var} \big( Y_i(1) - Y_i(0) \big)$ and the marginal distribution of the treatment effect $\Pr \{Y_i(1) - Y_i(0) \leq \delta \}$. Since $U_i$ is not observed, identifying the conditional densities of $Y_i(1), Y_i(0)$ given $U_i$ and the marginal density of $U_i$ will be the main challenge in the identification.
Assumptions (ref)-(ref) play a key role in the identification result. Below I present three examples of econometric frameworks that motivate Assumptions (ref)-(ref). The first example is rank invariance, which is widely used in the quantile treatment effect literature and the quantile IV literature: see CH2005,CH2006,AI2006,VX2017,CL2019,HX2023 and more.
The usage of this rank invariance assumption is mostly limited to the quantile treatment effect and not applied to the distributional treatment effect, due to the fact that it imposes excessive restriction on the joint distribution of the two potential outcomes. Under the rank invariance, $\mathrm{Var} \big( Y_i(1)|Y_i(0) \big)=0$ and vice versa. Thus, the treatment effect $Y_i(1) - Y_i(0)$ is also a deterministic function of $Y_i(1)$ and of $Y_i(0)$.
In this paper, I relax this deterministic relationship among the potential outcomes and the latent variable, by assuming the rank invariance not on the potential outcomes directly, but on some functional of the conditional distribution of the potential outcome given $U_i$. In this sense, the econometric framework of this paper is a relaxation of the rank invariance assumption in the quantile treatment effect literature and the quantile IV literature.
The second example is a panel data model with a latent state variable that is first-order Markovian, which is often referred to as a hidden Markov model. The hidden Markov model is widely discussed in the dynamic panel data model literature, especially in the context of the dynamic discrete choice and conditional choice probability estimation: see KS2009,AM2011,HS2012,HY2018 for more.
In this nonlinear panel data model, the potential outcome $Y_{it}(d)$ is a function of a latent variable $V_{it}$ and an error term $\varepsilon_{it}(d)$. Note that $V_{it}$ appears in the model twice; for $Y_{it}(1)$ and for $Y_{it}(0)$. In this sense, $V_{it}$ is a common shock to the potential outcomes where $\varepsilon_{it}(d)$ is a treatment-status-specific shock. The key elements of Example (ref) are that the common shock process $\{V_{it}\}_{t=1}^3$ and the treatment-status-specific shocks $\varepsilon_{i1}(1), \ldots, \varepsilon_{i3}(0)$ are all mutually independent and that dependence within $\left\lbrace V_{it} \right\rbrace_{t=1}^3$ themselves is restricted to be first-order Markovian given $D_i$. Thus, $V_{i2}$ has sufficient information on the dependence between $Y_{i2}(1)$ and $Y_{i2}(0)$ and the past and the future outcomes $Y_{i1}$ and $Y_{i3}$ can be used as proxies for $V_{i2}$.
The hidden Markov model is mostly applied to a single observed outcome setup. In this paper, I extend the hidden Markov model to a potential outcome setup so that there are two idiosyncratic error terms $\varepsilon_{it}(0)$ and $\varepsilon_{it}(1)$, specific to each treatment status. Moreover, I add one more conditional independence to the hidden Markov model by assuming that the two error terms are independent across the treatment status conditioning on $U_i$.
The second example closely relates to Section (ref) of this paper. In Section (ref), I revisit JMR and estimate the full distribution of treatment effect, in the empirical context of workplace wellness program as a treatment and monthly medical spending as an outcome. The dataset used in JMR contains short panel data on monthly medical spending, with one pretreatment time period. Thus, by assuming that the monthly medical spending is a function of two different types of random shocks\textemdash a transitory, idiosyncratic shock and a systemic health shock that is first-order Markovian\textemdash, the model described in Example (ref) can be applied to the dataset and we can use the pretreatment medical spending and the post-treatment medical spending as the two proxy variables.
The third example is where we have economic interpretation on the latent variable $U_i$ and therefore can find measurements on the latent variable. There are several notable papers in labor economics that adopts this approach: see CHH,CH2008,CHS and more.
The above model is a simplified version of the framework in CHS, applied to a potential outcome setup. In this example, an economic model gives us an interpretation on the latent variable $U_i$ and helps us find measurements on the latent variable. For example, in CHS, $U_{X,i}$ and $U_{Z,i}$ are assumed to be cognitive skill and noncognotive skill. Then, various measures on cognitive ability, temperament, motor and social developments and such are used as proxy variables. In this paper, I consider a more flexible outcome function than the CES function, at the cost of assuming a univariate $U_i$ and a strong independence assumption on the error terms.\footnote{The main focus of CHS is less on the outcome, but more on the skill formation. In the full framework of CHS, there exists time dimension and the skills vector $U_{it}$ is modeled with a dynamic process and the paper nonparametrically identifies the skill evolution process.} In this sense, this paper can also be thought of as nonparametric version of the CHH's framework.
The remainder of this section outlines the identification argument. For illustration purposes only, let $Y_i, X_i, Z_i, U_i$ be discrete: $Y_i \in \{y^1, \cdots, y^{M_Y}\}, X_i \in \{x^1, \cdots, x^{M_X}\}, Z_i \in \{z^1, \cdots, z^{M_Z}\}$ and $U_i \in \{u^1, \cdots, u^K\}$. With $M = M_Y \cdot M_X$, we can construct a $M \times M_Z$ matrix $\mathbf{H}_d$ of conditional probabilities as follows:
for each $d=0,1$. $\mathbf{H}_0$ is the conditional probability of $(Y_i, X_i)$ given $Z_i$ in the untreated subsample and $\mathbf{H}_1$ is the conditional probability in the treated subsample. From Assumptions (ref)-(ref), both $\mathbf{H}_0$ and $\mathbf{H}_1$ decompose into a multiplication of two matrices: for each $d=0,1$,
where
Note that the discreteness of $Y_i, X_i, Z_i$ is nonbinding; we can use partitioning on $\mathbb{R}$ when they are continuous.\footnote{Consider partitions on $\mathbb{R}$ such that $$ \left\lbrace \mathcal{Y}^m = \left(y^{m-1},y^m \right] \right\rbrace_{m=1}^{M_Y}, \hspace{5mm} \left\lbrace \mathcal{X}^m = \left(x^{m-1},x^m \right] \right\rbrace_{m=1}^{M_X}, \hspace{5mm} \left\lbrace \mathcal{Z}^m= \left(z^{m-1},z^m \right] \right\rbrace_{m=1}^{M_Z} $$ where $y^0=x^0=z^0=-\infty$ and $y^{M_Y} = x^{M_X} = z^{M_Z} = \infty$. Let $\mathcal{W}^1=\mathcal{Y}^1 \times \mathcal{X}^1, \mathcal{W}^2 = \mathcal{Y}^2 \times \mathcal{X}^1, \cdots, \mathcal{W}^M = \mathcal{Y}^{M_Y} \cdot \mathcal{X}^{M_X}$. $\left\lbrace \mathcal{W}^m \right\rbrace_{m=1}^M$ is a partition on $\mathbb{R}^2$. Then, $\mathbf{H}_d$ becomes
for each $d=0,1$. $\Gamma_d$ and $\Lambda_d$ are similarly constructed with partitioned $Y_i, X_i$ and $Z_i$.}The remaining discretization on $U_i$ is imposed only for the expositional brevity; the identification argument does not hinge on the discreteness of $U_i$. The continuous $U_i$ version of the identification follows the same argument and uses one additional assumption to find a labeling on the infinite number of functions: Assumption (ref). I present more discussion on Assumption (ref) later in this section and a full identification argument for continuous $U_i$ is provided in Subsection (ref) of Appendix.
The equation $\mathbf{H}_d = \Gamma_d \cdot \Lambda_d$ shows us that the conditional density model in (ref) is indeed a mixture model. For each subpopulation $\left\lbrace i: (D_i, Z_i) = (d,z)\right\rbrace$, there is a column in the matrix $\Lambda_d$ which denotes the subpopulation-specific distribution of $U_i$. Then, the density of $\left(Y_i, X_i \right)$ in that subpopulation admits a mixture model with the aforementioned columns of $\Lambda_d$ as mixture weights and the conditional density of $\left(Y_i(d), X_i \right)$ given $U_i$ as mixture component densities. The equation $\mathbf{H}_d = \Gamma_d \cdot \Lambda_d$ aggregates the finite mixture formulations across the subpopulations.
Note that from Assumption (ref), the joint distribution of $Y_i(1)$ and $Y_i(0)$ is identified if the conditional distribution of $Y_i(1)$ given $U_i$, the conditional distribution of $Y_i(0)$ given $U_i$, and the marginal distribution of $U_i$ are identified. The first two distributions correspond to $\Gamma_1$ and $\Gamma_0$ in the discretization. The last distribution is a function of $\Lambda_1, \Lambda_0$ and the distribution of $(D_i, Z_i)$, which is observed. Thus, to identify of the distributional treatment effect parameter is to identify $\Gamma_1, \Gamma_0, \Lambda_1$ and $\Lambda_0$.
To decompose $\mathbf{H}_d$ into $\Gamma_d$ and $\Lambda_d$, first fix $y \in \{y^1,\cdots,y^{M_Y}\}$ and extract rows of $\mathbf{H}_d$ and $\Gamma_d$ that correspond to $(y, x^1), \cdots, (y, x^{M_X})$:
for $d=0,1$. From Assumption (ref), the mixture component density matrix $\Gamma_d(y)$ can be further decomposed:
Now, sum $\mathbf{H}_d(y)$ across $y^1, \cdots, y^{M_Y}$:
Find that when $M_X=M_Z=K$ and both $\Gamma_X$ and $\Lambda_d$ have full rank,
Given a no repeated eigenvalue condition that for any $u \neq u'$ there exist some $(y, d)$ such that $\Pr \{Y_i(d)=y | U_i=u\} \neq \Pr \{Y_i(d)=y | U_i=u'\}$, diagonalization of $\mathbf{H}_d(y) \left( \sum_y \mathbf{H}_d(y) \right)^{-1}$ across different $y$ and $d$ identifies $\Gamma_X$ and $\{\Delta_d(y)\}_{y^1 \leq y \leq y^{M_Y}}$.\footnote{Eigenvalue decomposition on its own is not unique but we have sufficiently many constraints on $\Gamma_X$ for unqiueness; $\Gamma_X$, the eigenvector matrix, is nonnegative and its column-wise sums are one since they are conditional probabilities. See HS2008 for more.} Once $\Gamma_X$ is identified, the identification of $\Lambda_0, \Lambda_1$ follows from $\Gamma_X$ having full rank. When $M_X$ or $M_Z$ is bigger than $K$, we may stack some of the rows or the columns of $\sum_y \mathbf{H}_d(y)$ to make it into a square matrix.
Assumption (ref) formally states the full rank condition and the no repeated eigenvalue condition for discrete $U_i$.
Assumption (ref).b implicitly assumes that $M_X, M_Z \geq K$. The restriction that $M_X, M_Z \geq K$ is sensible since I use the variation in the conditional density of $X_i$ given $Z_i=z$ across $z$ to capture the variation in the latent variable $U_i$. The support for the two proxy variables has to be at least as rich as the support of the latent variable. Assumption (ref).c assumes that the eigenvalue decomposition does not have repeated eigenvalues.
Assumption (ref) reiterates Assumption (ref) for a setup where $U_i$ are continuous. Let $f_{Y(d)|U}$ denote the conditional density of $Y_i(d)$ given $U_i$, $f_{X|U}$ denote the conditional density of $X_i$ given $U_i$, and $f_{U|D=d,Z}$ denote the conditional density of $U_i$ given $D_i=d$ and $Z_i$, for $d=0,1$. Define integral operators $L_{X|U}$ and $L_{U|D=d,Z}$ that map a function in $\mathcal{L}^1(\mathbb{R})$ to a function in $\mathcal{L}^1(\mathbb{R})$: for $d=0,1$,
Assumption (ref).c corresponds to Assumption (ref).b and Assumption (ref).d to Assumption (ref).c.
When $U_i$ is continuous, we need an addtional assumption for the identification. This is because when $U_i$ is discrete and finite, a bijection between $u$ and $\Pr \{X_i= \cdot | U_i=u\}$ needs not be specified. However, when $U_i$ is continuous, we need an ordering on the infinite collection $\{f_{X|U}(\cdot | u) \}_u$ to connect $u$ to $f_{X|U}(\cdot | u)$.
The functional $M$ provides us an ordering on the infinite collection $\{ f_{X|U} (\cdot | u) \}_{u}$, by applying the functional to $\{ f_{Y(1)|U}(\cdot |u), f_{Y(0)|U}(\cdot |u) \}_u$. A simple example where Assumption (ref) fails is when $\mathcal{U} = [-1,1]$ and $Y_i(d) | U_i=u \sim \mathcal{N} ( u^2 + d, \sigma^2)$. Neither $f_{Y(1)|U}$ nor $f_{Y(0)|U}$ helps us find an ordering between $f_{X|U} (\cdot | u)$ and $f_{X|U}(\cdot | -u)$.
Along with Assumptions (ref)-(ref), Assumption (ref) is a key identifying assumption in the case of continuous $U_i$. As hinted by its label, Assumption (ref) draws the inspiration from the rank invariance assumption in Example (ref). Suppose that Assumption (ref) holds true for $\mathbf{E} \left[ Y_i(1)|U_i=u \right]$ and $\mathbf{E} \left[ Y_i(0)|U_i=u \right]$. Then, the two potential outcomes of a given unit have the same `latent rank' in the sense that their expected values $\mathbf{E} \left[ Y_i(1) | U_i \right]$ and $\mathbf{E} \left[ Y_i(0) | U_i \right]$ have the same rank in their respective distributions. A similar assumption can be made with other summary measures such as median or mode. Recall that Assumptions (ref)-(ref) is a relaxation of the rank invariance assumption. Assumption (ref) allows us to retain the rank interpretation on the latent variable $U_i$. Conditioning on the latent variable $U_i$, some summary measure applied to the conditional distributions of the potential outcomes has the same rank.
Theorem (ref) formally states the identification result.
The result of Theorem (ref) can be understood as applying the identification result of HS2008 twice, once to the treated population and again to the untreated population, and then connecting the two identification results. Also, when $U_i$ is finite, the result of Theorem (ref) can be understood as a point identification adaptation of the partial identifaction result from HKS2014; the additional identifying power comes from the conditional independence between $Y_i(d)$ and $X_i$ given $U_i$.
It directly follows Theorem (ref) that any functional of the joint distribution of $Y_i(1)$ and $Y_i(0)$ is identified: e.g., $\mathrm{Var} \big( Y_i(1) - Y_i(0) \big), \Pr \left\lbrace Y_i(1) \geq Y_i(0) \right\rbrace, \Pr \left\lbrace Y_i(1) \geq Y_i(0) | Y_i(0) \right\rbrace$ and etc. The rest of the section discusses the restrictions on the joint distribution of $Y_i(1)$ and $Y_i(0)$ implied by the identifying assumptions and a testable implication of the identifying assumptions which proposes a falsification test.
Assumption (ref) assumes that there exists a latent variable $U_i$ which contains sufficient information on the dependence between a treated potential outcome and an untreated potential outcome. Assumption (ref).b and Assumption (ref).c assume that the proxy variable $X_i$ and $Z_i$ create sufficient variation in the latent variable $U_i$. By assuming $X_i, Z_i$ and $U_i$ are scalar variables, I exclude the trivial case of $U_i = \big(Y_i(1), Y_i(0) \big)$ and impose implicit restrictions on the joint distribution of $Y_i(1)$ and $Y_i(0)$.
To discuss the implicit restrictions imposed by the identifying assumptions, let us consider a simple quantity of $\mathbf{E}[Y_i(1) Y_i(0)]$. $\mathbf{E}[Y_i(1) Y_i(0)]$ is a key ingredient in identifying $\mathrm{Var} \big( Y_i(1) - Y_i(0) \big)$, a measure of the treatment effect heterogeneity. In most econometric frameworks that identify ATE or QTE, $\mathbf{E}[Y_i(1) Y_i(0)]$ still remains unidentified. In this paper, using Assumption (ref).b or Assumption (ref).c, the conditional density of $Y_i(1)$ given $Y_i(0)$ is identified as a weighted average of the conditional densities of $Y_i$ given $(D_i=1,Z_i)$, identifying $\mathbf{E}[Y_i(1) Y_i(0)]$. The core idea in constructing the weights is that the conditional density $f_{Y(1)|U}$ is identified as a weighted average of $\{f_{Y|D=1,Z}(\cdot|z)\}_z$, from the completeness of $L_{U|D=1,Z}$. With $w(\cdot, \cdot)$ denoting the weighting function, $$ f_{Y(1)|Y(0)}(\cdot | y) = \int_{\mathbb{R}} \frac{w(y,z)}{f_{Y(0)}(y)} \cdot f_{Y|D=1,Z}(\cdot|z) dz $$ and thus $$ \mathbf{E}[Y_i(1)|Y_i(0)=y] = \int_{\mathbb{R}} \frac{w(y,z)}{f_{Y(0)}(y)} \cdot \mathbf{E}[Y_i |D_i=1,Z_i=z] dz. $$ $\mathbf{E}[Y_i(1) Y_i(0)]$ is identified as
Note that $\mathbf{E}[Y_i(1) Y_i(0)]$ is identified as an expected product of two random variables $Y_i(0)$ and $\mathbf{E}[Y_i|D_i=1,Z_i]$, reweighted with $\frac{w}{f_{Y(0),Z}}$. Even though we do not observe $Y_i(1)$ and $Y_i(0)$ simultaneously, the result above shows us that we can instead use $\mathbf{E}[Y_i|D_i=1,Z_i]$, a random variable that is observed for every untreated unit, in place of $Y_i(1)$ and reweight the joint density of $Y_i(0)$ and $Z_i$ with $w(\cdot, \cdot)$. Thus, the implicit restriction in identifying $\mathbf{E}[Y_i(1) Y_i(0)]$ is that the conditional expectation $\mathbf{E}[Y_i(1) |Y_i(0)=y]$ must be spanned by the observed conditional expectations $\{\mathbf{E}[Y_i|D_i=1,Z_i=z]\}_z$. Since the above identification argument can be rewritten with $\mathbf{E}[Y_i(0)|Y_i(1)=y]$, another implicit restriction is that the conditional expectation $\mathbf{E}[Y_i(0)|Y_i(1)=y]$ must be spanned by the observed conditional expectations $\{\mathbf{E}[Y_i|D_i=0,Z_i=z]\}_z$. The identifcation argument can also be extended to conditional densities, instead of conditional expectations; thus, more generally, the implicit restriction imposed on the joint distribution of $Y_i(1)$ and $Y_i(0)$ is that the conditional distribution of $Y_i(1)$ given $Y_i(0)$ must be spanned by the conditioanl distribution of $Y_i$ given $(D_i=1,Z_i)$ and vice versa.
When we extend Assumption (ref) so that both $u \mapsto M f_{Y(1)|U}(\cdot|u)$ and $u \mapsto M f_{Y(0)|U}(\cdot|u)$ are strictly increasing and continuously differentiable, we have a testable implication of Assumptions (ref)-(ref) and (ref)-(ref), from over-identification. Suppose that $\mathbf{E} \left[ Y_i(1)|U_i=u \right]$ and $\mathbf{E} \left[ Y_i(0)|U_i=u \right]$ are strictly increasing in $u$. Then, the conditional densities $\left( f_{Y(1)|U}, f_{X|U}, f_{U|D=1,Z} \right)$ are identified in the treated subsample and the conditional densities $\left( f_{Y(0)|U}, f_{X|U}, f_{U|D=0,Z} \right)$ are identified in the untreated subsample. Let $f_{X|D=1,U}$ denote the conditional density of $X_i$ given $U_i$, identified from the treated subsample and likewise for $f_{X|D=0,U}$. Then, Assumption (ref) imposes that
since $f_{X|D=1,U} = f_{X|D=0,U}$. In (ref), a monotone function $\tilde{g}$ is used to connect the identification result from the treated subpopulation to the untreated subpopulation, now that $f_{X|U}$ is not used to connect the two identification results. A test that uses (ref) as a null can be used as a falsification test on the framework proposed in this paper.
What does a test on the null (ref) exactly test? The mixture model on the conditional density $f_{Y,X|D=d,Z}$ assumes that conditioning on $U_i$, the potential outcome $Y_i(d)$ and the proxy variable $X_i$ are independent of each other. Recall that in Example (ref), the proxy variable $X_i$ is a past outcome. Thus, in the panel context, we can understand the falsification test as testing whether we can find a latent variable $U_i$ conditioning on which the outcomes are intertemporally independent. Note that the key identifying assumption is that the potential outcomes independent across the treatment status. While the conditional independence assumption across the treatment status remains untestable due to the limitation that we only observe either a treated potential outcome or a untreated potential outcome for a given unit, the falsification test in Example (ref) tests if the outcomes are intertemporally independent, conditioning on some latent variable.
In the case of discrete $U_i$, Assumption (ref) was not used in the identification. In fact, without introducing any further assumptions, we have a testable implication:
I develop an asymptotic theory in the next section under the finite support assumption on $U_i$, formally proposing a falsification test.
Based on the identification result for discrete $U_i$, I estimate the conditional density of $Y_i(1)$ and $Y_i(0)$ given $U_i$, by assuming a finite support for $U_i$ and solving a nonnegative matrix factorization (NMF) problem. The focus on the case of discrete $U_i$ has several reasons. Firstly, a dicretization is often used in econometric models with latent heterogeneity as an approximation to a continuous latent heterogeneity space: see bonhomme2022discretizing for more. Secondly, with parametrization, the estimation of infinite-dimensional objects such as conditional densities $f_{U|D=0,Z}$ and $f_{U|D=1,Z}$ becomes an estimation of finite-dimensional objects $\Lambda_0$ and $\Lambda_1$, giving us $\sqrt{n}$ rate. The $\sqrt{n}$ rate becomes helpful in deriving an asymptotic distribution for the distributional treatment effect estimators. Lastly, the linearity induced from discretization reduces the computational burden substantially. This does not mean that there is no feasible estimation method for continuous $U_i$. For a continuous latent variable case, we can construct a sieve maximum likelihood estimator, as suggested in the nonclassical measurement error literature. The specifics are discussed in the appendix subsection (ref).
The parameters of interest in this paper are the joint distribution of the potential outcomes $Y_i(1)$ and $Y_i(0)$ and the marginal distribution of the treatment effect $Y_i(1) - Y_i(0)$. To estimate these distributional treatment effect (DTE) parameters, I first estimate the conditional probabilites of $U_i$ given $Z_i$, namely the mixture weight matrices $\Lambda_0$ and $\Lambda_1$ in the finite mixture interpretation, by solving a nonnegative matrix factorization problem. Given the first step estimators on $\Lambda_0$ and $\Lambda_1$, I characterize the the joint distribution of $Y_i(1)$ and $Y_i(0)$ and the marginal distribution of $Y_i(1) - Y_i(0)$ as quadratic moments and estimate the distributions by plugging in the first step estimates to the induced $U$-statistics. In doing so, to account for the estimation error from the first step, I orthogonalize the score function. Neyman orthogonality makes the plug-in estimator robust to the first step estimation error and helps derive a limiting distribution for the estimator.
To estimate the mixture weight matrices $\Lambda_0$ and $\Lambda_1$ from (ref), I first let $M_Z=K$ by using a partition on $\mathbb{R}$ when the support of $Z_i$ has more than $K$ points and construct sample analogues of the conditional probability matrices $\mathbf{H}_0$ and $\mathbf{H}_1$ defined in the previous section: for $d=0,1$, let
Each column of $\mathbb{H}_0$ is a conditional empirical distribution function of $\left(Y_i, X_i \right)$ given $\big(D_i=0, Z_i \big)$ and each column of $\mathbb{H}_1$ is a conditional empirical distribution function of $\left(Y_i, X_i \right)$ given $\big(D_i=1, Z_i \big)$. As discussed in Section (ref), I use partitioning on $\mathbb{R}$ in constructing $\mathbb{H}_0$ and $\mathbb{H}_1$ when any of $Y_i, X_i$ and $Z_i$ is continuous.
To estimate $\Lambda_0$ and $\Lambda_1$, I formulate a nonnegative matrix factorization problem. Let $\iota_x$ be a $x$-dimensional column vector of ones. Then, the nonnegative matrix factorization problem is constructed as follows:
subject to linear constraints that
and quadratic constraints that
for each $(y,x)$. The linear constraints are probabilities being nonnegative and summing to one. The quadratic constraints are $X_i$ satisfying the exclusion restriction from Assumption (ref). When $\mathbb{H}_0$ and $\mathbb{H}_1$ are sufficiently close to $\mathbf{H}_0$ and $\mathbf{H}_1$, the identification result discussed in the previous section says that there is a unique decomposition of $\mathbb{H}_0$ and $\mathbb{H}_1$ which satisfies the linear and the quadratic constraints.
Note that the objective function in (ref) is quadratic when we fix either $\left(\Lambda_0, \Lambda_1 \right)$ or $\left(\Gamma_0, \Gamma_1 \right)$. Moreover, $\Gamma_0$ and $\Gamma_1$ can be further decomposed into three matrices $\Gamma_{X}, \Gamma_{Y(0)}, \Gamma_{Y(1)}$, each of which correponds to the conditional probabilities of $X_i$ given $U_i$, $Y_i(0)$ given $U_i$, and $Y_i(1)$ given $U_i$, respectively. Let $\Gamma_d(\cdot, \cdot)$ denote how $\Gamma_X$ and $\Gamma_{Y(d)}$ recover $\Gamma_d$: $\Gamma_d = \Gamma_d \Big( \Gamma_X, \Gamma_{Y(d)} \Big)$. The quadratic constraints are trivially imposed by optimizing over $\Gamma_X, \Gamma_{Y(0)}$ and $\Gamma_{Y(1)}$. Using these, I propose an iterative algorithm to solve the minimization problem.
Each step of the iteration is a quadratic programming with linear constraints, which can be solved with a built-in optimization tool in most statistical softwares. The stepwise optimization assures a convergence to a local minimum. To find the global minimum, I consider various initial values $\left( \Gamma_0^{(0)}, \Gamma_1^{(0)} \right)$.\footnote{To initialize $\Gamma_0^{(0)}, \Gamma_1^{(1)}$, I consider columns from $\mathbb{H}_d$ and weighted sums of columns of $\mathbb{H}_d$ with randomly drawn $K$ sets of weights that sum to one as initial values. Alternatively, we can select the eigenvectors associated with the first $K$ largest eigenvalues of ${\mathbb{H}_d}^\intercal \mathbb{H}_d$ as an initial value.}
Let $\widehat{\Lambda}_0$, $\widehat{\Lambda}_1$, $\widehat{\Gamma}_0$ and $\widehat{\Gamma}_1$ denote the solution to the minimization problem. Note that when $Y_i$ and $X_i$ are discrete, the estimates $\widehat{\Gamma}_0$ and $\widehat{\Gamma}_1$ directly estimate the conditional distribution of $Y_i(1)$ and $Y_i(0)$ given $U_i$. When $Y_i$ are $X_i$ are continuous and therefore partitioning was used in constructing $\mathbf{H}_0, \mathbf{H}_1$, we use $\widehat{\Lambda}_0$ and $\widehat{\Lambda}_1$ to estimate the distribution of $Y_i(1)$ and $Y_i(0)$ given $U_i$.
Given the estimates of the two mixture weights matrices $\Lambda_0$ and $\Lambda_1$, I construct an estimator for the joint distribution of $Y_i(1)$ and $Y_i(0)$ and the marginal distribution of $Y_i(1) - Y_i(0)$. Firstly, find that for any $y \in \mathbb{R}$,
Since $\Lambda_d$ is full rank, we have
The conditional distribution of $F_{Y(d)|U}(\cdot|u)$ is identified as a linear combination of the observed distributions $\{F_{Y|D=d,Z}(\cdot | z)\}_{z=1}^{K}$. Building on this, let $$ \tilde{\Lambda}_d = \left( \Lambda_d \right)^{-1} $$ for $d=0,1$. Let $\tilde{\lambda}_{jk,d}$ denote the $j$-th row and $k$-th column component of $\tilde{\Lambda}_d$. $\left(\tilde{\lambda}_{1k,d}, \cdots, \tilde{\lambda}_{K k,d} \right)^\intercal$, the $k$-th column of $\tilde{\Lambda}_d$, is a set of linear coefficients on $\{F_{Y|D=d,Z}(\cdot | z)\}_{z=1}^{K}$ to retrieve the conditional distribution of $Y_i(d)$ given $U_i=u^k$. Using the estimators on $\Lambda_0, \Lambda_1$ from the nonnegative matrix factorization, we estimate the linear coefficients as follows: $$ \widehat{\tilde{\Lambda}}_d = \left( \widehat{\Lambda}_d \right)^{-1} $$ for $d=0,1$.
Secondly, the distribution of $U_i$ is also identified from $\Lambda_0$ and $\Lambda_1$:
Let $p_{U}(k)$ denote $\Pr \{U_i = u^k\}$ for $k=1, \cdots, K$ and let $p_{D,Z}(d,j)$ denote $\Pr \{D_i=d,Z_i = z^j\}$ for $d=0,1$ and $j=1, \cdots, K$. Then, I estimate $p_U$ and $p_{D,Z}$ with $$ \hat{p}_{D,Z}(d,j) = \frac{1}{n} \sum_{i=1}^n \mathbf{1}\{D_i=d,Z_i=z^j\} $$ and $$ \hat{p}_U =
= \widehat{\Lambda}_0
+ \widehat{\Lambda}_1
. $$
By combining the two results, we get
Using this characterization, I estimate the joint distribution of $Y_i(1)$ and $Y_i(0)$ as a linear combination of $\{F_{Y|D=0,Z}(y|z^j) \cdot F_{Y|D=1,Z}(y'|z^{j'})\}_{j,j'}$ where the weights are computed with $\widehat{\Lambda}_0, \widehat{\Lambda}_1$ and $\{\hat{p}_{D,Z}(d,j)\}_{d,j}$. We can derive a similar result for the marginal distribution of $Y_i(1) - Y_i(0)$: for any $\delta \in \mathbb{R}$,
Both parameters of interest are identified as a weighted sum of quantities that are indexed by pairs of subpopulations $\{i: D_i=0, Z_i=z^j\}$ and $\{i:D_i=1, Z_i=z^{j'}\}$. As shown above, weights are estimated from the first step nonnegative matrix factorization and empirical measures of the subpopulations. It remains to estimate the quantities associated with each pair of subpopulations. I will discuss this for the marginal distribution of $Y_i(1) - Y_i(0)$; the case for the joint distribution of $Y_i(1)$ and $Y_i(0)$ follows naturally. For some $\delta$, let $$ \theta = F_{Y(1)-Y(0)}(\delta). $$
Firstly, find that $\theta$ is a summation over $K$ treated subpopulations and $K$ untreated subpopulations. Fix $j, j'$ and let
Find that
with $\left( Y_i, D_i, Z_i\right) \perp \!\!\! \perp \left( Y_{i'}, D_{i'}, Z_{i'}\right)$. Thus, $\theta_{jj'}$ is identfied from a quadratic moment $$ \mathbf{E} \left[m_{jj'} \left(W_i, W_{i'} ; \theta_{jj'},\tilde{\Lambda}_0, \tilde{\Lambda}_1, \{p_U(k)\}_k, \{p_{D,Z}(d,j)\}_{d,j} \right) \right]=0 $$ where $W_i= (Y_i, D_i, X_i, Z_i)$ and
By summing over $j$ and $j'$, we can construct a moment function $m = \sum_{j=1}^K \sum_{j'=1}^K m_{jj'}$ such that $$ \mathbf{E} \left[m \left(W_i,W_{i'} ; \theta, \tilde{\Lambda}_0, \tilde{\Lambda}_1, \{p_U(k)\}_k, \{p_{D,Z}(d,j)\}_{d,j} \right) \right]=0 $$ identifies $\theta$.
If the nuisance parameters $\tilde{\Lambda}_0, \tilde{\Lambda}_1, p_U, p_{D,Z}$ were known, the standard asymptotic theory of $U$ statistic would apply to the GMM estimator of $\theta$ using $\mathbf{E}[m(W_i, W_{i'}; \theta)]=0$ as the moment condition. However, in practice, we use first step estimates for the nuisance parameters. Thus, to account for the first step estimation error, we orthogonalize the moment function. Even though the NMF estimators $\left( \widehat{\Lambda}_0, \widehat{\Lambda}_1 \right)$ and the induced estimators $\left( \widehat{\tilde{\Lambda}}_0, \widehat{\tilde{\Lambda}}_1 \right)$ are complex nonlinear functions of the data matrix $\mathbb{H}_0$ and $\mathbb{H}_1$, $\left( \tilde{\Lambda}_0, \tilde{\Lambda}_1 \right)$ satisfy the following equations at their true values:
Equation (ref) corresponds to the conditional independence assumption that $$ \Pr \{Y_i(d)=y, X_i=x|U_i=u\} = \Pr \{Y_i(d)=y|U_i=u\} \cdot \Pr \{X_i=x|U_i=u\}. $$ and Equation (ref) corresponds to the law of iterated expectation that $$ \Pr \{X_i=x\} = \sum_{k=1}^K p_U(k) \Pr \{X_i=x|U_i=u^k\}. $$ Given $\{p_{D,Z}(d,j)\}_{d,j}$, Equation (ref) can be written as a quadratic moment condition and Equation (ref) as a linear moment condition. I use these additional moments in orthogonalizing the moment $m$ so that the Neyman orthogonality holds.
Let $\tilde{\lambda}$ and $p$ denote vectorizations of $\left( \tilde{\Lambda}_0, \tilde{\Lambda}_1\right)$ and $\left( \{p_U(k)\}_k, \{p_{D,Z}(d,j)\}_{d,j}\right)$. The orthogonalized score is constructed with the additional moment function
$\phi$ simply collects the quadratic moments from (ref) across $(y,d,x,k)$, the linear moments from (ref) across $(d,x)$, and the linear moment $$ p_{D,Z} (d,j) = \mathbf{E}[ \mathbf{1}\{D_i=d, Z_i=z^j\}] $$ across $(d,j)$. To complete the orthogonalization, I show that the Jacobian matrix of $\phi$ has full rank.
Then, we can construct an additional nuisance parameter
and the orthogonalized score
satisfies the Neyman orthogonality. $\mu$ is estimated by taking a sample analogue of the expression above. Given estimators $\left( \hat{\tilde{\lambda}}, \hat{p}, \hat{\mu}\right)$, I estimate $\theta$ with $$ \binom{n}{2}^{-1} \sum_{i<i'} \psi \left( W_i, W_{i'}; \hat{\theta}, \hat{\tilde{\lambda}}, \hat{p}, \hat{\mu} \right) = 0. $$ $\widehat{F}_{Y(0),Y(1)}$ and $\widehat{F}_{Y(1)-Y(0)}$ denote the distributional treatment effect estimators we obtain from this two-step procedure.
Theorem (ref) establishes the consistency of the mixture weight estimators $\hat{\Lambda}_0$ and $\hat{\Lambda}_1$.
A direct corollary of Theorem (ref) is that $\widehat{\tilde{\Lambda}}_0, \widehat{\tilde{\Lambda}}_1$ are consistent for $\tilde{\Lambda}_0$ and $\tilde{\Lambda}_1$ at the rate of $\frac{1}{\sqrt{n}}$. Theorem (ref) establishes the asymptotic normality of the distributional treatment effect estimators.
The asymptotic variance is computed from a projection of the orthogonal scores:
In Sections (ref)-(ref), the standard error is obtained with a plug-in estimator for the asymptotic variance.
In this section, I discuss Monte Carlo simulation results. I generated $B=200$ random samples from DGPs with discrete $Y_i(1), Y_i(0), X_i, Z_i$ and $U_i$ where $M_Y = 3, M_X = 6, M_Z = 3$ and $K=3$: $Y_i \in \{1,2,3\}$, $X_i \in \{1,2,3,4,5,6\}$ and $Z_i \in \{1,2,3\}$.\footnote{The specifics of the DGPs are as follows: $p_U = (0.286, 0.286, 0.438)$,
and $\Lambda$s in the order of decreasing smallest singular value are
} The treatment $D_i$ was drawn randomly, independent of $Y_i(1), Y_i(0), X_i, Z_i$. In the first step nonnegative matrix factorization, I collapsed the support of $X_i$ so that the effective number of points in the support of $X_i$ is three. Thus, the conditional probability matrix $\mathbb{H}_0$ and $\mathbb{H}_1$ were $9 \times 3$ matrices. Across difference DGPs, I varied $\Lambda$, the conditional probability of $U_i$ given $Z_i$ which is shared across treated and untreated subpopulation, to vary the informativeness of the proxy variable $Z_i$ with regard to the latent variable $U_i$.
Table (ref) contains the bias and the root mean squared error (rMSE) of the distributional treatment effect estimators $\widehat{F}_{Y(1)-Y(0)}$. As $\Lambda$ becomes less informative about the distribution of $U_i$, i.e. the smallest singular value $\sigma_{\min}(\Lambda)$ decreases, the rMSE goes up. This suggests that the first step nonnegative matrix factorization estimation quality depends on how informative the proxy variables $X_i$ and $Z_i$ are for the latent variable $U_i$. Additionally, Table (ref) contains the coverage probability of the confidence interval constructed with the asymptotic standard error and the type $I$ error of the falsification test proposed in Subsection (ref). The 95% confidence interval shows mostly correct coverage, sometimes slightly too conservative, and the falsification test is valid.
In this section, we revisit JMR and estimate the distributional treatment effect of workplace wellness program on medical spending. Jointly with the Campus Well-being Services at the University of Illinois Urbana-Champaign, the authors of JMR conducted a large-scale randomized controlled trials. The experiment started in July 2016, by inviting 12,459 eligible university employees to participate in an online survey. Of 4,834 employees who completed the survey, 3,300 employees were randomly selected into treatment, being offered to participate in a workplace wellness program names iThrive. The participation itself was not enforced; the treated individuals were merely financially incentivized to participate by being offered monetary reward for completing each step of the wellness program. Thus, the main treatment effect parameter of JMR is the `intent-to-treat' effect. The workplace wellness program consisted of various activities such as chronic disease management, weight management, and etc. The treated individuals were offered to participate in the wellness program starting the fall semester of 2016, until the spring semester of 2018.
One of the main outcome variables that JMR studied is the monthly medical spending. Since the authors had access to the university-sponsored health insurance data, they had detailed information on the medical spending behaviors of the participants. Taking advantage of the randomness in assigning eligibility to the participants, JMR estimated the intent-to-treat type ATE of the workplace wellness program on the monthly medical spending. The ATE estimate on the first-year monthly medical spending, from August 2016 to July 2017, showed that the eligibility for the wellness program raised the monthly medical spending by \$10.8, with $p$-value of 0.937, finding no significant intent-to-treat effect.
In JMR, the authors acknowledge that the null effect on the mean does not necessarily mean null effect everywhere, though they themselves do not explore the treatment effect heterogeneity in the paper.\footnote{In the original dataset used in JMR, the authors had connected the medical spending variables to additional survey variables such as age, health behavior, salary, etc. They did not explore how the treatment effect interacts with the additional characteristics, but they did add these additional control variables through double Lasso. Adding the control variables increased the point estimate for the ATE (\$34.9) but the estimate still remained insignificant, with $p$-value being 0.859.} On page 1890, JMR state “there may exist subpopulations who did benefit from the intervention or who would have benefitted hard they partcipated.”\footnote{Damon Jones, David Molitor, and Julian Reif, “What do workplace wellness programs do? Evidence from the Illinois workplace wellness study,” The Quarterly Journal of Economics, vol. 134, no.4 (2019): 1747-1791.} I build onto this observation and estimate the distributional treatment effect of the randomly assigned eligibility for the wellness program. By looking at the distribution, I find the proportion of the subpopulation among treated population that benefitted from the treatment.
The dataset built by the authors of JMR fits the context of the short panel model in Example (ref). For each individual, the dataset contains monthly medical spending records for the following three time durations: July 2015-July 2016, August 2016-July 2017 and August 2017-January 2019. Since the experiment started in the summer of 2016 and the treated individuals were offered to participate in the wellness program starting the fall semester of 2016, the monthly medical spending record for July 2015-July 2016 could be thought of as a `pretreatment' outcome variable. Thus, we could use the information from the distribution of the pretreatment outcome variable to connect the treated subsample and the untreated subsample. The followings are the variables taken from the dataset.
In this specific empirical context, the common shock $V_{it}$ could be thought of as underlying health status and the treatment-status-specific shocks $\left( \varepsilon_{it}(1), \varepsilon_{it}(0) \right)$ could be thought of as additional random shocks such as susceptibility to the workplace wellness program or transient health shock which does not persist over time. The first-order Markovian assumption in Example (ref) is consistent with the health economics literature and broader economics literature of modeling household choices regarding health expenditure: grossman1972concept,wagstaff1993demand,jacobson2000family,yogo2016portfolio and more. Applying assumptions in Example (ref), the treatment is allowed to affect the underlying health status in the post-treatment period of August 2017-January 2019, but is assumed to be independent of the underlying health status in July 2015-July 2017.
Before applying the DTE estimators to the dataset, I implemented the falsification test with $K=5$.\footnote{When constructing $\mathbb{H}_0$ and $\mathbb{H}_1$ to be used in the first step nonnegative matrix factorization, we used the quintiles of the marginal distributions: $\left( -\infty, F_{Y}^{-1}(0.2), F_{Y}^{-1}(0.4), F_{Y}^{-1}(0.6), F_{Y}^{-1}(0.8), \infty \right)$ and so on. Thus, the matrices $\mathbb{H}_0$ and $\mathbb{H}_1$ were $25 \times 5$ matrices.} The test statistic is computed with a $25 \times 1$ vector
Theorem (ref) can be easily extended to the marginal distribution of $X_i$ as well and therefore we test the null (ref) with $$ T_n = n {W_n}^\intercal Avar(W)^{-1} W_n, $$ from $\sqrt{n} W_n$ being asymptotically normal. In the dataset, $T_n$ was 16.435 and its $p$-value was 0.901, passing the falsification test.
Figure (ref) contains the estimated joint distribution of the two potential outcomes from the nonnegative matrix factorization algorithm with $K=5$. For visibility, I first partitioned the potential outcome variable with quantiles $F_Y^{-1}(1/7), \cdots, F_Y^{-1}(6/7)$ and plotted the joint distribution of partitioned potential outcomes. Since the treated potential outcomes are plotted on the vertical axis, higher mass on the left-upper triangle means that the treatment reduces the medical spending. Overall, there is no definitive pattern. One notable observation is that the joint density is higher where $F_Y(Y_i(1)) \approx F_Y(Y_i(0)) \approx 0$ and $F_Y(Y_i(1)) \approx F_Y(Y_i(0)) \approx 1$. This is intuitive since on the two ends of the underlying health status spectrum, the effectiveness of the workplace wellness program must be limited.
Figure (ref) contains the estimated marginal distribution of the treatment effect and its 95% pointwise confidence interval. Note that the point estimates are mostly upward-sloping and lie between zero and one. Though the quadratic moment representation used in the DTE estimators does not impose any monotoncity or nonnegativity restrictions, the estimated marginal distribution violates these constraints only on a small subset of the range $[-1000,1000]$. Overall, it is unclear if more than half of the people would be better off from the treatment; the confidence interval for $\Pr \{Y_i(1) - Y_i(0) \geq 0\}$ contains 0.5, not being able to reject the null $\Pr \{Y_i(1) - Y_i(0) \geq 0\} \leq 0.5$.
As comparison, estimates for the upper bound and the lower bound from M1982,FP2010 are also provided in Figure (ref), as green dotted lines. The point estimates are consistent with the partial identification result, lying between the lower bound and the upper bound. The comparison highlights the gain of the point identification result, at the cost of assuming stronger identifying assumptions. For $\delta \in [-500,600]$, the 95% confidence interval is included in the partially identified set, giving us much bigger power in inference.
Lastly, the point identification helps us analyze the pattern of the treatment heterogeneity. Recall that the ATE estimate was inconclusive about the effectiveness of the treatment. However, the DTE estimates on $\Pr \{Y_i(1) - Y_i(0) \leq \delta\}$ for $\delta \leq -600$ and the DTE estimates on $\Pr \{Y_i(1) - Y_i(0) \leq \delta\}$ for $\delta \geq 400$ shows us interesting treatment effect heterogeneity patterns, in favor of implementing the treatment. The negative impact of the treatment, i.e. how much more money you spend under the treatment, is capped at \$400: $\widehat{F}_{Y(1)-Y(0)}(400) \approx 1$. On the other hand, the left tail of the treatment effect distribution is thicker, implying that some people are greatly benefitted from participating in the program: $\widehat{F}_{Y(1)-Y(0)}(-600) \approx 0.15$.
This paper presents an identification result for the joint distribution of treated potential outcome and untreated potential outcome, given conditionally random binary treatment. The key assumptions in the identification are that there exists a latent variable that captures the dependence between the two potential outcomes and that there exist two proxy variables for the latent variable. By assuming strict monotonicity for some functional of the conditional distribution of potential outcomes given the latent variable, I interpret the latent variable as `latent rank' and strict monotonicty as `latent rank invariance.' In implementation, I propose a first step nonnegative matrix factorization and a second step plug-in GMM. $\sqrt{n}$-consistency of the first-step estimator and the asymptotic normality of the second step GMM estimator are established. Lastly, I apply the estimation method to revisit JMR and find that the potential medical spendings are positively correlated at the two ends of the support and the marginal distribution of the treatment effect has thicker left tail.