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.
59,836 characters · 12 sections · 50 citation commands
Lee bounds for random objects
\address[D. Kurisu]{Center for Spatial Information Science, The University of Tokyo, 5-1-5, Kashiwanoha, Kashiwa-shi, Chiba 277-8568, Japan.} \email{[email removed]}
\address[Y. Okamoto]{Graduate School of Economics, Kyoto University, Yoshida Honmachi, Sakyo, Kyoto 606-8501, Japan. } \email{[email removed]}
\address[T. Otsu]{Department of Economics, London School of Economics, Houghton Street, London, WC2A 2AE, UK. } \email{[email removed]}
Randomized controlled trials are among the most credible approaches for answering causal questions in economics, medical science, and many other disciplines. However, their validity can be undermined by endogenous sample attrition or non-response, which are ubiquitous in empirical studies. A common approach to address this issue, developed in a seminal work by Lee:2009, is to bound the average treatment effect of units whose outcomes are observed regardless of the treatment status. Recent empirical studies employing the Lee bounds include Alfonsi_etal:2020, Atal_etal:2024, and Dobbie_Fryer:2015, among many others.
The methodology of the Lee bounds has been extended in several directions to deal with noncompliance issues Chen_Flores:2015, continuous treatments Lee_liu:2025, multi-layered selection Kroft_etal:2025, and relaxed monotonicity assumptions as well as multiple outcomes Semenova:2025, for example. However, almost all existing works focus on outcomes that lie in a Euclidean space, even though many economically important outcomes do not live in Euclidean spaces.
One prominent class is compositional data, where components are nonnegative and sum to one. For example, labor and macro economists study intra-household time allocation (e.g., Cardia_Gomme:2018, Lise_Yamada:2019), while public economists analyze the composition of government expenditures (e.g., Brender_Drazen:2013, Shelton:2007). Other relevant examples include the composition of household consumption, energy mix across sources, agricultural land use, and portfolio shares in finance.
Another important class consists of distribution-valued outcomes. Economists have studied income and wage distributions (e.g., Battistin_etal:2009, DiNardo_etal:1996) as well as price dispersion (e.g., Cavallo:2017, Pennerstorfer_etal:2020). Related examples include the distributions of daily step counts, daily television viewing time, or hourly electricity consumption, which may be of interest in health, education, and public economics.
Beyond these classes, researchers also encounter network outcomes, such as cross-country trade and travel flows, as well as transaction networks among firms. Interval-valued data are likewise ubiquitous in econometric analysis---for example, when outcomes are reported in brackets for confidentiality, elicited as ranges to improve response rates, or observed with rounding error. Finally, correlation matrices are central objects for studying connectivity, including applications such as product recommendation and brain functional connectivity. We refer to Dubey:2024 and references therein for a general review of statistical analysis of such random objects, along with a range of applications.
Despite their prevalence, standard implementations of Lee bounds are not directly applicable to these cases, since such outcomes are not situated in Euclidean spaces. To overcome this limitation, we extend the scope of partial identification analysis via Lee bounds by adapting the methodology of metric statistics for random objects. In particular, we define the expected potential outcomes for random objects by introducing the notion of the Fr\'echet mean, a direct generalization of the conventional mean toward a general metric space,\footnote{Perhaps a popular example of the Fr\'echet mean in economics is the Aumann mean for random sets (e.g., BeMo08), which is shown to be a special case of the Fr\'echet mean (kuri:25), and our methodology also applies to random sets (see, Example (ref) below).} and then establish partial identification results on the expected potential outcome by exploiting a representation in an embedded Euclidean or Hilbert space, where the bounding argument in Lee:2009 may be adapted. Intuitively, our identified set is obtained by pulling back the Lee-type bounds on the embedded Euclidean or Hilbertian variables, and this approach enables analysis that takes into account the geometric structure of the metric space in which the outcome takes values. Based on our partial identification strategy, the treatment effect can be characterized by a difference in the embedded space, a difference in the projections, or geodesics in the metric space. Our estimator for the identified set is obtained by taking the sample analogs, and we present a valid bootstrap procedure for the confidence region of our identified set.
Causal inference on random object outcomes has been increasingly popular in recent econometrics and statistics literature, e.g., Gunsilius:2023, kurisu_etal:24, KuZhOtMu25a, KuZhOtMu25b, and ZhKuOtMu25. The present paper contributes to this literature by extending the method of Lee bounds to random objects. In contrast to the above papers, to the best of our knowledge, this is the first paper that conducts partial identification analysis on the Fr\'echet mean in a general metric space. In this sense, this paper opens a new avenue for applications of the methodologies of partial identification and metric statistics.
This paper is organized as follows. After closing this section with a recap of the conventional Lee bounds, Section (ref) showcases our partial identification analysis by focusing on compositional data. Then Section (ref) presents our general methodology for partial identification. In Section (ref), we discuss estimation and bootstrap inference on the proposed identified set. Section (ref) provides an empirical illustration. Finally, Appendix (ref) discusses additional examples of random objects, and Appendix (ref) contains the proofs of the theoretical results.
For the reader's convenience, we start with the standard Lee bounds for a scalar outcome. We also discuss potential issues that arise and are resolved in the subsequent sections. Consider a binary treatment $D\in \{0,1\}$ (one for the treatment and zero otherwise) with potential outcomes $Y_0$ and $Y_1$ that take values in $\mathcal{Y}\subset\mathbb{R}$. Let $S_{0}, S_{1}\in\{0,1\}$ denote potential selection indicators (one for non-missing and zero otherwise). As in Lee:2009, we assume that $D$ is independent from $(Y_1, Y_0, S_1, S_0)$ (random assignment) and $S_1 \ge S_0$ (monotonicity). Instead of the latent variables $(Y_1, Y_0, S_1, S_0)$, we observe $D$ and $(S, Y)$ that are generated from the general sample selection model:
When the outcome variable $Y$ is scalar-valued, we typically partially identify the average treatment effect for a certain subpopulation Lee:2009:
where we refer to the group of individuals with $S_{1} = S_{0} = 1$ as the always-observed group. The second term $\E{Y_{0} \mid S_{1} = S_{0} = 1}$ is point identified as $\E{Y \mid S = 1, D=0}$, but the first term $\E{Y_{1} \mid S_{1} = S_{0} = 1}$ is not. To partially identify the first term, observe that
where $p=\P{S=1 \mid D=0}/\P{S=1 \mid D=1}$. Then, by considering the worst (resp. best) case scenario where the smallest (resp. largest) $100\times p\%$ values of observed $Y$ (conditional on $S=1, D=1$) are entirely attributed to $Y_{1}$ of the always-observed group, we can obtain the lower (resp. upper) bound on the first term, leading to the following bounds:
where $Q_{Y}(\cdot)$ denotes the conditional quantile function of $Y$ given $D = 1, S = 1$. Based on these bounds, we can obtain the Lee bounds on the average treatment effect $\E{Y_{1} - Y_{0} \mid S_{1} = S_{0} = 1}$ for the always-observed group.
In the subsequent sections, we extend this partial identification strategy to outcomes that take values in more general metric spaces (i.e., random objects rather than Euclidean variables). One of the key issues is how to define and interpret a “mean” in such settings. For Euclidean outcomes, the usual expectation admits a geometric interpretation as the “center” of the distribution. This interpretation, however, is inherently metric-dependent, and as a result, when the outcome variable lies in a non-Euclidean space, the conventional Euclidean mean may ignore the geometric property and perhaps lead to a misleading notion of centrality.
To illustrate this issue, consider a three-part compositional outcome. Let $c_0=(1/3,1/3,1/3)$ denote the “barycenter” of the space of compositional data, and consider an interior point $c_1 = (2/3, 1/6, 1/6)$ and a boundary point $c_2 = (0, 1/2, 1/2)$. Under the Euclidean metric, as illustrated in Figure (ref), the distances from the barycenter to these two points are identical: $d(c_0, c_1) = d(c_0, c_2)$.
From a compositional perspective, however, these two displacements are qualitatively different even if they have the same Euclidean magnitude. The move from $c_0$ to $c_1$ remains within the interior of the simplex and corresponds to a finite reallocation of mass among components. In contrast, the move from $c_0$ to $c_2$ pushes the first component to zero, i.e., to the boundary (or an extreme point) of the simplex. This latter behavior can be considered more extreme. Accordingly, a metric that respects this kind of compositional geometry, such that $d(c_0, c_1) < d(c_0, c_2)$ as in Figure (ref) based on the Aitchison metric, is more natural.
The same consideration applies in a more general setting considered in Section (ref). Therefore, to extend Lee's (2009) bounding procedure to random object outcomes, we begin by introducing an appropriate metric that allows statistical analysis to reflect the geometry of the space in which outcomes take values. We then define an extended notion of the mean (i.e., Fr\'echet mean) with respect to the metric and develop the Lee-type identified region for the Fr\'echet mean of random object outcomes. Another key issue is how to operationalize the bounding approach. We address this issue by exploiting a well-behaved embedding that reduces our partial identification problem in a metric space to the one in a familiar Euclidean or Hilbert space.
Compositional data play a central role in economic analysis, as many key economic objects are inherently defined as allocations across mutually exclusive categories that sum to a whole. For example, household time allocation among labor, leisure, and housework is a fundamental concept for understanding labor supply decisions and intrahousehold bargaining. Analogously, the allocation of total expenditure across goods underpins classical consumer choice theory. In evaluating cash transfer programs, researchers likewise examine how households allocate consumption across categories; for example, the share devoted to health-related versus non-health food expenditures.
To accommodate these important scenarios, we now extend the scope of the Lee bounds to a compositional data scenario. Recall the treatment indicator ($D$) and selection indicators ($S_1, S_0, S$) defined as before. We retain the random assignment and monotonicity assumptions. Suppose that the potential outcomes $Y_0$ and $Y_1$ take values in the space for compositional data with positive components:
We impose the strict positivity condition $y_j > 0$ only to simplify the discussion that follows. Cases with zero components are accommodated within our general framework in Section (ref).
To introduce a notion of the expectation of the potential outcome $Y_t \in \mathcal{Y}$ for $t\in\{0,1\}$, let us consider the Aitchison metric $d$ on $\mathcal{Y}$ defined as
where $\|\cdot\|$ is the Euclidean metric on $\mathbb{R}^k$ and
The Aitchison metric is commonly applied for compositional data analysis (Aitchison:1982). As illustrated in Figure (ref), with the Aitchison metric, the resulting distances reflect the underlying geometric structure. Based on the metric space $(\mathcal{Y},d)$, we introduce the Fr\'echet mean of $Y_t\in\mathcal{Y}$, which is defined as
The Fr\'echet mean is a direct generalization of the conventional mean for the Euclidean space toward a general metric space, and its statistical analysis has been increasingly popular in recent literature (see, e.g., Dubey:2024). Analogously, the conditional Fr\'echet mean is defined as $\mathbb{E}_\oplus[Y_t \,| \,\cdot\,]=\mathrm{arg\,min}_{y \in \mathcal{Y}}\E{d^2(y, Y_t)\,|\,\cdot\,}$, and the average treatment effect for the always-observed group is characterized by some contrast of $\mu_{1\oplus}$ and $\mu_{0\oplus}$, where $\mu_{t\oplus}=\mathbb{E}_\oplus[Y_t\mid S_1=S_0=1]$.
As in the scalar outcome case, the assumptions of random assignment and monotonicity imply that $\mu_{0\oplus}$ is point identified as
In contrast, $\mu_{1\oplus}$ is not point identified as in the scalar outcome case. Thus, our main task is to characterize the sharp identified set for $\mu_{1\oplus}$. To this end, we note that (i) the map $\Psi:\mathcal{Y} \to \mathbb{R}^k$ is an isometric and injective map from $\mathcal{Y}$ to $\mathbb{R}^k$, and (ii) the image space $\Psi(\mathcal{Y}) = \{x \in \mathbb{R}^k: \sum_{j=1}^k x_j =0\}$ is a closed convex subset of $\mathbb{R}^k$. Due to these properties, $\mu_{1\oplus}$ can be written as
Note that the expectation $\mathbb{E}[\cdot]$ on the right-hand side is the one for the conventional Euclidean vector, and this representation motivates the following two-step procedure to construct an identified set for $\mu_{1\oplus}$:
For Step 1, a similar reasoning to the scalar-valued case applies, and the sharp identified set of $\E{\Psi(Y_1)\mid S_1=S_0=1}$ can be written as
where $\mathscr{A}$ is the collection of all subsets $\mathcal{A}$ such that $\P{\Psi(Y)\in \mathcal{A} \mid D = 1, S = 1} = p$ (the set in ((ref)) coincides with the set $\mathcal{I}_1$ in our general result in Section (ref)). Furthermore, a standard result in convex analysis implies that the identified set in ((ref)) admits an alternative tractable characterization $\mathcal{S}_1$. For Step 2, we can show that $\Psi^{-1}(\mathcal{S}_1)$ is indeed the sharp identified set of $\mu_{1\oplus}$, so the main result of this section is summarized as follows.
This proposition is obtained as a special case of the general result in the next section. Before turning to our general result, we illustrate the proposed strategy and its appeal using a numerical example.
This subsection illustrates our proposed procedure. We consider a hypothetical randomized controlled trial setting with sample attrition using data from the 2024 American Time Use Survey. We restrict attention to respondents aged 25-55 who are employed and report positive time spent on leisure, labor-market work, and housework.\footnote{We define leisure as the sum of first-tier codes 1 (Personal Care), 11 (Eating and Drinking), 12 (Socializing, Relaxing, and Leisure), and 13 (Sports, Exercise, and Recreation). Labor-market work corresponds to first-tier code 5 (Work & Work-Related Activities). We define housework as the sum of first-tier codes 2 (Household Activities), 3 (Caring for & Helping Household Members), 4 (Caring for & Helping Nonhousehold Members), and 7 (Consumer Purchases). See \url{https://www.bls.gov/tus/lexicons/lexiconwex2024.pdf} for the list of detailed activities within each category.} The resulting sample size is $1,397$. We then randomly assign individuals to two groups, treating the split as if it were a treatment-control assignment ($711$ in the control group and $686$ in the treatment group). To mimic attrition, we randomly drop observations according to group-specific selection probabilities. Specifically, we set $\P{S=1\mid D=1}=0.90$ and vary the control-group retention rate over $\P{S=1\mid D=0}\in\{0.85,0.70\}$. Since $\P{S=1\mid D=1}$ is held fixed, a smaller $\P{S=1\mid D=0}$ implies a smaller $p$, which in turn increases the degree of “contamination” $1-p$ and limits identifying power. Because there is, in fact, no treatment, we expect no treatment effect and hence assess whether the estimated bounds $\widehat{\mathcal{S}}_1$ contain the point-identified Fr\'echet mean for the controlled group $\widehat{\mu}_{0\oplus}$.
Figure (ref) presents the estimated bounds $\widehat{\mathcal{S}}_1$ (inner, purple region) and point-identified Fr\'echet mean $\widehat{\mu}_{0\oplus}$ (square dot), as well as the 95% bootstrap confidence region for ${\mathcal{S}}_1$ (outer, gray region), whose validity is shown in Section (ref). We can see that the bounds and confidence sets contain $\widehat{\mu}_{0\oplus}$, suggesting the validity of our proposed bounds. In the low-attrition case (Figure (ref)), the identified set is reasonably tight and informative. When the attrition rate is higher (Figure (ref)), the identified set expands, but it still provides non-trivial identifying power.
Table (ref) shows the projections of the estimated bounds and confidence region, as well as the naively applied standard Lee bounds. In Panel A, we can find that the componentwise Lee bounds do not contain the Fr\'echet mean for the control group (which coincides with the counterfactual Fr\'echet mean as well in our scenario) in any dimensions. This pattern suggests that the standard Lee bounds may fail to deliver a valid identified set for the center of the counterfactual outcome distribution (as measured by the Fr\'echet mean). By contrast, the componentwise projections of our proposed bounds do cover $\widehat{\mu}_{0\oplus}$. We also note that these confidence intervals are simultaneously valid, as they are obtained by projecting the confidence region for $\Psi^{-1}(\mathcal{S}_1)$.
In Panel B, although the naive Lee bounds become wider, the bounds for the leisure component still fail to cover $\widehat{\mu}_{0\oplus}$. In contrast, the componentwise projections of our proposed bounds remain valid, although they are less informative than in the low-attrition case shown in Panel A.
This section extends the benchmark example in the last section and presents a general framework to conduct partial identification analysis for the sample selection model in ((ref)), where the outcome variable $Y$ resides in a metric space $(\mathcal{Y},d)$.\footnote{Although we focus on the sample selection model, our analysis also applies to the contaminated data model studied by Horowitz_Manski:1995; see Remark (ref) below.} To begin with, we redefine the notation and formally introduce our setup.
This setup is identical to the one in Lee:2009 except that the support $\mathcal{Y}$ can be non-Euclidean, and we are concerned with causal inference by a randomized experiment in the presence of sample selection. The monotonicity assumption $S_1 \ge S_0$, which rules out “defiers” ($S_1=0, S_0=1$) from the population, requires that the assignment $D$ can affect the sample selection only in one direction. This assumption is standard in many econometric models and aligns well with the threshold-crossing model Vytlacil:2002.
In this setup, the average treatment effect of interest is characterized by some contrast of the Fr\'echet means \[ \mu_{t\oplus}=\mathbb{E}_\oplus[Y_t\mid S_1=S_0=1]=\mathop{\rm arg~min}\limits_{y \in \mathcal{Y}}\E{d^2(y,Y_{t})\mid S_{1}=S_{0}=1}, \] for $t \in \{0,1\}$. Again, $\mu_{0\oplus}$ is point identified as in ((ref)) under Assumption (ref). Thus, we focus on partial identification of $\mu_{1\oplus}$, and impose the following assumptions on the metric space $(\mathcal{Y},d)$.
Assumption (ref) (ii) is a regularity condition requiring the existence of the second moment and continuity of the distribution. The key assumption in our analysis is hence Assumption (ref) (i). Assumption (ref) (i) enables us to characterize the (conditional) Fr\'echet mean of the random variable $Y$ using the map $\Psi$. In fact, the (conditional) Fr\'echet mean of $Y$ is obtained by pulling back, via $\Psi^{-1}$, the (conditional) expectation of the embedded variable $\Psi(Y)$ in the Hilbert space $\mathcal{H}$, which is uniquely determined by the Riesz representation theorem. We refer to kuri:25 for further details on this point. Many of the common metric spaces studied in the literature of metric statistics are known to satisfy Assumption (ref) (i). Here, we list some examples of metric spaces that are commonly used in statistical analysis of metric space-valued data. Additional examples are contained in Appendix (ref).
Now, we are in a position to derive the sharp identified set of $\mu_{1\oplus}=\mathbb{E}_\oplus[Y_1\mid S_1=S_0=1]$ (see Definition (ref) for the formal definition of the sharp identified set). Under Assumption (ref), the potential Fr\'echet mean can be expressed as $\mu_{1\oplus}=\Psi^{-1}(\E{\Psi(Y_{1})|S_{1}=S_{0}=1})$. Thus, as in the last section, we begin by characterizing an identified set for $\E{\Psi(Y_{1})|S_{1}=S_{0}=1}$. Let $p= \P{S=1 | D=0}/\P{S=1 | D=1}$ and observe that \[ \P{\Psi(Y)\in\mathcal{A}|D = 1, S = 1} = p\P{\Psi(Y_{1})\in\mathcal{A}|S_{1} = S_{0} = 1} + (1-p)\P{\Psi(Y_{1})\in\mathcal{A}|S_{1} = 1, S_{0} = 0}, \] where $\mathcal{A}$ is an arbitrary measurable subset of $\mathcal{H}$. Hence, a $100\times p\%$ of the observed subpopulation with $D = 1, S = 1$ corresponds to the always-observed group. Then, by letting $\mu(\mathcal{A})\coloneqq \P{\Psi(Y)\in\mathcal{A}|D = 1, S = 1}$, the sharp identified set for $\E{\Psi(Y_1) \mid S_1=S_0=1}$ is obtained as
where $\int_{\mathcal{H}} \cdot \,d\nu$ means the Bochner integral, and the sharpness of $\mathcal{I}_1$ is shown in Appendix (ref). Our partial identification results are presented as follows.
Theorem (ref) establishes a Lee-type identified set for the conditional Fr\'echet mean $\mu_{1\oplus}$ of general random objects. Theorem (ref) (iii) shows that, when $\Psi(Y)$ is finite-dimensional, the proposed procedure delivers a sharp identified set for $\mu_{1\oplus}$. Representing a closed convex set as in (ref) is standard in econometrics (e.g., BeMo08), and this representation is computationally feasible in practice. Further details are discussed in Section (ref). We note that incorporating covariates is straightforward (see Remark (ref) below).
When $\Psi(Y)$ is infinite-dimensional, evaluating (ref) may be infeasible, and $\mathcal{S}_1$ need not coincide with the sharp identified set for $\mu_{1\oplus}$. Nevertheless, one may still be able to obtain a sharp identified set over a finite grid of points by applying the same procedure to a finite-dimensional projection of $\Psi(Y)$, which is useful to approximate $\mathcal{S}_1$ or to obtain a feasible, valid identified set. See Remark (ref) below for the case of distributional outcomes.
We conclude our identification analysis with a remark on the contaminated and corrupted data models in the sense of Horowitz_Manski:1995.
Once the sharp identified set for $\mu_{1\oplus}$ is obtained, the remaining question is how to summarize the treatment effect. Unlike the Euclidean-outcome case, when outcomes take values in a general metric space $\mathcal{Y}$, one cannot in general define treatment effects by algebraic differencing of potential means, since $\mathcal{Y}$ need not admit a linear structure. Nevertheless, researchers can define treatment effects in ways that are tailored to their substantive objectives. Indeed, there can be several approaches to defining causal effects, which we list below.
Difference in the embedded space. The first possibility is to define the causal effect through differences in an embedded space. For a general metric space that can be embedded into a Hilbert space $\mathcal{H}$, we can define the contrast $\Psi(\mu_{1\oplus})-\Psi(\mu_{0\oplus})$ as the causal effect of interest, provided that the embedded object $\Psi(Y)$ has a reasonable interpretation.
A leading example is the case in which $\mathcal{Y}$ is the $2$-Wasserstein space for distributional outcomes. In this setting, it is common to work with quantile-function representations and to define a quantile-treatment-effect-type parameter as the difference in the embedded space of quantile functions $\Psi(\mathcal{Y})$ (e.g., LiKoWa23 and va:25):
In this case, the identified set for the embedded contrast can be written as the Minkowski difference
Difference in the projections. A second option is to focus on a real-valued contrast obtained by projecting the mean objects $\mu_{0\oplus}$ and $\mu_{1\oplus}$ onto $\mathbb{R}$. For example, suppose $Y$ is a compositional outcome and the researcher's interest is in the percentage-point change in the first component. Let $\Pi:\mathcal{Y}\to\mathbb{R}$ denote the projection onto the first coordinate, i.e., $\Pi(y)=y_1$ for $y=(y_1,\ldots,y_k)\in\mathcal{Y}$. Then a natural causal estimand is $\Pi(\mu_{1\oplus})-\Pi(\mu_{0\oplus})$. More generally, the same construction applies whenever $\Pi:\mathcal{Y}\to\mathbb{R}$ has a meaningful interpretation in the application at hand. In this case, an identified set for the projected contrast is given by
Geodesics. As a third option, one may employ the geodesic average treatment effect (GATE) proposed by kurisu_etal:24. The GATE is defined as the entire geodesic path from $\mu_{0\oplus}$ to $\mu_{1\oplus}$:
Intuitively, for $\alpha, \beta\in\mathcal{Y}$, the geodesic $\gamma_{\alpha,\beta}$ is defined as the shortest path from $\alpha$ to $\beta$. For a precise definition of geodesics in a general metric space and for the relationship between GATE and the conventional ATE defined for Euclidean outcomes, see kurisu_etal:24. The full geodesic includes both the direction and the magnitude of the change induced by the treatment. This provides insight into how the underlying object evolves along the shortest intrinsic path. Based on the identified set $\Psi^{-1}(\mathcal{S}_1)$ for $\mu_{1\oplus}$, the identified set of the GATE $\gamma_{\mu_{0\oplus},\mu_{1\oplus}}$ is given by the following set of geodesics: \[ \left\{ \gamma_{\mu_{0\oplus},\mu} : \mu\in\Psi^{-1}(\mathcal{S}_1) \right\}. \] Provided that $\mathcal{S}_1$ is sharp, the induced identified sets for these causal estimands are sharp as well.
This section provides detailed estimation and inference procedures for finite-dimensional outcomes. In the infinite-dimensional case, the same approach applies after projecting the object of interest onto a finite-dimensional parameter.
We begin with the estimation procedure. Recalling that
the Fr\'echet mean $\mu_{0\oplus}$ can be estimated by its sample analog, or $\Psi^{-1}(\widehat{\mathbb{E}}[\Psi(Y) \mid S=1, D=0])$, where $\widehat{\mathbb{E}}[\Psi(Y) \mid S=1, D=0])$ is the standard sample conditional mean in the Euclidean space. The sharp identified set of $\mu_{1\oplus}$ can be estimated by
where $\mathbb{S}^{d-1}=\{x \in \mathbb{R}^d:\|x\|=1\}$ is the $(d-1)$-dimensional unit sphere,
$\widehat{p} = \{(\sum_{i=1}^n S_i(1-D_i))/(\sum_{i=1}^n (1-D_i))\}/\{(\sum_{i=1}^n S_i D_i)/(\sum_{i=1}^n D_i)\}$, and $\widehat{Q}_{\langle u, \Psi(Y)\rangle}(q)$ denotes the $q$-th sample conditional quantile of $\langle u, \Psi(Y)\rangle$ given $D = 1, S = 1$. In practice, the intersection “$\bigcap_{u\in \mathbb{S}^{d-1}}$” in (ref) can be approximated arbitrarily well by evaluating (ref) on a sufficiently fine grid of points on the unit sphere $\mathbb{S}^{d-1}$. When $d$ is small, one may use a grid with (approximately) equal angular spacing. For example, see BeMoMo10 for the case $d=2$; for $d=3$, one may employ a Fibonacci lattice. In higher-dimensional cases, an approximation based on random number generation is often more convenient. For instance, if $Z\sim \mathcal{N}(\mathbf{0}, \mathbf{I}_d)$, then $Z/\|Z\|$ is uniformly distributed on the $(d-1)$-dimensional unit sphere Devroye:1986, and thus can be used to form an approximation of the intersection.
This subsection discusses a method for uncertainty quantification. For $\mu_{0\oplus}$, the quantity $\E{\Psi(Y)\mid S=1,D=0}$ is point identified and corresponds to a standard Euclidean conditional mean. Accordingly, conventional procedures (e.g., the bootstrap) apply, and one can construct a confidence region for $\E{\Psi(Y)\mid S=1,D=0}$, denoted by $\widehat{\mathcal{R}}_0$, in a standard manner. The corresponding confidence region for $\mu_{0\oplus}$ is then given by $\Psi^{-1}(\widehat{\mathcal{R}}_0)$.
We therefore focus on uncertainty quantification for the sharp identified set $\mathcal{S}_1$. Let $\alpha\in(0,1)$ be a researcher-specified significance level. For each $b\in\{1,\ldots, B\}$, compute
where $\widetilde{\sigma}_{\mathcal{I}_1}^{(b)}(u)$ is the bootstrap counterpart of $\widehat{\sigma}_{\mathcal{I}_1}(u)$, and $\widehat{V}_{\sigma}(u)$ is the estimator of the asymptotic variance $V_{\sigma}(u)$ of $\widehat{\sigma}_{\mathcal{I}_1}(u)$, which is given by Lee:2009. In practice, this supremum is approximated by the maximum over some fine grid points. Then, we compute the critical value
and obtain the confidence region of $\mathcal{S}_1$ as
where $\mathbb{S}^{d-1}$ is the $(d-1)$-dimensional unit sphere. The following proposition establishes the validity of this procedure.
In view of the proof of Proposition (ref), the joint confidence region for $(\Psi(\mu_{0\oplus}), \mathcal{S}_1)$ is also readily available. Instead of (ref), define
where $\widehat{\xi}(j)$ and $\widetilde{\xi}^{(b)}(j)$ are the sample- and bootstrap- estimators of $\xi(j)\coloneqq$ the $j$-th element of $\E{\Psi(Y)\mid S=1,D=0}$, and ${V}_{0}(j)$ is the asymptotic variance of $\widehat{\xi}(j)$ with a consistent estimator $\widehat{V}_{0}(j)$. Then, using the critical value $\widehat{\text{j.c.v.}} =$ the $(1-\alpha)$-th quantile of $\{T_{\mathtt{joint}}(b) : b=1,\ldots,B\}$, one can construct a joint confidence region such that
where
and $\widehat{\mathcal{R}}_{1,\mathtt{joint}}$ is constructed as before but replacing $\widehat{\text{c.v.}}$ with $\widehat{\text{j.c.v.}}$.
This general construction of the joint confidence region will be useful for conducting a hypothesis test of the treatment effects (see Section (ref)). When the treatment effect is defined as the difference in the embedded space $\Psi(\mathcal{Y})$, the $1-\alpha$ confidence region for the treatment effect is given by $\widehat{\mathcal{R}}_{1,\mathtt{joint}} \ominus \widehat{\mathcal{R}}_{0,\mathtt{joint}}$. If the treatment effect is measured by the projection difference, $\Pi(\Psi^{-1}(\widehat{\mathcal{R}}_{1,\mathtt{joint}})) \ominus \Pi(\Psi^{-1}(\widehat{\mathcal{R}}_{0,\mathtt{joint}}))$ provides a valid confidence interval. When the geodesic treatment effect is employed, $\{ \gamma_{\nu,\mu} : \nu \in \Psi^{-1}(\widehat{\mathcal{R}}_{0,\mathtt{joint}})\,\text{and}\, \mu\in\Psi^{-1}(\widehat{\mathcal{R}}_{1,\mathtt{joint}})\}$ has a correct asymptotic coverage.
This section illustrates our proposed method using a distributional outcome: sleep duration (hours slept per night). In particular, we study how individuals' sleep responds to financial incentives. Sleep is a central component of time allocation and a potentially important input into economic performance (e.g., Rao_etal:2021, Streatfeild_etal:2021), and recent work in economics has examined whether incentives can shift sleep behavior (e.g., Avery_etal:2025, Bessone_etal:2021). We use a subset of the data from Bessone_etal:2021 and reanalyze the experiment with our method, focusing on distributional changes in sleep hours.
In Bessone_etal:2021, low-income adults in India were employed in a paid data-entry job and wore an actigraphy-based wearable device throughout the study. After eight baseline days, participants were randomly assigned to one of three arms (control, sleep-aid devices$+$incentives, and sleep-aid devices$+$encouragement) and were further cross-randomized to a workplace nap offer. In the incentives arm, participants received sleep-improving devices and were offered monetary incentives tied to increases in actigraph-measured sleep relative to baseline. In this illustration, we focus on the incentives-versus-control comparison and pool over the nap assignment. As a result, the sample sizes are 152 in the control group and 150 in the treatment (incentives) group.
As the distributional sleep outcome, we use the nightly sleeping hours from days 9-28 and treat the resulting set of hours for each individual as a distributional outcome. For illustration, we classify an individual as a non-respondent if sleep data are missing for at least three nights out of these 20 nights. An assumption here is that, when the number of missing nights is below the threshold, missingness is at random and therefore does not distort the distribution; otherwise, the missingness is potentially endogenous. Under this setup, the attrition rate is 12.5% in the control group and 10.7% in the treatment group.
We estimate the identified set following the procedure described in Remark (ref) and Section (ref). We choose $k=17$ evaluation points, with $q_1=0.10, q_2=0.15, \ldots, q_{17}=0.90$, and estimate $\mathcal{S}_{1,\Pi}$. We then randomly draw $200{,}000$ points $\{(f(q_1),\ldots,f(q_{17}))_j\}_{j=1}^{200{,}000}$ from $\widehat{\mathcal{S}}_{1,\Pi}$, and construct $200{,}000$ linearly interpolated curves connecting $f(q_1),\ldots,f(q_{17})$. The 95% joint confidence region is computed following Corollary (ref) and visualized in the same manner as $\widehat{\mathcal{S}}_{1,\Pi}$.
The estimation results are shown in Figure (ref). We observe several notable features. First, the estimated identified set (inner, purple band) is tight and informative. This is because $\widehat{p}=(1-0.125)/(1-0.107)\approx 0.98$ is close to one, implying that the observed distribution is only mildly “contaminated” and the identified set is therefore close to point identification. Second, the identified set lies above the control mean (red curve), and the joint confidence region shows no overlap between the treatment and control groups. These findings suggest that providing incentives shifts the distribution of sleep upward. Moreover, the observations are consistent with $\mu_{1\oplus}$ first-order stochastically dominating $\mu_{0\oplus}$, which is a richer implication than a comparison of average sleep duration alone.