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.
115,043 characters · 17 sections · 68 citation commands
TITLE
\setcounter{footnote}{0}
\raggedbottom
Randomization or unconfoundedness is probably the most important assumption for causal inference: see e.g., rubin1974estimating,rubin1978bayesian,lalonde1986evaluating,hahn1998role,hirano2003efficient,ix25lalonde. However, it can be too strong or even unrealistic in observational studies, in which case it is a common approach to rely on instrumental variables. But, it can be difficult to find credible instruments, and more importantly, they often require the researcher change the target population of interest, depending on the specific instruments available: see e.g., imbens1994late. Therefore, focusing on the case of a binary treatment, we develop a new methodology to learn about the average treatment effect (ATE) without relying on unconfoundedness or instrumental variables.
Our starting point is to note that the presence of an unobserved confounder is a key source of the problem. One possibility to proceed is to exploit a proxy variable for the unobserved confounder as in chalak2019identification. But, proxies are not always available, and therefore we focus on the case where we observe only the outcome and treatment, possibly along with exogenous covariates.
Following imbens2003sensitivity, we start with assuming that there exists an unobserved random variable $G$ such that the treatment assignment and potential outcomes become independent once we condition on $G$: there may be additional observed covariates that need to be controlled for, but we suppress them from our discussion for simplicity. So, $G$ represents an unobserved confounder. Our key assumption is then that $G$ has finite support, where the number $\bar G$ of the mass points of $G$ is unknown: as $\bar G$ increases, the population exhibits a higher degree of unobserved heterogeneity. Therefore, our work builds upon the sensitivity analysis of imbens2003sensitivity, as well as modern techniques that exploit discretization of unobserved heterogeneity as in e.g., bonhomme2022discretizing.
Sensitivity analyses for causal inference in this context have been discussed by several authors: e.g., rosenbaum2002sensitivity,imbens2003sensitivity,masten2018identification,yadlowsky2022bounds,bonvini2022sensitivity. The ideas are all similar: i.e., the distribution of the potential outcomes and treatment assignment is allowed to deviate from complete independence, and the identified set for ATE is traced out as a function of the degree of the deviation. The last four references are particularly relevant for us.
imbens2003sensitivity and yadlowsky2022bounds formulate sensitivity analysis through explicit models of unobserved confounders. By contrast, masten2018identification and bonvini2022sensitivity remain agnostic about the underlying confounding mechanism. Specifically, masten2018identification bound the discrepancy between the marginal probability of treatment and its conditional counterpart given the potential outcomes. Meanwhile, bonvini2022sensitivity posit that unconfoundedness may fail for a subset of the population and treat the proportion of such individuals as a sensitivity parameter.
masten2018identification give suggestions for how to interpret their sensitivity parameter, and by following the view advocated in rosenbaum2002rejoinder, they also argue that interpretations of sensitivity parameters are not always necessary. In this view, economic interpretability or specific reasons for deviations from the baseline model are not essential for sensitivity analysis; what matters is simply that the baseline model may be false for some reason. We take a different stance. We believe that it is important to understand what a particular deviation represents. A researcher may want to assess how plausible a deviation is, or how seriously she should take potential deviations in a given empirical setting. Such reasoning is facilitated when deviations from the baseline model have a clear economic interpretation. Moreover, an interpretable framework enables the researcher to gauge the relative likelihood of different deviations and, consequently, to “weight” the resulting causal inferences accordingly. Other authors have also advocated the importance and usefulness of interpretability of sensitivity analysis. For instance, heckman2000causal wrote, “The bounding and sensitivity analysis movement is likely to be more influential if it relies on explicit economic models and uses economically interpretable models to conduct semiparametric bounding and sensitivity analyses.” See also cinelli2020making from the statistics community.
We therefore prefer an explicit approach on the presence of unobserved confounders, as it makes it easier to provide an economic interpretation of deviations from the baseline model. However, both of the two explicit frameworks we mentioned earlier have their own limitations. imbens2003sensitivity adopts a fully parametric approach, where the complete joint distribution of the potential outcomes, treatment assignment, and unobserved confounder needs to be specified. This can be restrictive. yadlowsky2022bounds's approach is nonparametric but is less transparent since it restricts the odds-ratio of the treatment assignment given the unobserved confounder. Therefore, our goal is to develop an explicit and easy-to-interpret framework that is less restrictive than imbens2003sensitivity.
Similarly to imbens2003sensitivity, we begin by assuming Gaussianity for the potential outcome distribution, but our setup is more flexible. Specifically, we assume that the potential outcome $Y(d)$ is Gaussian with mean $\mu_g(d)$ conditional on $G=g$, where $d\in\{0,1\}$ indicates the treatment status, and $G$ is the only source of potential confounding. The variable $G$ is assumed to be discrete, but the number $\bar G$ of its support points is unknown. To ensure that the group label $g$ represents a well-defined type, we assume that $\bigl( \mu_g(0), \mu_g(1) \bigr) \neq \bigl( \mu_{\tilde g}(0), \mu_{\tilde g}(1) \bigr)$ whenever $g \neq \tilde g$, where the ordering of the group label $g$ is normalized by the lexicographic ordering of the mean vectors of the potential outcomes; this ordering remains well-defined even if $G$ is multi-dimensional, provided that its elements are all discrete. For notational convenience, we label the support of $G$ as $\{1,2,\cdots, \bar G\}$. Finally, we leave the dependence between $G$ and the treatment assignment unspecified.
In this formulation, heterogeneity in treatment effects arises from two components: mean heterogeneity across latent types $G$ and an idiosyncratic residual component, although only the former is relevant for ATE. While this structure is shared with imbens2003sensitivity, we do not impose functional-form restrictions on the mean heterogeneity. In particular, the location component of $Y(d)$ given $G$, i.e., $\sum_{g=1}^{\bar G}\mathds{1}(G=g)\mu_g(d)$, is unrestricted apart from its finite support. As a result, the Gaussianity assumption on the residual component imposes relatively weak restrictions on the marginal distribution of $Y(d)$. Indeed, finite Gaussian mixtures are dense in a broad class of distributions, so they can approximate a wide range of distributions with arbitrary accuracy as $\bar G$ becomes sufficiently large.
It is an old idea to decompose heterogeneity in treatment effects into two components. For example, imbens2003sensitivity incorporates this through a fully parametric structure for sensivity analysis, while abadie2024instrumental use a similar decomposition to improve the asymptotic mean squared error of linear instrumental variable estimators. We adopt this decomposition for causal inference in the absence of instrumental variables. An important consequence of our formulation is that the resulting version of the “Manski bounds” on ATE remains finite even when the outcome variable has unbounded support.
We note that gardner2020identification also exploits discrete unobserved confounders for causal inference. However, gardner2020identification implicitly assumes that the number of latent types in the treatment group always coincides with that in the control group, and that the researcher knows how the types revealed in one group correspond to those revealed in the other; see Footnote (ref) for details. As a result, gardner2020identification point-identifies ATE regardless of how large $\bar G$ is. We view these assumptions as too strong. Instead, we derive sharp identified sets for ATE as a function of a candidate value of $\bar G$ in its identified set. We then characterize how these sets vary with $\bar G$ and compare them with the Manski bounds.
The exact economic interpretation of $G$ depends on the context, but we can generally view $\bar G$, the number of mass points of $G$, as the number of “unobserved types” in the population that may affect both the treatment assignment and the potential outcomes: e.g., different skill levels or various industries. In essence, $\bar G$ is a key parameter that describes the unknown degree of heterogeneity in the population. Our sensitivity analysis then involves expressing the sharp identified set for ATE as a function of $\bar G$. However, we do not need to trace out all integers for a complete sensitivity analysis. Indeed, we show that $\bar G$ is partially identified, and therefore there are only finitely many values of $\bar G$ that are consistent with the data. Further, we show that there is a special cutoff value in the identified set of $\bar G$ such that whenever $\bar G$ is larger than the cutoff, the sharp identified set for ATE becomes a version of the Manski bounds manski2003partial,manski2010partial: we say that it is a version of the Manski bounds because we do not assume that the support of the outcome is bounded, but our bounds have the same structure as the Manksi bounds except that they use the means of the mixture component distributions of the outcome given the treatment status.
The fact that the unobserved confounder $G$ has finitely many mass points leads to finite mixture models of the observed outcome given the treatment status. Finite mixtures have been frequently used in econometrics: e.g., ichimura1998maximum,arci03fin,kasahara2009nonparametric,hen14par,compiani2016using. In our framework, however, they serve as a tool for causal inference. To be more specific, let $D$ denote the observed treatment status, and let $Y = D Y(1) + (1-D) Y(0)$ be the observed outcome. Then, the distribution of $Y$ given $D=d$ is a finite Gaussian mixture. Consequently, the data generate two finite mixture models, one for each treatment status, whose parameters can be identified using standard methods. The Gaussianity assumption is introduced solely to guarantee identification of the mixture parameters, although similar identification results can be obtained under alternative assumptions: see e.g., yak68id,mp00finmix. Indeed, we show how the baseline model can be extended to incorporate exogenous covariates, account for sample selection, and accommodate non-Gaussian distributions of $Y(d)$ given $G$. The combination of distributional assumptions and discrete unobserved heterogeneity has also been employed by bonhomme2022discretizing in the context of panel data models. Our contribution is to exploit this combination for causal inference and sensitivity analysis.
One of the challenges for causal inference here is that different types in the population do not always correspond to distinct component distributions of $Y$ given $D=d$. Therefore, neither the number $\check G(1)$ of mixing components in the treatment group, nor that in the control group denoted by $\check G(0)$, generally reveals the true value of $\bar G$. Nevertheless, we show that the sharp identified set for $\bar G$ can be expressed as a function of $\check G(1)$ and $\check G(0)$, both of which can be estimated at an arbitrarily fast rate: see e.g., chen09order.
A complete sensitivity analysis can then be conducted by tracing out the sharp identified set for ATE for each value of $\bar G$ in its identified set. The number of admissible values of $\bar G$ increases quadratically in general as $\check G(1)$ and $\check G(0)$ increases. However, our identification results show that there exists a cutoff value $\bar G_C$ such that $\bar G_C$ increases linearly as $\check G(1)$ and $\check G(0)$ increase, and that the sharp identified set for ATE corresponds to (a version of) the Manski bounds whenever $\bar G$ is larger than $\bar G_C$. Therefore, the total number of admissible values of $\bar G$ that need to be checked for a complete sensitivity analysis is finite and increases only linearly as $\check G(0)$ and $\check G(1)$ increase.
For the purpose of estimation, we take a plug-in approach. Specifically, we first estimate the mixture orders of the observed outcome for the treatment and control groups, for which we suggest using the method of chen09order. We then estimate the mixture parameters by MLE, replacing the unknown mixture orders with their first-step estimates. Finally, for each admissible value of $\bar G$, we plug these estimators into the expression for the identified set for ATE. In the supplement, we establish the asymptotic properties of the proposed estimators, including $\sqrt{n}$-consistency, and develops valid inference procedures.
The proposed method provides applied researchers with a tool for program evaluation that does not depend on either the unconfoundedness assumption nor the availability of good instruments: see, e.g., the surveys by iw09survey,ix25lalonde. For instance, even when the treatment group comes from experimental data and the control group from observational data, our approach continues to provide a valid method for learning about ATE, although we refer interested readers to yang2025cross for a more systematic approach to integrating experimental and observational data for causal inference. In addition, our approach is not tied to particular research designs such as difference-in-differences or regression discontinuity designs, nor does it require the design-specific identification assumptions underlying those approaches, such as parallel trends or exogenous running variables.
Finally, we illustrate the practical value of our approach by applying it to the widely studied dataset of lalonde1986evaluating. The experimental control group in this dataset contains many observations with zero earnings. Excluding these observations raises concerns about selection bias. To address this issue, we apply both our baseline model and its modified version for selection. The results are not overly different. The benchmark Heckit regression results, which can be interpreted as the case of $\bar G=1$, show that the coefficient of the training indicator is approximately 0.07. However, when we allow the presence of an unobserved confounder, the admissible set of $\bar G$ is $\{2,3,4\}$. Among these values, $\bar G\geq 3$ produces the Manski bounds, whereas $\bar G= 2$ yields sharp identified $\{ 0.28, 0.56 \}$ for the average effect of job training.
The rest of the paper is organized as follows. Section (ref) describes the setup, and introduces the identifiable mixture parameters and group clustering issues. Section (ref) presents our identification results, characterizing how the identified set of ATE changes as a function of $\bar G$. Section (ref) describes the estimation procedure. Section (ref) develops extensions to allow for covariates, selection, and non-Gaussian distributions. Section (ref) provides an empirical illustration and, finally, Section (ref) concludes. All proofs are collected in the Appendix. In addition, we provide a Supplement that includes the proofs of the auxiliary lemmas used in the Appendix and the asymptotic properties of the estimator together with the inference procedures.
Before we proceed, we make brief comments on our notation.
An array with its elements separated by commas such as $v = (v_1,\cdots, v_m)$ is always considered a column vector, as well as any element of $\mathbb{R}^m$. Given two vectors $v \in \mathbb{R}^{m_1}$ and $w \in \mathbb{R}^{m_2}$, the array $( v , w)$ is considered a (column) vector in $\mathbb{R}^{m_1 + m_2}$. For any subset of $\mathbb{R}^m$, the lexicographical order is our default order relation, and we use the symbol $\prec$ to denote it. For any subset of matrices, its elements are lexicographically ordered by using the row-wise vectorization. Given an ordered set $S = \{ s_1 \prec \cdots \prec s_m \}$, we write $( v_s )_{s \in S} = ( v_{s_1}, \cdots , v_{s_m} )$. The super-script $^{\mathpalette\raiseT{\intercal}}$ means transpose. In addition, we say that a tuple $(S_1, S_2, \cdots, S_k)$ of subsets of $S$ is an ordered partition of a set $S$ if it satisfies: (i) $S_i \neq \emptyset$ for all $i=1,\dots,k$; (ii) $S_i \cap S_j = \emptyset$ for all $i \neq j$; and (iii) $\bigcup_{i=1}^k S_i = S$.
We employ the usual notation $N ( a , b)$ to represent a normal distribution with mean $a \in \mathbb{R}$ and variance $b >0$, while $N ( a , 0 )$ denotes a degenerate distribution that assigns probability one to the value $a$. With a slight abuse of notation, for $A \in \mathbb{R}^m$ and a positive semi-definite matrix $B \in \mathbb{R}^{m \times m}$, we write $Y \sim N ( A , B )$ to mean that $v^{\mathpalette\raiseT{\intercal}} Y \sim N ( v^{\mathpalette\raiseT{\intercal}} A , v^{\mathpalette\raiseT{\intercal}} B v )$ for all $v \in \mathbb{R}^m$. The symbols $\stackrel{p}{\rightarrow}$ and $\stackrel{d}{\rightarrow}$ denote convergence in probability and distribution, respectively, while w.p.a.1 abbreviates “with probability approaching one.”
In this section, we show how causal effects can be identified when treatment assignment may be correlated with potential outcomes, but only through discrete confounding factors. Specifically, we assume that there are only finitely many unobserved innate types, where the number of types is unknown. Further, we assume that the potential outcomes of interest follow Gaussian mixtures, where the number of mixing components, which is unknown as well, is constrained by that of the innate types. Gaussianity is not essential for our discussion, but it is convenient because identifiability of Gaussian mixtures is well understood yak68id. An extension to a selection model is discussed in Section (ref), while extensions to other non-Gaussian mixture models are presented in Section (ref). The treatment assignment can be correlated with the unobserved types, and therefore the types are unobserved confounding factors. The causal parameter of interest is the average treatment effect (ATE).
Let $D \in \{ 0 , 1\}$ be a binary treatment, and let $Y(d)$ for $d\in \{0,1\}$ be potential outcomes. The observed outcome is $Y = DY(1) + (1-D) Y(0)$. For the purpose of identification analysis, we assume that the joint distribution of $(Y,D)$ is known. We assume that there are finitely many innate types, which are represented by $G$. The econometrician does not observe $G$: all she knows about it is that it is discrete with finite support: i.e., $G$ takes a value from $\{1,2,\cdots, \bar G\}$, where $\bar G \geq 2$ is unknown. There may be observable covariates $X$, but we suppress them throughout the discussion for expositional simplicity. Therefore, all statements should be understood as conditional on $X=x$ when such covariates are present. We also discuss an alternative approach to incorporating covariates in Section (ref), which is more restrictive but more pragmatic for implementation.
We now assume that the potential outcomes can be modeled by Gaussian mixtures.
Assumption (ref) implies that $Y(d)$ is independent of $D$ once we condition on $G$. Therefore, $G$ plays the role of the unobserved confounder. A similar idea has been employed in imbens2003sensitivity for a sensitivity analysis. However, Assumption (ref) provides a more flexible environment than the fully parametric approach of imbens2003sensitivity. All that is assumed for the distribution of $G$ is that $\bar G$ is finite. Assumption (ref) is silent about dependence between $D$ and $G$: the only restriction on their joint distribution is that all types and treatment statuses are realized with positive probability. Also, the conditional mean of $Y(d)$ given $G$ is nearly unrestricted: $\sum_{g=1}^{\bar G}\mathds{1}(G=g)\mu_g(d)$ has a multinomial distribution with an unknown number of mass points, where the multinomial distribution can approximate any distribution arbitrarily well if $\bar G$ is sufficiently large. Finally, it is worth noting that $\bar G$ often admits an economic interpretation, although the precise interpretation depends on the context. This feature facilitates sensitivity analyses in the spirit of imbens2003sensitivity and may help motivate additional restrictions that the researcher may want to impose on $\bar G$, or more generally on the distribution of $G$.
We emphasize a few important features of Assumption (ref). Assumption (ref) does not require that all $\mu_g(d)$'s be distinct for a given $d\in\{0,1\}$, albeit that $\bm{\mu}_1,\cdots, \bm{\mu}_{\bar G}$ must be distinct as vectors. In other words, the mean vectors are all distinct across different types, but that does not necessarily imply that all types are revealed by the marginals. In fact, Assumption (ref) allows $Y(1)$ and $Y(0)$ to have different numbers of mixing locations, where the number of distinct values in $\mu_1(d),\cdots, \mu_{\bar G}(d)$ for a given $d\in \{0,1\}$ can be strictly smaller than $\bar G$. Assumption (ref) imposes lexicographic ordering on the mean vectors such that
This is purely about how we label the innate types $1,2,\cdots, \bar G$, and there is no loss of generality. Finally, Assumption (ref) focuses on Gaussian location mixtures for the distribution of $Y(d)$, but this imposes only few restrictions. Specifically, as we commented earlier, the mixture mean $\sum_{g=1}^{\bar G}\mathds{1}(G=g)\mu_g(d)$ can approximate any distribution arbitrarily well as long as $\bar G$ is sufficiently large. Then, the distribution of $Y(d)$ is obtained by adding an independent Gaussian noise to it, where the Gaussian noise is homoskedastic across different types. The homoskedastic error across types is to rule out the possibility that two distinct types share the same mean for both $d=0$ and $d=1$.
We present a simple example of the data-generating process (DGP) described in Assumption (ref).
For a given $d\in\{0,1\}$, not all of $\mu_1(d),\cdots, \mu_{\bar G}(d)$ need to be distinct. Let $\check{G}(d)$ be the number of distinct elements among $\mu_1(d),\cdots, \mu_{\bar G}(d)$: note that Assumption (ref) allows either $\check{G}( 0) = 1$ or $\check{G}( 1) = 1$, but not both at the same time. We will denote the distinct elements by $\check{\mu}_1(d) < \cdots < \check{\mu}_{\check G (d)}(d)$. The relationship between $\mu_g(d)$'s and $\check{\mu}_g(d)$'s determines an ordered partition on $\{1,2,\cdots, \bar G\}$ of the form
noting that $ \bigcup_{j =1}^{\check{G}(d)} \mathcal{G}_j (d) = \{ 1 , \cdots , \bar{G} \} $. Essentially, each $\mathcal{G}_j(d)$ is a cluster of the types that correspond to the distinct value $\check{\mu}_j(d)$ and the resulting order of the tuple $ ( \mathcal{G}_{1} (d) , \cdots , \mathcal{G}_{\check{G}(d)} (d) )$ is given by such distinct values. We remark that, in contrast to an unordered partition, the order in which the subsets $\mathcal{G}_j (d) $ are listed matters.
An important implication of Assumption (ref) is that, for both $d=0$ and $d=1$, there exist ordered partitions satisfying Equation (ref). In particular, every type $g$ must be represented in both the treatment and control groups. If the treatment and control groups were drawn from different populations, then comparing them for causal inference would be simply nonsensical.
The ordered partitions of the innate types satisfy several properties as a consequence of the normalization imposed by the lexicographic ordering of the distinct mean vectors $\bm{\mu}_1,\cdots, \bm{\mu}_{\bar G}$. First, we always have $\max\ \mathcal{G}_j(0) < \min\ \mathcal{G}_k(0)$ whenever $j<k$. In words, this means that $\mathcal{P} (0)$ partitions $\{1 , \cdots , \bar{G}\}$ into $\check{G}( 0 )$ consecutive blocks. However, this is not necessarily the case for the treatment group: i.e., for $\mathcal{P} ( 1 ) $, $j<k$ does not necessarily imply that $\max \mathcal{G}_{j} (1) < \min \mathcal{G}_{k} (1)$, nor even $\max \mathcal{G}_{j} (1) \leq \min \mathcal{G}_{k} (1)$. Finally, we remark that $\mathcal{G}_j(0)\cap \mathcal{G}_k(1)$ must be either empty or a singleton for any $j, k$.
We continue to consider Example (ref) to illustrate the ideas.
Define $\pi_{j} (d , d^\prime) := {\mathbb{P}}\{ G \in \mathcal{G}_{j} (d) \mid D = d^\prime \}$ for $( d , d^\prime ) \in \{ 0 , 1 \}^2$. The distribution of the observed outcome, conditional on $D$, will follow a Gaussian mixture model such that
By Proposition 2 in yak68id, the parameters that are directly identified from the Gaussian mixtures are
Moreover, ${\mathbb{P}} (D = 1)$ is trivially identified. Therefore, we will assume that the parameters in (ref) and (ref) as well as the treatment probability ${\mathbb{P}}(D=1)$ are known when we discuss identification of the average treatment effect. In other words, the researcher knows that there are ordered partitions $\mathcal{P} (d) = ( \mathcal{G}_{1} (d) , \cdots , \mathcal{G}_{\check{G}(d)} (d) )$ of the population, $d=0,1$, where each $\mathcal{G}_j(d)$ corresponds to the mixture component mean $\check{\mu}_j(d)$, but the exact compositions of $\mathcal{G}_j(d)$'s are unidentified, and thus the researcher does not know which of the mixture components a specific type $g$ belongs to.\footnote{gardner2020identification considers a similar setup and approach with a focus on panel data. However, his analysis implicitly relies on strong assumptions. In particular, it assumes that the types revealed by mixtures in the treatment and control groups are exactly matched. For illustration, suppose that the mixture model of the treatment group reveals two distinct types, say $\{1,2\}$. gardner2020identification then implicitly assumes that the mixture model of the control group also reveals exactly two distinct types, say $\{a,b\}$, and that the econometrician knows that $\{1\}$ and $\{a\}$ represent the same type in the population, and likewise for $\{2\}$ and $\{b\}$: i.e., it is implicitly assumed that $\check G(0) = \check G(1)$, and more importantly, the types in the population that are revealed by $\check\mu_j(1)$ and $\check\mu_j(0)$ always coincide, from which it is guaranteed that $\mu_j(d) = \check\mu_j(d)$ for both $d\in \{0,1\}$. Consequently, gardner2020identification obtains point identification of the average treatment effect regardless of the number of mass points of the unobserved confounder. This “exact matching” assumption appears to be too restrictive and arbitrary in realistic settings.}
Using Example (ref), the situation is illustrated below.
Unlike the parameters in (ref) and (ref), the clustering partitions $\mathcal{P} (0 ) $ and $\mathcal{P} (1 ) $ are unidentified. Without additional assumptions, we do not know, e.g., whether any $\mathcal{G}_{j} (0)$ and $\mathcal{G}_{k} (1)$ share any elements (types) in common. Formally, this means that the composition of each set $ \mathcal{G}_{j} (d) $ is unknown to the researcher. Furthermore, $ \pi_{j} (0 , 1 )$ and $ \pi_{j} (1 , 0 )$ are unidentified as well. However, the total number of innate types $\bar G$ is partially identified, as we show in Lemma (ref) below.
To facilitate the subsequent identification analysis, we show in the next subsection that there is a succinct way of representing the admissible partitions of innate types, i.e., the ordered partitions that are consistent with Assumption (ref). Specifically, we will show that, under Assumption (ref), each potential candidate for the latent pair of ordered partitions $\bigl( \mathcal{P} (0), \mathcal{P} (1) \bigr)$ can be uniquely represented by a $\check{G}(0)\times \check{G}(1)$ matrix.
For $(m_1 , m_2) \in \mathbb{N}^2$, let $\{ 0 , 1\}^{m_1\times m_2}$ denote the set of $m_1\times m_2$ binary matrices, where each entry is either zero or one.
Condition (1) mandates that each row of $\bm{\iota} \in \mathscr{I}$ contains at least one non-zero entry, while (2) requires the same for each column. The set $\mathscr{I}$ is clearly nonempty, and it is a singleton if and only if either $\check{G}(0) = 1$ or $\check{G}(1) = 1$, in which case the only element of $\mathscr{I}$ is a matrix with ones in all its entries.
The idea behind Definition (ref) is as follows. Let $( \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}p}} ( 0 ) , \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}p}} ( 1 ) )$ be an admissible candidate for the (true) latent pair of underlying partitions $( \mathcal{P}(0) , \mathcal{P} ( 1 ) )$, where $ \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}p}} ( 0 ) = ( \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{1} ( 0 ) , \cdots , \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{ \check{G}(0) } ( 0 ) )$ and $ \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}p}} ( 1) = ( \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{1} ( 1 ) , \cdots , \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{\check{G}(1)} ( 1 ) )$. Note first that the intersections $\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_j(0)\cap \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_k(1)$ for $j\in\{1,\cdots, \check{G}(0)\}$ and $k\in \{ 1,\cdots, \check{G}(1)\}$ are all that matters to completely characterize the candidate $( \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}p}} ( 0 ) , \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}p}} ( 1 ) )$ because $\cup_{j=1}^{\check{G}(0)} \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_j(0)\cap \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_k(1) = \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_k(1)$ and $\cup_{k=1}^{\check{G}(1)} \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_j(0)\cap \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_k(1) = \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_j(0)$. However, under Assumption (ref), $\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_j(0)\cap \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_k(1)$ is always either empty or a singleton, and Conditions (1) and (2) in Definition (ref) trace the non-emptiness patterns in those intersections. Specifically, the matrix $\bm{\iota} \in \{ 0 , 1\}^{\check{G}(0) \times \check{G}(1) }$ associated with the candidate $ ( \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}p}} ( 0 ) , \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}p}} ( 1 ) )$ is constructed such that $\iota_{j,k}$ equals one if and only if $\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_j(0)\cap \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_k(1) $ is nonempty.
The following lemma formalizes the preceding discussion.
Lemma (ref) follows as a direct consequence of the lexicographic ordering and the distinct mean vector condition specified in Assumption (ref). A specific closed-form expression for the bijection is provided in Section (ref) in the Appendix. Due to the existence of such a bijection, hereafter, we will use $\bm{\iota} \in \mathscr{I}$ to refer to a candidate for the (true) latent pair of ordered partitions $( \mathcal{P}(0) , \mathcal{P} ( 1 ) )$.
By using Lemma (ref), we can show that $\bar G$, the number of innate types, is partially identified.
Lemma (ref) shows that $\bar G$ is point identified if and only if $\min\{ \check{G}(0), \check{G}(1) \}$ is equal to $1$. In general, $\bar G$ belongs to a set of finitely many integers. The sharp identified set for $\bar G$ is obtained by noting that the number of types associated with the partition represented by $\bm{\iota} = ( \iota_{j,k} ) \in \mathscr{I}$ is given by the sum of its entries, i.e.\ $\sum_j \sum_k \iota_{j,k}$. Hence, the sharp identified set for $\bar G$ becomes $\{ \sum_j \sum_k \iota_{j,k} : \bm{\iota}= ( \iota_{j,k} ) \in \mathscr{I} \}$. Intuitively, the lower bound $\bar G_L$ is trivial, and it corresponds to the case of no hidden type in either the treatment or the control group: the total number of innate types cannot be smaller than the number of distinct types in either the treatment or control group. The upper bound is a consequence of the normalization of the lexicographic ordering, and it corresponds to the case of maximal clustering. Specifically, given $\check{G}(0)$ and $\check{G}(1)$, we know that $\check{\mu}_1(0)<\dots<\check{\mu}_{\check{G}(0)}(0)$ and $\check{\mu}_1(1)<\dots<\check{\mu}_{\check{G}(1)}(1)$. Therefore, the largest number of different ways of lexicographically pairing them is obtained when all the means of the treatment group are paired with each of the means in the control group.
Before proceeding to the next section, we revisit Example (ref) to illustrate the ideas discussed above. In particular, we emphasize that it is necessary to consider ordered partitions to derive the bijection established in Lemma (ref).
Our parameter of interest is ATE. Since
where ${\mathbb{P}} \{ G \in \mathcal{G}_{j} (d) \} = \pi_{j} (d , 0 ) {\mathbb{P}} ( D = 0) + \pi_{j} (d , 1 ) {\mathbb{P}} ( D = 1)$ by the law of total probability, we can express ATE as follows:
where $\{\pi_k (1 , 0 ): k=1,\cdots,\check{G}(1)\}$ and $\{ \pi_j ( 0 , 1 ): j=1,\cdots, \check{G}(0)\}$ are unidentified conditional probabilities. Assumption (ref) does not directly impose any restrictions on these unidentified parameters except that they are non-degenerate probabilities such that $\sum_{j=1}^{\check{G}(1)}\pi_j(1,0) = \sum_{j=1}^{\check{G}(0)} \pi_j(0,1) = 1$. Therefore, it is not too difficult to obtain the sharp identifiable set of $\tau$ by searching over the right simplexes. Specifically, let $\Pi_j = \{ \textbf{p} \in \mathbb{R}_{++}^j : \sum_{l=1}^j p_l= 1 \}$ for $j \in \mathbb{N}$, and define
Theorem (ref) is straightforward and easy to understand: estimating $\mathcal{T}$ does not involve any optimization beyond estimating the identified mixture parameters. It is worth noting that $\mathcal{T}$ is reminiscent of the Manski bounds manski2003partial,manski2010partial: indeed,
whereas, for $d=0,1$, we have
for which we recall that either $\check{G}(0)=1$ or $\check{G}(1)=1$ is allowed. The fact that we have open intervals in the case of $\check{G}(d) > 1$ is because all of $\check{\mu}_j(d)$'s for $d=0,1$ and $j=1,2,\cdots, \check{G}(d)$ have non-zero weights under Assumption (ref). Since the support of $Y(d)$ for $d=0,1$ is unrestricted under Assumption (ref), the na\"{i}ve Manski bounds would be unbounded. However, Assumption (ref) does impose finite mixture structures on $Y(d)$, which are exploited in Theorem (ref). Specifically, the distributions of $Y$ given $D=d$ for $d=0,1$ are also finite mixtures, and Theorem (ref) formally establishes sharp bounds based on the mixing components that are identified from the distribution of $Y$ given $D=d$.
For the sake of illustration, we have computed $\tau$ and $\mathcal{T}$ in the setup of Example (ref) below.
Theorem (ref) takes a completely agnostic stance about $\bar G$, which can be extreme. For instance, the researcher may want to conduct a sensitivity analysis for various values of $\bar G$, or she may be willing to impose further restrictions on $\bar G$. Theorem (ref) does not provide a useful path for these purposes. $\bar G$ is a key parameter that summarizes how much heterogeneity there is in the population, and the researcher knows that it is an integer between $\bar G_L$ and $\bar G_U$ by Lemma (ref). As $\bar G$ moves away from $\bar G_L$, the situation becomes more challenging in that the number of types that are not revealed in either the treatment or the control group becomes larger. The researcher may be willing to rule out extreme possibilities of “too many clusterings,” and to focus on the case where the total number of innate types in the population is revealed by either the treatment or the control group. Alternatively, the researcher may consider using a weighted-averaging scheme for different possibilities for the value of $\bar G$. Theorem (ref) is not convenient to open up those possibilities. Therefore, in the following subsection, we will characterize the sharp identified set for ATE as a function of $\bar G$.
In order to describe the sharp identified set for $\tau$ for a given value of $\bar G$, we first characterize the identified set for a given $\bm{\iota}\in\mathscr{I}$. Here, $\bm{\iota}\in\mathscr{I}$ represents a pair of specific partitions on $\{1,2,\cdots, \bar G\}$ that are allowed under Assumption (ref), whereas $\bar G$ is an economic parameter given by nature. We start with introducing the following set of matrices.
The sets $\mathcal{M}_{\bm{\iota}} (0)$ and $\mathcal{M}_{\bm{\iota}} (1)$ are collections of nonnegative matrices, where the nonnegative entries are in accordance with those in $\bm{\iota}\in\mathscr{I}$, and their columns and rows are restricted to sum to $\pi_j(0,0)$ and $\pi_k(1,1)$, respectively. The purpose of these sets is to provide a systematic way to track all possible values that the conditional probability ${\mathbb{P}}\{ G \in \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{\bm{\iota}, j} (0) \cap \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{\bm{\iota}, k} (1) \mid D = d \}$ can take for $d\in \{0,1\}$ and $(j , k ) \in \{ 1, \cdots , \check{G}(0) \} \times \{ 1 , \cdots , \check{G}(1) \}$, where $\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}p}}_{\bm{\iota}}(d) := \bigl( \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{\bm{\iota}, 1}(d),\cdots, \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{\bm{\iota}, \check{G}(d)} (d) \bigr )$ denotes the ordered partition associated with $\bm{\iota}$ in group $d$.
To be specific, for a given matrix $\text{$\mathbf{p}$} = ( p_{j,k} ) \in \mathcal{M}_{\bm{\iota}} (0)$, the entry $p_{j,k}$ represents a possible value of ${\mathbb{P}} \{ G \in \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{ \bm\iota , j} (0) \cap \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{\bm\iota , k} (1) \mid D = 0 \}$ if $( \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}p}}_{\bm{\iota}} (0) , \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}p}}_{\bm{\iota}} (1) )$ were the underlying partition. The restrictions on the column-wise and row-wise sums are because $\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}p}}_{\bm{\iota}}(d)$ is a partition on $\{1,\cdots,\bar G\}$: therefore, taking the union of $\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{ \bm\iota , j} (0) \cap \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{\bm\iota , k} (1)$ across $j$ (or $k$) leads to $\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{\bm\iota , k} (1)$ (or $\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{\bm\iota , j} (0)$). For example, Condition (b) of Definition (ref).(1) ensures that the values of ${\mathbb{P}}\{ G \in \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{\bm{\iota}, j} (0) \cap \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{\bm{\iota}, k} (1) \mid D = 0 \}$, $k =1 , \dots , \check{G}(1)$ are logically consistent with that of the identified conditional probability $\pi_{j} ( 0 , 0) $ as
We also note that the zero entries of the matrix $\mathbf{p}$ are associated with the case where $\operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{\bm{\iota},j} (0) \cap \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{\bm{\iota}, k} (1) $ is empty. In addition, we trivially have that $\sum_j \sum_k p_{j,k} = 1$ due to Condition (b). A symmetric argument can be applied to the components of $\mathbf{q}\in \mathcal{M}_{\bm{\iota}}(1)$, which are associated with possible values of ${\mathbb{P}}\{ G \in \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{\bm{\iota},j}(0) \cap \operatorname{\text{\usefont{U}{BOONDOX-cal}{m}{n}g}}_{\bm{\iota}, k} (1) \mid D = 1 \}$.
For $\bm{\iota}\in\mathscr{I}$, define
and note that $\mathcal{T}_{\bm{\iota}}$, which is a subset of $\mathcal{T}$, is always an interval because $\mathcal{M}_{\bm{\iota}} (d)$'s are convex. Let $\boldsymbol{\iota}_{\mathcal{P}} \in \mathscr{I}$ represent the unknown latent pair of ordered partitions $( \mathcal{P} (0) , \mathcal{P} (1) )$, i.e., the true pair.
An immediate consequence of Theorem (ref) and Lemma (ref) is that $\cup_{\bm{\iota} \in \mathscr{I}} \mathcal{T}_{\bm{\iota}} = \mathcal{T}$. In fact, we can divide the set $\mathcal{T}$ by using potential values of the economic parameter $\bar G$ instead of $\bm{\iota}$. We will discuss this issue in detail later in this section.
Before we proceed, we revisit Example (ref) to illustrate the ideas.
All the $\mathcal{T}_{\bm{\iota}}$'s described in Example (ref) are either a singleton or an open interval. It can be shown that this is always the case in general: see Lemma (ref) in the Appendix. Specifically, $\mathcal{T}_{\bm{\iota}}$ is a singleton if and only if $\sum_{j,k} \iota_{j,k} = \check{G} (0) = \check{G} (1)$. We also note that Lemma (ref) in the Appendix provides closed-form expressions for the infimum and supremum of $\mathcal{T}_{\bm{\iota}}$. Therefore, there is no need to solve an optimization problem to compute $\mathcal{T}_{\bm{\iota}}$.
We now consider grouping $\bm{\iota}\in \mathscr{I}$ by using the potential values of $\bar G$, i.e., the number of innate types. As we argued earlier, $\bar G$ is an economic parameter summarizing the amount of unobserved heterogeneity in the population, whereas $\bm{\iota}\in\mathscr{I}$ is a consequence of Assumption (ref) given $\bar G$.
For an integer $m$ satisfying $ \bar G_L \leq m \leq \bar G_U$, we define the sets
The next theorem establishes that $ \mathcal{T} ( \bar{G} ) $ would be the sharp identified set for $\tau$ if $\bar{G}$ were known, and it also characterizes the shape of $ \mathcal{T} ( m ) $ for different values of $m$. Define $\bar G_C := \check{G}(0) + \check{G}(1) - 1$.
Theorem (ref) immediately implies that
Also, the cut-off value $\bar G_C$ is noteworthy. For all $m\geq \bar G_C$, $\mathcal{T}(m)$ is essentially $\mathcal{T}$: they are equal, or the difference is at most finite when $\bar G_C \leq m < 2\bar G_C - 2$. The case where $\mathcal{T}\backslash \mathcal{T}(m)$ is a nonempty finite set should be considered a technical nuisance that arises only because $\mathcal{T}(m)$ may be a union of disjoint open intervals, whereas $\mathcal{T}$ is always an open interval. Even in this case, it is guaranteed that $\mathcal{T}$ and $\mathcal{T}(m)$ share the same closure so that $\inf\mathcal{T} = \inf\mathcal{T}(m)$ and $\sup\mathcal{T} = \sup\mathcal{T}(m)$ for all $m\geq \bar G_C$. For all $m < \bar G_C$, $\mathcal{T}(m)$ is a strictly smaller set than $\mathcal{T}$, and $\mathcal{T}(m)$ strictly shrinks as $m$ becomes smaller. Therefore, if we restrict $\bar G$ such that $\bar G<\bar G_C$, then the sharp identified set for $\tau$ will be strictly smaller than $\mathcal{T}$.
The cutoff value $\bar G_C$ naturally arises as the largest value allowed for $\bar G$ in some situation. One such case is when there is rank invariance in the two groups: i.e., for any $g, \tilde g\in \{1,2,\cdots, \bar G\}$, we have
In other words, the treatment does not work against any particular type: if the high-talent type does better than the low-talent type on average without the treatment, then the high type does better on average with the treatment as well, and vice versa. We note that the rank invariance does not rule out the possibility that the treatment may be more beneficial on average to one type than to another: it is the relative rankings that do not change. It can be shown that if the rank invariance condition holds, then the sharp identified set for $\bar G$ becomes $\bigl\{m\in\mathbb{N}: \bar G_L \leq m \leq \bar G_C \bigr\}$.\footnote{To see this, note that under rank invariance, for both $d = 0 ,1$, $\mathcal{P} ( d ) $ partitions $\{ 1 , \cdots , \bar{G} \}$ into $\check{G}(d)$ consecutive blocks such that $\max\ \mathcal{G}_j(d) < \min\ \mathcal{G}_k(d)$ whenever $j<k$. This immediately rules out $\bar{G} >\bar{G}_C$ because $\mathcal{G}_j(0)\cap\mathcal{G}_k(1)$ is either empty or a singleton for all $j,k$ by Assumption (ref). However, $\bar{G} = \bar{G}_C$ can be achieved, e.g., by choosing $\iota$ so that all entries in the first row and the last column are equal to one, and all remaining entries are zero.} Therefore, $\bar G = \bar G_C$ corresponds to the worst-case scenario under the rank invariance.
The causal effect $\tau$ is not point-identified in general because of potential confounding from unobserved heterogeneity. However, Theorem (ref) suggests that there are several conservative measures of interest. For example, the researcher may present $\ell_m := \inf \mathcal{T}(m)$ for all integers $m$ between $\bar G_L$ and $\bar G_C$, which will provide a complete picture on the lower bounds on the causal effect $\tau$. We may do the same for the upper bounds if aggressive measures are desired. If more succinct summaries are wanted, then we may consider
where $w_m$'s are weights chosen by the researcher. Observe that $\ell_{\bar{G}_C} $ represents the Minimum Lower Bounds (MLB) on $\tau$, whereas $\tau_{ALB}$ shows the Average Lower Bounds (ALB). For example, \[ \tilde{\tau}_{ALB} := \frac{1}{\bar G_C - \bar G_L + 1}\sum_{m=\bar G_L}^{\bar G_C} \ell_m \] is the Maximum Entropy ALB for $\tau$ in that it uses the uniform weights on the sharp identified set of $\bar G$ that reflects the rank invariance condition in (ref). For the use of the maximum entropy principle under partial identification in different contexts, see e.g., jun2024information. We can similarly define the Maximum Upper bounds (MUB) and the Average Upper Bounds (AUB).
The MLB and ALB have their own advantages and disadvantages. For instance, the MLB on $\tau$ is easier to interpret, and it does not require the researcher choose a prior on $\bar G$. But it may be overly conservative. The ALB, which is the average of the lower bounds on $\tau$ for different values of $\bar G$, is always less conservative than MLB since it reflects the researcher's attitude about uncertainty on $\bar G$. We do not advocate a particular stance here.
Before proceeding, we provide the value of each $\mathcal{T}(m)$ in Example (ref) as an illustration.
To conclude this section, we highlight that additional restrictions on the partitions can be systematically incorporated into our framework to construct the sharp identified set for $\tau$. For instance, if we invoke rank invariance, the sharp identified set becomes \[ \mathcal{T}_{\mathrm{ri}} : = \bigcup_{\bm{\iota} \in \mathscr{I}_{\mathrm{ri}}} \mathcal{T}_{\bm{\iota}} , \] where $ \mathscr{I}_{\mathrm{ri}} \subset \mathscr{I}$ denotes the subset of admissible partitions under Assumption (ref) and rank invariance. Specifically, it can be shown that $\mathscr{I}_{\mathrm{ri}} = \{ \bm{\iota} \in \mathscr{I}: \ \iota_{j^\prime m^\prime} = 0 \ \text{if} \ \iota_{j m} = 1 \ \text{for some} \ j < j^\prime , m < m^\prime \}$ from inspecting the proof of Lemma (ref). In addition, for a given $m \in \{ \bar G_L, \cdots , \bar G_C \}$, we can construct the sharp identified set for $\tau$ under the rank invariance as follows: \[ \mathcal{T}_{\mathrm{ri}} (m): = \bigcup_{\bm{\iota} \in \mathscr{I}_{\mathrm{ri}} ( m ) } \mathcal{T}_{\bm{\iota}} \quad \text{with} \ \ \mathscr{I}_{\mathrm{ri}} ( m ) = \mathscr{I} ( m ) \cap \mathscr{I}_{\mathrm{ri}}. \]
We remark that if $\bar G=\check G(0)=\check G(1)$, then $\mathscr{I}(\bar G)\cap \mathscr{I}_{\mathrm{ri}}$ contains only the identity matrix. Consequently, $\mathcal{T}_{\mathrm{ri}}(\bar G)$ is a singleton. This observation helps explain the point-identification result of gardner2020identification. Reformulating the discussion in Footnote (ref), his framework implicitly imposes both rank invariance and the restriction $\bar G=\check G(0)=\check G(1)$, which together imply point identification of $\tau$.
While extending Theorem (ref) to the sets $\mathcal{T}_{\mathrm{ri}}$ and $\mathcal{T}_{\mathrm{ri}}(m)$ is beyond the scope of this paper, the following example provides an illustrative special case.
In this section, we propose an estimator of $\ell:=( \ell_{\bar G_L}, \cdots, \ell_{\bar G_C} )$, while we provide its asymptotic properties in the Supplement together with the inference procedures. We do not discuss the cases where $m>\bar G_C$, but there is no loss of generality here because Theorem (ref) implies $\ell_{\bar G_C} = \ell_m$ whenever $\bar G_C \leq m \leq \bar G_U$. Also, we focus on the lower bounds $\ell_m$, but the methodology and results can be symmetrically applied to the upper bounds $\sup \mathcal{T}(m)$ as well.
We consider a random sample $\{ (Y_1, D_1), \dots, (Y_n, D_n) \}$ drawn from the distribution of $(Y, D)$ under Assumption (ref). Estimating $\ell$ requires preliminary estimators of the identified parameters, for which we present the following two-step procedure.
In the first step, for each $d=0,1$, we estimate the mixture order $ \check G(d)$ by using any method that produces an estimator $\hat G(d) \in \mathbb{N}$ such that $\hat G(d) = \check{G}(d)$ w.p.a.1, i.e.,
For instance, one can use chen09order's estimator, which involves the use of penalized maximum likelihood with a cluster-based penalty function that follows fan01var. In the second step, for each $d=0,1$, we estimate the vector of the mixture parameters \[ \gamma(d) := \bigl( \pi_1(d,d),\cdots, \pi_{\check{G}(d)-1}(d,d), \check{\mu}_1(d),\cdots, \check{\mu}_{\check{G}(d)}, \sigma^2(d)\bigr) \] by maximum likelihood, using $\hat{G}(d)$ in lieu of $\check{G}(d)$, and we denote the resulting estimator by \[ \hat \gamma(d) := \bigl( \hat\pi_1 ( d , d ),\cdots, \hat\pi_{\hat{G}(d)-1}(d,d) , \hat{\mu}_1(d) , \cdots , \hat{\mu}_{\hat G(d)} (d) , \hat\sigma^2 (d) \bigr) , \] emphasizing that $\hat{\mu}_{j} (d)$ is an estimator of $\check{\mu}_{j} (d)$, not ${\mu}_{j} (d)$.\footnote{When $\check{G}(d)=1$, $\gamma(d)$ reduces to $ ( {\mu}_1(d) , \sigma^2(d) )$, and its estimator $\hat\gamma(d)$ coincides with the usual maximum likelihood estimator under normality.} To obtain $\hat \gamma(d)$ and ensure its existence w.p.a.1, we remark that the maximization is carried out over a compact subset $\Gamma ( \hat{G} (d) ) \subset \mathbb{R}^{2 \hat{G} (d)} $ specified in the Supplement. We refer to mp00finmix for further discussion on implementation aspects.
In order to construct the estimator of $\ell$ and to derive its asymptotic properties in the next subsection, we introduce the following notation. Let $\Lambda := \bigl( \mathbb{P}(D=0) , \gamma(0) , \gamma(1) \bigr)$ denote the vector of parameters characterizing the joint distribution of $(Y,D)$ and write $\mathscr{I}_{LC} : = \cup_{m=\bar{G}_L}^{\bar{G}_C} \mathscr{I} (m)$. Then, define the vector-valued function $\mathcal{L}(\Lambda) := \bigl( \mathcal{L}_{\bm\iota} (\Lambda) \bigr)_{\bm\iota \in \mathscr{I}_{LC}}$ by
where $\bar{j}_k (\bm{\iota}) := \max \bigl\{ j \in \{ 1, \cdots , \check{G} (0) \} : \iota_{jk} = 1 \bigr\}$ and $ \ushort{k}_j (\bm{\iota}) := \min \bigl\{ k \in \{ 1, \cdots , \check{G} (1) \} : \iota_{jk} = 1 \bigr\}$. Moreover, for a given real vector $v = (v_{\bm\iota} )_{\bm\iota \in \mathscr{I}_{LC}}$ indexed by $\mathscr{I}_{LC}$, define the function $\mathcal{M}$ by
From Lemma (ref) in the Appendix, we can now write
Now, we consider a plug-in approach to estimate $\ell$ by using (ref) and the preliminary estimators, i.e., $(\hat G(0), \hat G(1))$ and $(\hat\gamma (0), \hat\gamma (1))$. Let $\hat\Lambda := \bigl( \hat{\mathbb{P}}(D=0) , \hat\gamma(0), \hat\gamma(1) \bigr)$, where $\hat{{\mathbb{P}}} ( D = d) : = n^{-1} \sum_{i=1}^n \mathds{1} ( D_i = d)$. We then define the estimator of $\ell$ by
where $\hat{\mathscr{I}}_{LC}$, $\hat{\mathcal{L}}_{\bm\iota}$, and $\hat{\mathcal{M}}$ are defined analogously to $\mathscr{I}_{LC}$, $\mathcal{L}_{\bm\iota}$, and $\mathcal{M}$, respectively, with $(\hat G(0) ,\hat G(1))$ replacing $(\check G(0) ,\check G(1))$. In this regard, note that $\hat{\mathscr{I}}_{LC} = \mathscr{I}_{LC}$, $\hat{\mathcal{L}} = \mathcal{L}$ and $\hat{\mathcal{M}} = \mathcal{M}$ hold w.p.a.1 due to Eq.\ (ref). We remark that computing $\hat{\mathcal{L}}(\hat{\Lambda})$ can be demanding when $\hat{G}(0)$ and $\hat{G}(1)$ are large, since the cardinality of $\hat{\mathscr{I}}_{LC}$ increases rapidly with these quantities. In such cases, one may instead employ approximation methods such as resampling elements from $\hat{\mathscr{I}}_{LC}$. The resulting approximation error can be made statistically negligible by increasing the number of replications.
The function $\mathcal{M}_m$ generally requires that we consider all entries $\mathcal{I}(m)$ except the one corresponding to the MLB $\ell_{\bar{G}_C}$. To see this point, let $\ushort{\bm\iota}$ denote the $\check G(0) \times \check G(1) $ matrix that has all ones in the first column and in the last row, and zeroes everywhere else, noting that $\ushort{\bm\iota} \in \mathscr{I} ( \bar{G}_{C} )$. Theorems (ref) and (ref), as well as Lemma (ref) in the Appendix, imply
As a result, letting $\hat{G}_{C} = \hat G(0) + \hat G(1) - 1$, $\hat{\ell}_{\hat{G}_{C}}$ can be computed simply by replacing the parameters in (ref) with their corresponding estimators. Of course, this is equivalent to applying the plug-in principle to the formula of the Manski bounds in (ref).
We remark that ${\ell}_{\bar{G}_{C}}$ can also be estimated via an overfitted MLE, i.e., by specifying a mixture order larger than $\check{G}(d)$ rather than estimating it in a first step. To see this, observe that
An alternative plug-in estimator can therefore be constructed by replacing ${\mathbb{P}} (D =d)$ and $\mathbb{E} (Y\mid D =d)$ with their sample analogs, and $\max_{j} \check\mu_j (0)$ and $\min_j \check\mu_j (1)$ with the corresponding estimates obtained from the overfitted MLE.
Because the overfitted MLE yields consistent estimates of $\max_j \check\mu_j (0)$ and $\min_j \check\mu_j (1)$ -- a consequence of consistency in Wasserstein distance -- this alternative estimator is also consistent, though converging at a slower rate of $n^{-1/4}$ up to a polylogarithmic factor ho16strong. Thus, we do not consider this approach in the asymptotic analysis provided in the Supplement. In practice, however, if the primary objective is to obtain a consistent estimate of the Manski bound ${\ell}_{\bar{G}_{C}}$, this alternative avoids the need to estimate the mixture orders.
We have been implicit about covariates $X$, and therefore, our discussion so far can be understood as conditional on $X=x$: i.e., the total number $\bar G$ of types is allowed to be heterogeneous across different subpopulations, and our results should be understood as bounding the conditional ATE (CATE). However, when it comes to estimation and inference, this approach of completely localizing at (or around) $X=x$ may not be practically attractive, especially when the dimension of $X$ is high. Further, censoring or selection is an issue that frequently arises in applications, and it generally makes na\"{i}ve Gaussian models unrealistic. Therefore, we now discuss extensions to address covariates and selection; we will focus on the simplest case of Type I censoring, but extensions to other types of selection are straightforward.
Let $X$ be a random vector of covariates, which does not include a constant, and consider the shifted potential outcome
where $\beta(d)$ is a vector of slope coefficients. In this setting, $Y(d) := \sum_{g=1}^{\bar G} \mathds{1}(G=g) Y_g (d)$ represents the potential residual after partialling out the effects of the covariates, where $Y_g (d)$ can be interpreted as a group-specific “residual.”
We consider the following assumption about a possible selection mechanism and the residual distributions. Let $\tilde{Y} : = D \tilde{Y}(1) + (1-D) \tilde{Y}(0) $ denote the shifted outcome, and let $\lambda$ be the inverse Mills ratio.
By condition (1), $X$ cannot include a constant. Conditions (1) and (2) guarantee identification of $(\eta_0 , \eta_1)$ via probit, while (3) ensures identification of $\beta(d)$ and $\sigma_{\epsilon\nu}(d)$ via the least squares of $Y$ on $X$ and $\lambda(\eta_0 + X'\eta_1)$ conditional on $D=d$. So, hereafter, $\bigl( \eta_0, \eta_1, \beta(d), \sigma_{\epsilon\nu}(d) \bigr)$ will be treated as known parameters for identification purposes.
In this setup, CATE can still be represented by $\tau$ as defined in equation (ref): i.e., \[ \text{CATE}(x) = x'\{ \beta(1) - \beta(0)\} + {\mathbb{E}}\{ Y(1) - Y(0) \} = x'\{ \beta(1) - \beta(0)\} + \tau. \] Hence, constructing identified sets for CATE reduces to constructing identified sets for $\tau$. In the special case of $\sigma_{\epsilon\nu}(d) = 0$, selection is exogenous, and it is innocuous to focus on “uncensored observations.” Then, we can obtain identified sets for $\tau$ by applying our previous results to the homogenized outcome $Y = \tilde Y - X'\beta(d)$ given $D=d$. More generally, however, selection necessitates further modification, which we elaborate on below.
The potential residuals $(Y(0), Y(1))$ satisfy Assumption (ref), and therefore the conditional distribution of $Y = D Y(1) + (1-D)Y(0) = \tilde Y - X'\beta(d)$ given $D=d$ (but unconditional on $S=1$) remains a Gaussian mixture as specified in equations (ref)-(ref). Therefore, Theorem (ref) would still apply if the number $\check G(d)$ of mixing components, and \[ \gamma(d) = \bigl( \pi_1(d,d),\cdots, \pi_{\check{G}(d)-1}(d,d), \check{\mu}_1(d),\cdots, \check{\mu}_{\check{G}(d)}, \sigma^2(d)\bigr) \] were known. Therefore, all we need to establish here is that these parameters are still identified when we condition on $S=1$. In the following lemma, we show that conditioning on $S=1$ yields a non-Gaussian mixture likelihood, from which we can identify and estimate $\check G(d)$ and $\gamma(d)$. Let $\Phi$ denote the standard normal distribution function, and let $\phi := \Phi'$ be its density.
From this lemma, estimating the identified sets for $\tau$ is just a matter of estimating the mixture parameters $\check G(d)$ and $ \gamma(d)$; see Section (ref). With this aim, considering a random sample
we suggest two possibilities.
The first one consists in a three-step procedure. In the first two steps, we estimate $\eta_0, \eta_1, \beta(d)$, and $\sigma_{\epsilon\nu}(d)$ by a combination of probit and least squares, often referred to as Heckit. In the third, we first estimate $\check G(d)$, for which we suggest using chen09order's penalization method. Then, we estimate $\gamma(d)$ by maximum likelihood using the selected sample and treating $\check G(d)$ as the true mixture order. Specifically, this procedure can be described as follows.
The second approach consists in estimating all model parameters jointly via (conditional) maximum likelihood, except for the mixture orders $(\check{G}(0), \check{G}(1))$. Specifically, we first estimate these mixture orders (e.g., by chen09order's method) and then, treating such estimates as the true mixture orders, we jointly estimate all other parameters via (conditional) maximum likelihood. In this regard, note that the joint log-density of $(S \tilde{Y} , S D , X , S)$ evaluated at $(s \tilde{y} ,s d,x,s)$ is given by
while the closed-form expression for $f_{\tilde Y\mid D,X,S}$ can be obtained from Lemma (ref).
The three-step procedure is attractive for its computational simplicity, while the alternative joint estimation approach is asymptotically efficient and attains the efficiency bound. Details on asymptotic theory and inference are provided in the Supplement.
Since we use finite mixtures to deal with unobserved discrete types, we start by introducing an identifiable mixture class. We suppress covariates again for simplicity. Consider a parametric family of univariate distributions $\mathscr{F} := \{ F(\cdot ; \theta):\ \theta \in \Theta \}$, where $\Theta \subseteq \mathbb{R}^K$, for some $K\in \mathbb{N}$, and $F$ satisfies $\int |y| \ dF(y;\theta) <\infty$ for all $\theta \in \Theta$. Let $\mathscr{M}$ be the class of finite mixtures based on $\mathscr{F}$ such that
where $\Pi_m = \{ \textbf{p} \in \mathbb{R}_{++}^m : \sum_{l=1}^m p_l= 1 \}$ and $\Omega_m\subseteq \Theta^m$.
The class $\mathscr{M}$ depends on $\mathscr{F}$ as well as $\{\Pi_m:\ m\in\mathbb{N}\}$ and $\{\Omega_m: m\in\mathbb{N}\}$, but we suppress this dependence for the sake of notational simplicity. The parameter space for the class $\mathscr{M}$ is $\cup_{m=1}^\infty \Pi_m\times \Omega_m$. This formulation does not require that the true number of mixing components be known: e.g., if $m_0$ is the true number of component distributions, then the true parameters $(p_1,\cdots, p_{m_0})$ and $(\theta_1,\cdots, \theta_{m_0})$ are in $\Pi_{m_0}$ and $\Omega_{m_0}$, respectively. Also, it is worth noting that $\{\Omega_m:\ m\in\mathbb{N}\}$ need not be the same as $\{\Theta^m:\ m\in\mathbb{N}\}$. For example, certain elements of $\theta_m$ may remain constant across $m$, as in the Gaussian location mixture models discussed earlier, in which the variance does not vary across components.
Following mb88mix, we say that $\mathscr{M}$ is an identifiable class when the following condition is satisfied: if
where $\theta_j \neq \theta_{j'}$ and $\tilde\theta_j \neq \tilde\theta_{j'}$ for all $j\neq j'$, then it follows that $\tilde m = m$ and $\bigl( ( \tilde{p}_1 , \tilde{\theta}_1 ) , \cdots , ( \tilde{p}_{\tilde{m}} , \tilde{\theta}_{\tilde{m}} ) \bigr)$ is a permutation of $\bigl( ( p_1 , \theta_1 ) , \cdots , ( p_m , \theta_m ) \bigr)$. There are many well-known examples of identifiable mixture classes. For instance, the class of Poisson mixtures is identifiable; see kar01rob and chen09order. Likewise, exponential mixture models are also identifiable. Additional examples and discussions can be found in eve13fin and mp00finmix.
Now, we discuss extending Assumption (ref) beyond Gaussian mixtures. Let $Y_g(d)$ represent the potential outcome for the treatment status $d\in \{0,1\}$ for an individual who belongs to type $g\in \{1,2,\cdots, \bar G\}$, where $\bar G$ is unknown to the researcher. We assume that the marginal distribution of $Y_g(d)$ is represented by $F\bigl(\cdot; \theta_g(d) \bigr)\in \mathscr{F}$, and denote the mean vector of $\bigl( Y_g(0), Y_g(1)\bigr)$ by $\bm{\mu}_g : = \bigl( \mu_g(0), \mu_g(1) \bigr) : = \bigl( \mathbb{E}\{ Y_g(0)\}, \mathbb{E}\{ Y_g(1)\} \bigr)$, where $\mu_g(d)$ depends on $\theta_g(d)$ because of $\mu_g(d) = \int y\ dF\bigl(y;\theta_g(d) \bigr)$. We now state a generalization of Assumption (ref).
Assumption (ref) can be easily nested into Assumption (ref) by setting
with $\epsilon(d) \sim N(0, \sigma^2(d) )$ being independent of $(D,G)$ so that $\theta_g (d) = ( \mu_g (d), \sigma^2(d) )$. Condition ((ref)) is simple regularity: there are only finitely many types that are relevant. Condition ((ref)) is an innocuous normalization, provided that the mean vectors are all distinct across different types: i.e., $\bigl( \mu_g(0), \mu_g(1) \bigr) \neq \bigl( \mu_{g'}(0), \mu_{g'}(1) \bigr)$ whenever $g\neq g'$. Condition ((ref)) allows $Y(d)$ to depend on $G$, but $Y(d)$ becomes independent of $D$ once $G$ is controlled for. Therefore, $G$ is the only source of confounding. It is worth noting that condition ((ref)) does not impose any restrictions on the dependence between $G$ and $D$, or between $Y_g(d)$ and $Y_{g'}(d')$ as long as $g\neq g'$. Condition ((ref)) requires that the class $\mathscr{M}$ be identifiable, but, it does not require that all $\theta_g(d)$'s be distinct for different values of $g$ given a fixed value of $d=0,1$. For example, if $\theta_1(0) = \theta_2(0)$, then the number of mixing components that is identified from the control group will be strictly smaller than $\bar G$. In fact, neither the treatment group nor the control group may reveal the number $\bar G$ of types in the population. Condition ((ref)) ensures that different component distributions have distinct means, meaning that we only focus on type heterogeneity in the mean.
As we commented above, for any given $d\in \{0,1\}$, not all of $\theta_g(d)$'s need to be distinct for different values of $g$. However, conditions ((ref)) and ((ref)) impose some restrictions. For instance, if $\theta_g(0) = \theta_{g'}(0)$, then we must have $\theta_g(1) \neq \theta_{g'}(1)$. Otherwise, we would have $\bigl( \mu_g(0), \mu_g(1) \bigr) = \bigl( \mu_{g'}(0), \mu_{g'}(1)\bigr)$, in which case there is no strict lexicographic ordering between the means of types $g$ and $g'$, and therefore, in view of condition ((ref)), $g$ and $g'$ should simply be treated as the same type.
Let $\check{\theta}_1(d),\cdots, \check{\theta}_{\check{G}(d)}(d)$ be distinct parameter values among $\theta_1(d),\cdots, \theta_{\bar G}(d)$. Let $\check{\mu}_j(d)$ be the mean corresponding to $\check{\theta}_j(d)$ so that $\check{\mu}_1(d) < \cdots < \check{\mu}_{\check{G}(d)}(d)$, and define \[ \mathcal{G}_j(d) := \bigl\{ g\in \{1,\cdots, \bar G\}:\ \check{\theta}_j(d) = \theta_g(d) \bigr\} = \bigl\{ g\in \{1,\cdots, \bar G\}:\ \check{\mu}_j(d) = \mu_g(d) \bigr\}, \] where the second equality is by condition ((ref)) in Assumption (ref). We define \[ \pi_j(d,d') := \mathbb{P}\{ G \in \mathcal{G}_j(d)\mid D = d'\} \] for $(d,d')\in \{0,1\}^2$ and $j = 1,\cdots, \check{G}(d)$, and we obtain the following lemma. Let $Y = DY(1) + (1-D) Y(0)$ be the observed outcome as before.
As in Section (ref), we suggest first estimating $\check G(d)$ by any procedure that yields a consistent estimator in the sense of equation ((ref)). For instance, one may apply chen09order's estimator, or makha21's penalized-likelihood method when the mixing distribution is multidimensional.
Given the estimates of the mixture orders, the remaining mixture parameters $(\pi_j(d,d), \check{\theta}_j(d))$ can then be estimated by MLE, replacing the true mixture orders with their estimated counterparts, as in Section (ref). We refer to the Supplement for a detailed discussion on estimation and inference for general mixture models.
We now illustrate our methodology by using the data from lalonde1986evaluating and imbens2003sensitivity. We take two approaches: one is to ignore zero earnings to use the Gaussian mixture method of Section (ref), and the other is to explicitly address the censoring issue by using the model we presented in Section (ref). For this exercise, we focus on the experimental treatment and control groups of the data, both of which were drawn from the same population with relatively poor labor market prospects such as ex-drug addicts, ex-criminal offenders, and high-school dropouts. Using the PSID control as a comparison group does not seem to be reasonable because we do not believe that the condition in (ref) holds even after controlling for covariates such as age, age-squared, and years of education.
The treatment $D$ indicates whether the agent attended a training program or not. The (pre-homogenized) outcome $\tilde Y(d)$ is the logarithm of the earnings in 1978, and the convariates we use for homogenization consist of age, age-squared, and the years of education; these are the traditional variables we use in a Mincer equation. We control for the same covariates for the selection equation based on the Probit.
The (experimental) control group contains a total of 260 observations, which is reduced to 168 after removing observations with zero earnings. The treatment group has a total of 185 observations, and the number becomes 140 after removing observations with zero earnings. Therefore, we have a total of 308 observations to estimate the treatment effect of the job training on the log-earnings, where we control for age, age-squared, and years of education.
Given the limited sample size, we impose the restrictions that $\beta(0) = \beta(1)$ and $\sigma_{\epsilon\nu}(0) = \sigma_{\epsilon\nu}(1)$: i.e., the coefficients of the interactions of the treatment with the covariates and the inverse Mills ratio are all equal to zero. As a benchmark, we first estimate the treatment effect by a simple linear model with and without selection. Estimation results by simple ordinary least squares (OLS) and by two-step Heckit are presented in the following table.
All the coefficients (except for the intercept) are individually insignificant at a 10% significance level; the standard errors are not reported then. However, the point estimates still appear to be reasonable for a Mincer equation. Therefore, we rely on them for the purpose of homogenization. According to these estimates, receiving the job training increases earnings by about 7%.
We now estimate the mixture models for both the treatment and control groups following the three-step procedure from the previous section. For this purpose, we first apply chen09order's method with and without selection to obtain $\hat G(0) = \hat G(1) = 2$. We use the minimum suggested value for the tuning parameter (i.e., the smallest penalization for increasing the value of $\hat{G}(d)$), which is suggested in chen09order's Monte Carlo section. Despite that, our estimates $\hat G(d)$ are small. Before we proceed, we present the estimates of the identified parameters from the mixture models.
The estimates reported in Table (ref) can be used to show how the sharp bounds for the training effect will change as the number of unobserved types changes. Since $\hat G(0) = \hat G(1) = 2$, we have $\bar G\in \{2,3,4\}$. However, it suffices to check $\bar G \in \{2,3\}$, because the Manski bounds will be the best we can obtain whenever $\bar G\geq \check{G}(0) + \check{G}(1) - 1 = 3$. The sharp identified sets for the training effect, under $\bar G = 2$ and $\bar{G} \in \{3,4\}$, are presented in the following table.
Tables (ref) and (ref) show a complete picture of the effect of the job training and its robustness. The coefficients of “training” are the estimates of the effect of the job training with and without addressing selection when we assume that there is no unobserved heterogeneity, which can be interpreted as $\bar G = 1$. In contrast, fitting the mixture models suggests that $\bar G$ may be as large as $4$. If $\bar G \in \{3,4\}$, then we already hit the Manski bounds, in which case we do not learn much about the training effect. If $\bar G = 2$, then the regression estimates in Table (ref) may be substantially underestimating the training effect. For instance, under rank invariance, the estimated effect is $0.439$ in the no-selection model and $0.564$ when selection is accounted for.
We have shown how to learn about the average treatment effect without assuming unconfoundedness by assuming that unobserved confounders are discrete with the number of mass points unknown. As concluding remarks, we make two comments.
First, it is important for our methodology that both the treatment and control groups are influenced by every latent type in the population. If the treatment and control groups are sampled from the same population, then this requirement is generally plausible. However, it may become problematic if the researcher combines data from different sources. For example, in lalonde1986evaluating, the experimental treatment group and the PSID control group are drawn from markedly different populations, making it difficult to justify the assumption that they share the same set of latent types.
Second, we have focused on the average treatment effect. Extending our methodology to other causal parameters, such as quantile treatment effects, appears to be a nontrivial question. This difficulty is partly due to the absence of a law of iterated quantiles analogous to the law of iterated expectations. More recently, masten2025relaxing extended the c-dependence framework of masten2018identification to sensitivity analysis for a broader class of causal parameters, including quantile treatment effects. Nevertheless, our approach remains unique in that it formulates sensitivity analysis directly through restrictions on latent confounders. Consequently, any substantive prior knowledge or assumptions about the confounders can be incorporated in a transparent and systematic manner.