EconBase
← Back to paper

Prediction Sets and Conformal Inference with Interval Outcomes

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.

94,625 characters · 11 sections · 63 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Prediction Sets and Conformal Inference with Interval Outcomes

\if10 \fi \if00\fi

abstract\fontsize{11}{12}\selectfont Given data on a random variable \(Y\), a prediction set with miscoverage level \(\alpha \in (0,1)\) is a set that contains a new draw of \(Y\) with probability \(1-\alpha\). Among all prediction sets satisfying this coverage property, the oracle prediction set is the one with minimal volume. The oracle prediction set offers a complementary view of the distribution of \(Y\), beyond point estimators such as the mean and quantiles, and has attracted considerable interest recently. This paper develops methods for estimating such prediction sets conditional on observed covariates when \(Y\) is censored or interval-valued. We characterise the oracle prediction set under partial identification induced by interval censoring and propose consistent estimators for both oracle prediction intervals and more general oracle prediction sets consisting of multiple disjoint intervals. In addition, we apply conformal inference to construct finite-sample valid prediction sets for interval outcomes that remain consistent as the sample size grows, using a conformity score tailored to interval data. The proposed procedure accounts for irreducible prediction uncertainty due to the stochastic nature of outcomes, modelling uncertainty arising from partial identification, and sampling uncertainty that vanishes as sample size increases. We conduct Monte Carlo simulations and two empirical applications using UK job postings data and the US Current Population Survey. The results demonstrate the robustness and efficiency of the proposed methods. {\bf Keywords:} {Prediction sets; Partial identification; Nonparametric estimation; Conformal Prediction.}

\onehalfspacing

Introduction

Interval data are pervasive. Surveys often resort to “bracketing” strategies to address item non-response. For example, respondents are often asked whether their money in savings accounts, in investments, or their income falls into a sequence of brackets such as $[\$10K, \$20K]$, $[\$20K, \$50K]$, etc, while also being given the option to provide a point answer to such questions (see heeringa1995unfolding and moore2001using). Interval-valued data also accommodate common patterns such as right censoring, where an outcome variable is greater than a specified value (i.e., $Y \in [Y^L, \infty)$); left-censoring, where an outcome variable is less than a certain value (i.e., $Y \in (-\infty, Y^U]$); and data from competing risks, where one typically observes the $\max$ or $\min$ of variables of interest. Such data typically involves outcomes of interest such as income, wealth, unemployment duration, and so on. This paper provides a new method for predicting an outcome variable of interest $Y$ when we observe a sample of intervals $[Y^L, Y^U]$ such that $Y \in [Y^L, Y^U]$ rather than direct observations of $Y$.

More specifically, suppose we have a random variable $Y$ defined on $\mathcal Y \subset \mathbb{R}$ with distribution $P_Y$. One parameter of interest is a {\it prediction set} $C\subset \mathcal{Y}$ with miscoverage level $\alpha \in (0,1)$, which is any measurable subset of $\mathcal{Y}$ such that $P_Y( Y \in C) \geq 1-\alpha$. The oracle prediction set is the set with the minimal volume that has this coverage property, which is related to the level set of the density of $Y$ when this density exists. With covariates, we have access to a sample $(X_{i}, Y_{i}) \in \mathcal{X} \times \mathcal{Y}$, $i = 1,\dots,n$, with supports $\mathcal{X} \subset \mathbb{R}^{d}$ and $\mathcal{Y}\subset \mathbb{R}$. We use these data to construct prediction sets $C_n(x)\equiv C(x; X_{1}, Y_{1}, \dots, X_{n},Y_{n}) \subset \mathcal{Y}$ for a given miscoverage level $\alpha$ at $x \in \mathcal{X}$. It is desirable that such a prediction set $C_n$ is both consistent for the oracle prediction set and has a {\it finite-sample coverage property:}

equation[equation omitted — 93 chars of source]

where $\Prob$ is the distribution over $(X_{i}, Y_i)$ for $i= 1,\dots, n+1$. This paper constructs such a set when $Y$ is interval-valued. The exact finite-sample coverage property that our construction provides is given in more detail in (ref).

Such an “oracle prediction” exercise is familiar in the literature, especially in the context of forecasting in linear regressions (see diebold2015ForecastingEconomics). For example, assuming that the outcome $Y$ satisfies $Y_{i}\mid X_{i} \sim \mathcal N(X_{i}'\beta, \sigma^2)$ with known $\sigma^{2}$, the population oracle prediction set is $X_{n+1}'\beta \pm 1.96\sigma$. Given data $\mathcal{D} = \cqty{(X_1, Y_1), \ldots, (X_n, Y_n)}$, a consistent estimator for such a set is $X_{n+1}'\hat\beta \, \pm \, 1.96\sigma$, where $\hat \beta$ is the least squares estimator of $\beta$. On the other hand, “an operational density forecast that accounts for parameter uncertainty” is (diebold2015ForecastingEconomics): \[ N\left(X_{n+1}' \hat{\beta}, \ {\sigma}^2 \left(1 + X_{n+1}' (\mathbf{X}'\mathbf{X})^{-1} X_{n+1} \right)\right). \] This “operational density” summarises both the sampling uncertainty in estimating $\hat \beta$ and the disturbance uncertainty arising from the intrinsic randomness in the conditional distribution $Y_{n+1}\mid X_{n+1}$. One can use it to derive a forecast density with a particular coverage property.

The prediction set is closely related to quantile methods if we have direct observations of the outcome variable $Y$. If the oracle prediction set is a single interval, it is given by $C(X_{n+1}) \equiv [\tau_{0}(X_{n+1}), \tau_{1}(X_{n+1})]$, where $\tau_0(X_{n+1})$ and $\tau_1(X_{n+1})$ are appropriately chosen quantiles for $Y$ given $X$ so that one can obtain a coverage rate $(1-\alpha)$ with minimal volume. Taking the Gaussian specification above, with $\alpha=5\%$ and assuming $\beta$ is known, $\tau_0(X_{n+1}) = X_{n+1}'\beta - 1.96\sigma$ and $\tau_1(X_{n+1}) = X_{n+1}'\beta + 1.96\sigma$, where $1.96\sigma$ is the $95\%$-quantile of $|\varepsilon_{n+1}| \equiv |Y_{n+1}-X_{n+1}'\beta|$. When $\sigma$ is known as presumed above, but $\beta$ is not, estimates for these can be obtained as $\hat \tau_0(X_{n+1}) = X_{n+1}'\hat \beta - 1.96\sigma$ and $\hat \tau_1(X_{n+1}) = X_{n+1}'\hat \beta + 1.96\sigma$, where $\hat \beta$ is the OLS or ML estimate of $\beta$. This nonetheless ignores the parameter uncertainty in estimating $\beta$. This is not a problem if we are just interested in interval consistency, but as noted in romano2019ConformalizedQuantile, replacing $\beta$ with $\hat \beta$ “is not guaranteed to satisfy the coverage statement (\dots) when $C(X_{n+1})$ is replaced by the estimated interval $\hat C(X_{n+1})$.” To address this shortcoming, romano2019ConformalizedQuantile employ conformal inference and modify the interval $\hat C(X_{n+1}) = [\hat \tau_0(X_{n+1}),\hat \tau_1(X_{n+1})]$ by adjusting the endpoints to $\hat \tau_0(X_{n+1})-\vartheta_{1-\alpha}$ and $\hat \tau_1(X_{n+1})+\vartheta_{1-\alpha}$. Here, $\vartheta_{1-\alpha}$ is the $(1-\alpha)(1+1/n_2)$-th empirical quantile of the set of conformity scores $\cqty{\max\{\hat \tau_0(X_i)-Y_i,Y_i-\hat \tau_1(X_i)\}: i\in \mathcal{I}_2}$ estimated in a calibration sub-sample $\mathcal{I}_2$ with sample size $n_2$. This corrects for under- or over-coverage and allows for the construction of prediction intervals satisfying desired coverage levels in finite samples using quantile estimates.

However, several questions and limitations surface when the outcome variable $Y$ is interval-censored. First, the definition of the oracle prediction set needs to be extended in partially identified models. We also note that although we can construct a valid prediction interval using quantile estimates, such intervals are not necessarily minimal in volume without assumptions on the shape of the distribution (such as the Gaussian assumption in the example above), and this approach is not applicable when the oracle prediction set consists of multiple disjoint intervals. Furthermore, a new conformity score is needed to accommodate interval-valued observations. Our method addresses these limitations and offers a complementary perspective on the distribution of outcomes beyond existing quantile-based methods.

We first define the oracle prediction sets under partial identification. We then characterise the oracle prediction sets under interval censoring via a feasible optimisation problem that relates the oracle prediction sets to the conditional distribution of the observed upper and lower brackets, $(Y^L, Y^U)$. We use a kernel-based conditional distribution function estimator as a basis for the construction of a consistent set estimator for the oracle prediction set. The oracle prediction set may be the union of disjoint intervals, for instance, when the prediction density is multimodal. We provide an approach that allows for the prediction set estimator to be a union of disjoint intervals with the prescribed coverage. In cases with a multimodal density, a union of intervals can have a smaller volume than one interval. We also allow for both random censoring with continuously distributed brackets and fixed censoring with discretely distributed brackets, such as those arising from “unfolding brackets” or specific survey schemes, see heeringa1995unfolding and moore2001using. See also the bracketing approach to outcomes using the Health and Retirement Survey in manski2002InferenceRegressions.

We use recent developments on conformal inference to propose a conformal prediction set $\tilde{C}_{n}$ using a sample of exchangeable data, which has finite-sample validity while also enjoying the asymptotic property of being consistent for the oracle prediction set under partial identification. Our conformal procedure is designed for interval-valued outcomes, which requires the design of a particular conformity score function to handle interval data. In particular, the score function is non-positive if and only if the constructed prediction set contains the interval outcome; this conformity score function therefore accounts for both over- and under-coverage, as in romano2019ConformalizedQuantile. As a competing alternative prediction set estimator, we show in the appendix that a conformal quantile regression method provides a valid prediction interval for interval-valued outcomes, which can be easier to implement in practice, but in general will not be efficient in terms of the volume of the prediction set. In addition, quantile regression-based methods are not suitable when the prediction set is not a single interval.

We illustrate our construction using a set of Monte Carlo experiments, where we consider both fixed and random censoring. Our proposed prediction set is shown to have robust finite sample coverage properties and to have a smaller volume compared with the alternative prediction sets constructed with quantile regression methods.

We apply our methods to two real-world datasets. The first is the UK job advert dataset from Adzuna, an online UK job advert aggregator. Adzuna data are widely used as an indicator of UK economic activity. One feature of these data is that many job ads report wages as ranges. We use these data to construct predictive wage intervals as a function of UK location and job categories. The second application uses the Current Population Survey (CPS) data. Among people who report their income, a significant fraction of the observations are interval-valued. Interval-valued observations are sometimes discarded or handled by imputation methods that rely on assumptions such as missing at random and may cause bias if misspecification happens. Our proposed prediction set is robust in that we do not require a specification of the censoring mechanism.

Literature review

The paper combines insights from the partial identification literature and conformal inference, which has recently been gaining significant interest. The partial identification literature (see manski2003PartialIdentification) takes the view that one should try to make inferences under minimal assumptions, even if these assumptions do not allow for point identification of the parameter of interest. The motivation for such a view is that inferences are less credible and conclusions are suspect when they rely on non-testable or non-credible assumptions. In our context for example, if we are interested in $\beta = E[Y]$ and we observe $[Y^L, Y^U]$ rather than $Y$, then a partial identification approach would just say that, given the data, the identified set for $\beta$ (or all we can learn about $\beta$ given the data) is the interval $\left[E[Y^L], E[Y^U]\right].$ Imputation based inference uses models to predict a $Y^P$ within $[Y^L, Y^U]$ and takes $E[Y^P]$ as a proxy for $\beta = E[Y]$. The credibility of this exercise depends on the imputation model, i.e. the set of assumptions it relies on. In the context of interval data, manski2002InferenceRegressions study identification of linear regression coefficients when the data (on regressors or outcomes) is interval valued using a partial identification approach. See tamer2010PartialIdentification, molinari2020Chapter5, and kline2023RecentDevelopments for reviews on the developments in the partial identification literature.

In order to obtain a prediction interval that has finite-sample coverage guarantees, we use conformal inference, which is a distribution-free method that builds upon the idea of exchangeability and permutation tests (vovk2005AlgorithmicLearning, vovk2009OnlinePredictive). The idea of conformal inference has become increasingly popular in the machine learning literature since it requires no assumptions on the underlying distribution other than the observations being exchangeable, and it also allows the use of any prediction protocol that treats the data symmetrically and hence preserves exchangeability. barber2023ConformalPrediction discusses the validity of conformal inference when the data is not i.i.d. and the predictor is not symmetric. Conformal inference has been applied to obtain valid distribution-free prediction sets lei2013DistributionfreePrediction, in non-parametric regression lei2014DistributionfreePrediction and in high-dimensional settings lei2018DistributionfreePredictive. romano2019ConformalizedQuantile considers conformal quantile regression to construct prediction intervals instead of conditional mean regression. sadinle2019LeastAmbiguous proposes a way to construct set-valued classifiers. chernozhukov2021DistributionalConformal considers the construction of valid prediction intervals with the distributional regression method. We add to this literature by allowing the outcome observations to be interval-censored or interval-valued. Relatedly,candes2023ConformalizedSurvival considers conformal inference for survival analysis. Using our notation, this amounts to observing either $Y_i$ or the interval $(Y_i^L, \infty)$ for each $i$. They propose a weighted conformal inference procedure tibshirani2019ConformalPrediction under the assumption that censoring is conditionally independent of her outcome given predictors $X_i$. We take a different approach and ask how much information we can extract from a sample with interval-valued outcomes. For example, our procedure allows for the case when all outcomes are interval censored (i.e., $Y_i$ is never observed) and the procedure uses the information from the upper bracket $Y^U_i$ when it is not $\infty$. Also, our procedure does not require the conditional independence assumption, which can be violated in many empirical contexts.

There are several related areas connected to our theoretical development. The classical way to construct estimators for the oracle prediction set relies on the level set estimation (see, for example, wilks1941DeterminationSample, polonik1995MeasuringMass and samworth2010AsymptoticsOptimal). The uniform consistency of the nonparametric estimator of the conditional distribution is an extension of classical results on the uniform convergence rate of the kernel-based density and regression estimators, such as einmahl2000EmpiricalProcess, gine2002RatesStrong and hansen2008UniformConvergence. This paper is also related to the literature on quantile regression (see koenker2017QuantileRegression for a review).

Outline

This paper is organised as follows. In (ref), we provide the setup and definitions of optimal sets and approaches to obtaining estimators of such sets with censoring. (ref) provides a study of the consistency of such sets using kernel estimators of the conditional distribution. (ref) provides results on finite sample coverage using a conformal inference procedure under interval censoring. (ref) provides some Monte Carlo evidence, and (ref) applies our methods to UK job postings data as well as the US CPS (Current Population Survey) Data. (ref) concludes. The proofs for the theorems in the main text, as well as some further results and comments, are collected in the Appendix.

Setup and definitions

Suppose there is a random sample \(\bar{\mathcal{D}} = \{\bar{Z}_{i} = (X_{i}, Y_{i}, Y_{i}^{L}, Y_{i}^{U}) \in \mathbb{R}^{d+3}: i \in \mathcal{I}\}\) with sample size \(\abs{\mathcal{I}} = n\) drawn independently from a joint distribution \(P\). We only observe the sub-vector \(Z_{i} = (X_{i}, Y_{i}^{L}, Y_{i}^{U})\), and \(Y_{i}\) is not directly observed but assumed to satisfy the condition \(Y_{i}^{L} \leq Y_{i} \leq Y_{i}^{U}\). For each $i\in \mathcal{I}$, the predictor \(X_i\) is a \(d\)-dimensional random vector with support \(\mathcal{X} \subset \mathbb{R}^{d}\). The latent outcome \(Y_{i}\) and the observed lower and upper bounds \(Y_{i}^{L}\) and \(Y_{i}^{U}\) are one-dimensional. Let \(\bar{Z} = (X, Y, Y^L, Y^{U})\) denote a generic sample from the joint distribution \(P\). We will use notation such as $P_Y$ to denote the marginal distribution of $Y$, and $P_{Y\mid X}$ to denote the distribution of $Y$ conditional on $X$. Note here that the joint distribution \(P_{X, Y^L, Y^U}\) of $(X,Y^{L}, Y^{U})$ is identified given the observed sample \(\mathcal{D} = \cqty{Z_ i = (X_i, Y_i^L, Y_i^U) \in \mathbb{R}^{d+2}: i \in \mathcal{I}}\), but the joint distribution \(P_{X, Y}\) of predictors $X$ and latent outcome $Y$ is not point identified given the observed sample $\mathcal D$ without restrictions on the censoring mechanism.

A prediction set is defined as a set $C\subset \mathcal{X} \times \mathcal{Y}$, with the interpretation that, given a realisation of the predictor $X = x$, we predict that $Y$ is likely to fall inside the section $C(x) = \cqty{y : (x, y) \in C}$. In the absence of interval censoring, estimation of a prediction set that satisfies certain optimality criteria (i.e. the minimal prediction set with coverage guarantees) has been studied in the literature; see lei2014DistributionfreePrediction, barber2021LimitsDistributionfree among others. As we discuss later in this section, the oracle prediction set without interval censoring is related to the upper level set of the density function of $Y$ given $X$ under point identification.

The situation is significantly different when we are faced with interval censoring and $P_{X,Y}$ is partially identified. In the rest of this section, we first define an optimality criterion for a prediction set under censoring, and analyse why the classical oracle prediction set (which corresponds to an upper level set) is no longer applicable under censoring. We then propose a feasible estimation procedure for the oracle prediction set. We discuss the potential reasons why the procedure can result in conservative prediction sets and show that the proposed procedure does not incur any unnecessary conservativeness.

There are two criteria for an oracle prediction set: validity and efficiency vovk2009OnlinePredictive,lei2013DistributionfreePrediction. Let $\mathcal{P}_{I}$ be the identified set of the joint distribution of $(X,Y)$ under censoring. That is, $\mathcal{P}_{I}$ is the set of joint distributions of $(X,Y)$ that are compatible with the joint distribution $P_{X,Y^{L},Y^U}$ and some censoring mechanism. A prediction set $C$ is {valid} if it satisfies the following coverage properties for a new sample $(X_{n+1}, Y_{n+1})$ drawn from the same distribution $P$ at a given miscoverage level $\alpha$.

definition[Validity under Partial Identification] We say a prediction set \(C\) is marginally valid under partial identification with miscoverage level $\alpha$, if, \begin{equation} \inf_{P\in \mathcal{P}_I} P \pqty{Y_{n+1} \in C(X_{n+1})} \geq 1 - \alpha, \end{equation} and \(C\) is conditionally valid under partial identification with miscoverage level $\alpha$ if for almost everywhere \(x\in \mathcal{X}\), \begin{equation} \inf_{P\in \mathcal{P}_I} P \pqty{Y_{n+1} \in C(x) \mid X_{n+1} = x} \geq 1 - \alpha. \end{equation}

When the model is point identified, \(\mathcal{P}_{I}=\cqty{P_{X,Y}}\) is a singleton set, the conditions in Equations (ref) and (ref) reduce to the standard conditions of validity, see for example, lei2014DistributionfreePrediction. We will simply refer to these two conditions as {marginal validity} and {conditional validity} henceforth.

It is natural to require a prediction set to be valid, so that we have confidence it will contain the true value of the outcome with a certain probability. However, the requirement of validity by itself is vacuous. For example, the prediction set $C$ such that $C(x) =(-\infty, \infty)$, for all $x \in \mathcal{X}$, is technically valid at any miscoverage level, but it is also uninformative. Therefore, our target parameter of interest is an oracle prediction set, denoted by $C^{*}_{\mathcal{P}_I}$, which is defined as the minimal-volume prediction set that satisfies the validity conditions in Definition (ref). For simplicity, we will omit the subscript and write $C^* =C^{*}_{\mathcal{P}_I}$ whenever the dependence on the identified set ${\mathcal{P}_I}$ is clear. Since the minimal-volume prediction set is intuitively the most informative among all valid prediction sets, this property is also called efficiency.\footnote{This notion of efficiency, vovk2009OnlinePredictive,lei2013DistributionfreePrediction,lei2014DistributionfreePrediction differs from the usual definition, where an estimator is considered efficient if it has the smallest variance within a class of estimators. However, the concepts are related, as an estimator with a smaller variance generally leads to the construction of tighter and more informative confidence intervals.}

For interpretability, we will focus on the prediction sets such that $C(x)$ takes the form of a union of a (potentially unknown) number of disjoint intervals for each $x\in \mathcal{X}$. Specifically, we assume that $C(x) \in \mathcal{C}$, $P_X$-almost surely, where

equation*[equation* omitted — 108 chars of source]

denotes the collection of all finite disjoint unions of closed intervals, and the symbol $\sqcup$ indicates that the union is disjoint. The following oracle prediction set under partial identification will be our estimation target.

definition[Oracle prediction Set under Partial Identification] The oracle prediction set under partial identification for a given miscoverage level $\alpha$ is defined as \(C^{*} = C^{*}_{\mathcal{P}_{I}}\subset \mathcal{X}\times \mathbb{R}\) which solves the following problem for $P_{X}$-almost surely \(x\in \mathcal{X}\), \begin{equation} C^*(x) = \arg \min_{C(x)\in \mathcal{C}} \mu(C(x)) \mathtext{s.t.} \inf_{P\in \mathcal{P}_I} P\pqty{Y\in C(x)\mid X = x} \geq 1 - \alpha, \end{equation} where \(\mu\) denotes the Lebesgue measure. Since conditional validity implies marginal validity, the oracle prediction set will also be marginally valid under partial identification.

Before proposing a feasible estimation procedure for the oracle prediction set $C^{*}$, we explain why the existing estimation methods are not directly applicable for prediction sets in the case of interval outcome observations. This analysis highlights the challenges posed by partial identification and the difference between our approach and the classical level set method studied in the literature.

When the joint distribution of $(X,Y)$ is point identified, and $\mathcal{P}_{I} =\cqty{P_{X,Y}}$, it is well known in the literature that the oracle prediction set is related to the upper level set of the conditional density of $Y$ assuming this density is well defined. Let $p(y| x)$ denote the conditional density of $Y$ at $y$ given $X = x$. An upper level set for $p(y\mid x)$ with threshold $\lambda$ at $x \in \mathcal{X}$ is defined as,

equation*[equation* omitted — 68 chars of source]

Let $\lambda_{x}^{\alpha}$ be chosen such that $$\int \indicator{p(y\mid x) \geq \lambda_{x}^{\alpha}} p(y\mid x) \dd y = 1 -\alpha. $$ The oracle prediction set under point identification is $C^{*}_{P_{Y|X}}=\cqty{(x, y): x\in \mathcal{X}, y\in L(x, \lambda_{x}^{\alpha})}$, which is both conditionally and marginally valid (see lei2014DistributionfreePrediction). For a given $x$, it collects the values of $y$ with sufficiently high density and intuitively assembles the set of $y$'s with the highest probability for a given volume. Once this is calibrated to yield the desired miscoverage level, one obtains a minimal-volume set. Estimating the oracle prediction set $C^{*}_{P_{Y|X}}$ therefore reduces to the problem of estimating the upper level sets $L(x, \lambda)$ of the conditional density $p(y \mid x)$ and selecting an appropriate threshold parameter.

However, this characterisation of the oracle prediction set in terms of level sets is no longer applicable under interval censoring, as $P_{Y\mid X}$ is not identified. In fact, the next proposition states that no meaningful bounds can be obtained for the level set of the conditional density of $Y$ given $X$ under interval censoring. For brevity, we temporarily drop the conditioning predictors $X$ in the following proposition. (It can be interpreted as conditional on $X =x$ for $x \in \mathcal{X}$.)

propositionSuppose $Y$ is interval valued, for any $\lambda <\infty$ and any $y\in [a,b]$, where $a<b$ and $(a,b) \in \mathbb{R}^{2}$ is in the support of $(Y^{L}, Y^{U})$, there exists some $P'\in \mathcal{P}_{I}$, such that $y\in L'(\lambda)$, where $L'(\lambda)$ is the upper level set for the density $p'(y)$ of $P'$.

Proposition (ref) states that for any point $y \in [a, b]$ such that $(Y^{L}, Y^{U})$ has positive density at $(a, b) \in \mathbb{R}^2$, we can find a random variable $Y'$ that is observationally equivalent to $Y$ under censoring, and this $Y'$ admits a level set that contains $y$ for any threshold $\lambda$. This is possible because one can adjust the density of $Y'$ at any point $y$ satisfying the stated conditions to be arbitrarily large while ensuring that this density remains within the identified set $P_I$.

The following two lemmas characterise the identified set of the conditional distribution of $Y$ given $X$ under interval censoring. Lemma (ref) relates the sharp identified set of the conditional distribution of $Y$ given $X$ to the {conditional distribution} of $(Y^L, Y^U)$, which can be directly identified and estimated from the data, and provides a feasible lower bound for the restriction in Equation ((ref)). Lemma (ref) shows that the lower bound is tight, and we are not unnecessarily conservative in constructing the oracle prediction set.

lemma[Theorem SIR-2.3 in molinari2020Chapter5] The sharp identified set for the conditional distribution of \(Y\) given \(X\), under interval censoring is given by \begin{equation} \cqty{P_{Y\mid X}: P_{Y\mid X}(Y \in [t_{0},t_{1}]\mid X = x) \geq P\pqty{[Y^{L},Y^U] \subset [t_{0},t_{1}] \mid X = x}, \forall t_{0} \leq t_{1}}. \end{equation}

Lemma (ref) suggests that we could implement a feasible version of Equation (ref), by replacing the optimisation constraint $\inf_{P\in \mathcal{P}_I} P\pqty{Y\in C(x)\mid X = x} \geq 1 - \alpha$ with $$P\pqty{[Y^{L},Y^U] \subset C(x)\mid X = x} \geq 1 - \alpha.$$ This implies the original constraint, since for $C \in \mathcal C$, we can write $C = \sqcup_{m=1}^M [a_{m}, b_{m}]$, and

align*[align* omitted — 250 chars of source]

for a disjoint union of intervals $C\in \mathcal C$. As a result, a prediction set that satisfies the feasible version of the constraint will have a coverage guarantee since

align*[align* omitted — 165 chars of source]

A prediction set $C$ is conservative if one of the three inequalities is strict. The first inequality is due to the partial identification of the joint distribution of $(X,Y)$ and the last inequality is determined by the distribution of $(X, Y^{L}, Y^{U})$. Both are intrinsic properties and cannot be improved without imposing additional restrictions. Lemma (ref) shows that for $C \in \mathcal{C}$, the second inequality is an equality and hence implementing the feasible version of the optimisation problem induces no additional conservativeness.

lemmaFor $P_X$-almost everywhere $x\in \mathcal{X}$, and any $C \in \mathcal{C}$, we have under interval censoring, \begin{equation*} \inf_{P\in \mathcal{P}_I} P\pqty{Y\in C\mid X = x}= P\pqty{[Y^{L},Y^U] \subset C\mid X = x}. \end{equation*}

The above lemma shows that even though the oracle prediction set can be conservative, i.e. overcovering $Y$ for the distribution actually generating the data, the conservativeness is only due to partial identification and the conditional distribution of $(Y^L, Y^U)$ and not because the constraint we focus on differs from the original one.

In fact, given Lemmas (ref) and (ref), we could replace the condition in Equation ((ref)) which depends on the (unobserved) conditional distribution of the latent outcome $Y$ with a condition that depends on the conditional distribution of $(Y^{L}, Y^U)$, which can be estimated given the observed sample $\mathcal{D}$. We then arrive at the following feasible estimation problem:

proposition[Feasible Estimation] The oracle prediction set $C^*$ is the solution to the following feasible optimisation problem for $x$ in the support of $X$, \begin{equation} C^*(x) = \arg \min_{C(x)\in \mathcal{C}} \mu(C(x)) \mathtext{s.t.} P\pqty{[Y^{L},Y^U] \subset C(x)\mid X = x} \geq 1 - \alpha. \end{equation} In the special case when $C^*$ is a single interval, we will denote it as $C^*_{I}$.

A natural estimator $\hat{C}$ for the oracle prediction set can be obtained by solving the sample version of the optimisation problem above. The next two sections study the consistency of such estimators for the oracle prediction set and provide a modification of the estimated prediction set that has finite sample coverage guarantees using conformal inference.

Consistency

We first consider the estimation of the oracle prediction interval $C^*_{I}= [\tau_{0}(x), \tau_{1}(x)]$, for $x\in \mathcal{X}$. This analysis highlights the conditions required to ensure the consistency of the prediction set under different schemes of interval censoring. Subsequently, we will extend our analysis to the estimation of oracle prediction sets $C^* \in \mathcal{C}$ consisting of multiple intervals. To measure the difference between two sets $A, B \subset \mathbb{R}$, we use the volume of their symmetric difference, denoted by $\mu(A \bigtriangleup B)$. The symmetric difference is defined as $A \bigtriangleup B = (A \setminus B) \cup (B \setminus A)$. This is commonly used in the literature on level set estimation and conformal inference. Another commonly used notion of set distance is the Hausdorff distance. In Appendix (ref), we discuss the relationship between these two metrics and why the Hausdorff distance is unsuitable in certain pathological cases.

Let \(\hat{C}_{I}(x)= [\hat{\tau}_{0}(x), \hat{\tau}_{1}(x)]\) denote the solution to the empirical analogue of the optimisation problem defined in (ref) given a random sample $\mathcal{D}$

equation[equation omitted — 122 chars of source]

Here we use the shorthand notation \(P(t_{0},t_{1};x) = P(t_{0}\leq Y^{L} \leq Y^{U} \leq t_{1}\mid X = x)\), and \(P_{n}(t_{0},t_{1};x)\) is an estimator for \(P(t_{0}, t_{1};x)\) based on the random sample \(\mathcal{D} = \{(X_{i}, Y_{i}^{L}, Y_{i}^{U}): i\in \mathcal{I}\}\) with sample size \(n\).

assumption[Estimation] The estimator \(P_{n}(t_{0},t_{1};x)\) is consistent for \(P(t_{0},t_{1};x)\) uniformly over \((t_{0}, t_{1}, x)\in \mathbb{R}^{2}\times \mathcal{X}_{n}\), for some \(\mathcal{X}_{n}\subset \mathcal{X}\), \begin{equation*} \sup_{x\in \mathcal{X}_{n}}\sup_{t_{0}\leq t_{1}} \abs{P_{n}(t_{0},t_{1};x) - P(t_{0},t_{1}; x)} = o_{p}(1). \end{equation*}

This assumption does not specify the estimator \(P_{n}\) we use. For example, a candidate estimator \(P_{n}(t_{0},t_{1};x)\) based on kernel smoothing is

equation[equation omitted — 211 chars of source]

where \(K_{h}(x) = \frac{1}{h} K\pqty{\frac{x}{h}}\) for a choice of kernel smoothing function \(K(\cdot)\) and bandwidth parameter \(h\), see tsybakov2009IntroductionNonparametric and hansen2008UniformConvergence. In this case, $\mathcal{X}_{n}$ can be an expanding subset of $\mathcal{X}$ at a suitable rate to achieve uniform consistency as in hansen2008UniformConvergence, Theorem 8 (see discussion in Appendix (ref)). Another candidate is the distributional regression estimator foresi1995ConditionalDistribution, chernozhukov2013InferenceCounterfactual,chernozhukov2021DistributionalConformal.

Given a miscoverage level $\alpha \in (0,1)$, the next assumption we make is that \(C_{I}^*(x) = [\tau_{0}(x), \tau_{1}(x)]\) is the unique interval that satisfies the $1-\alpha$ coverage condition with the shortest length.

assumption[Identification] There exists a unique solution \(C^*_{I}(x)=[\tau_{0}(x), \tau_{1}(x)]\) to the optimisation problem in (ref): that is, \(P(\tau_{0}(x),\tau_{1}(x); x) \geq 1 -\alpha\), and, for any \((t_{0}, t_{1}) \neq (\tau_{0}(x), \tau_{1}(x))\) such that \(t_{1}- t_{0} \leq \tau_{1}(x) - \tau_{0}(x)\), then \(P(t_{0},t_{1}; x) < 1 - \alpha\). Additionally, for any \(\epsilon > 0\), there exists \(\delta > 0\) such that for all \(x\), if \(t_{1}(x) - t_{0}(x) \leq \tau_{1}(x) - \tau_{0}(x)\) and \(\mu([t_{0}(x), t_{1}(x)] \bigtriangleup [\tau_{0}(x), \tau_{1}(x)]) > \epsilon\), then $P(t_{0}, t_{1}; x) < P(\tau_{0}, \tau_{1}; x) - \delta$.

Assumption (ref) states that we cannot achieve the same coverage probability with a different interval $(t_{0}(x),t_{1}(x))$ that is not larger than $(\tau_{0}(x), \tau_{1}(x))$. Figure (ref) shows a specific joint density of \(Y^{L}\) and \(Y^{U}\). The set \(\{(a,b): t_{0}\leq a \leq b \leq t_{1}\}\) is a triangular region in \(\mathbb{R}^{2}\). In Figure (ref), the plotted triangular region corresponds to $t_{0}= -0.75$ and $t_{1}= 0.5$. $P(t_{0},t_{1};x)$ will be the integral of the joint density over the triangular region. The first part of Assumption (ref) thus requires that if either the triangular region is smaller or if we move the triangular region away from the optimal position, we will strictly lose coverage probability of \((Y^{L}, Y^U)\); and the second part is a regularity condition ensuring that the “coverage loss” is bounded away from $0$ over the support of $X$. This rules out the case when \(P(t_{0},t_{1};x)\) gets flatter around \((\tau_{0}(x),\tau_{1}(x))\) for some \(x\).

figure[figure omitted — 342 chars of source]

In the absence of fixed censoring, we assume the following regularity condition for the conditional distribution \(P(t_{0},t_{1};x)\). The case with fixed censoring will be discussed later.

assumption[Regularity] There exists an \(\epsilon_{0}>0\) and positive constants \(c_{1}, \gamma_{1}\), such that for all \(x\in \mathcal{X}\) and \(\tau_{0}(x),\tau_{1}(x)\) as defined in Assumption (ref), \begin{equation*} P\pqty{\tau_{0}(x)- \epsilon, \tau_{1}(x)+ \epsilon; x} - P\pqty{\tau_{0}(x), \tau_{1}(x);x} \geq c_{1}\epsilon^{\gamma}, \end{equation*} for any \(0 < \epsilon < \epsilon_{0}\). In addition, for all \(\epsilon < \epsilon_{0}\) and \(t_{0}\leq t_{1}\), \(P(t_{0}- \epsilon,t_{1}+ \epsilon;x) - P(t_{0},t_{1};x)\leq c_{2}\epsilon^{\gamma_{2}}\) for some constant \(c_{2}, \gamma_{2}\).

Assumption (ref) is a smoothness assumption on \(P(t_{0},t_{1};x)\). In particular, we are assuming that by considering a slightly larger interval \([\tau_{0}- \epsilon, \tau_{1}+ \epsilon]\) we can strictly increase our coverage probability. Assumption (ref) also assumes an upper bound on the density, which can be violated if the conditional distribution of \((Y^{L}, Y^{U})\) is discrete, which might happen under a fixed censoring scheme. This assumption is related to the \(\gamma\)-exponent condition in the level set literature, see polonik1995MeasuringMass and lei2014DistributionfreePrediction.

The following theorem shows that the estimated prediction interval is consistent for the oracle prediction interval under random censoring.

theoremUnder Assumptions (ref), (ref), and (ref), the prediction interval estimator \(\hat{C}_{I}(x)\) defined in Equation (ref) is uniformly consistent for the oracle prediction interval \(C_{I}^{*}(x)\), that is, \begin{equation*} \sup_{x\in \mathcal{X}_{n}}{\mu\pqty{\hat{C}_{I}(x)\bigtriangleup C_{I}^{*}(x)} } = o_{p}(1) \end{equation*} as the sample size \(n\to \infty\).

Assumption (ref) can be violated when there is fixed censoring, i.e., when there are some fixed brackets, and we observe only the bracket in which an outcome falls. Suppose the conditional distribution of $(Y^L, Y^U)$ can be decomposed into a continuous part and a discrete part, where the discrete part has support over the set of points \(\{(L_{k}, U_{k})\in \mathbb{R}^2: 1 \leq k \leq K\}\).

In the case of a mixture of random and fixed censoring, the regularity condition in Assumption (ref) is violated. We impose the following assumption, which slightly modifies the estimator to allow for non-smoothness of the conditional distribution.

assumption[Approximate solution] The estimator for the oracle prediction interval \(\hat{C}_{I}(x) = [\hat{\tau}_{0}(x), \hat{\tau}_{1}(x)]\) is constructed by solving the following optimization problem, for $x \in \mathcal{X}$, \begin{equation} \min_{t_{0} < t_{1}} \, \, t_{1}- t_{0} \mathtext{s.t.} P_{n}(t_{0}, t_{1}; x) \geq 1 - \alpha - \psi_{n}, \end{equation} for some $\psi_{n}>0$ and $\psi_{n}\to 0$, such that for some $\mathcal X_n$, \begin{equation} P\pqty{\sup_{x\in \mathcal{X}_n} \sup_{t_{0}\leq t_{1}}\abs{P_{n}(t_{0},t_{1};x) - P(t_{0},t_{1};x)} > \psi_{n}}=o(1). \end{equation}

To see the reason why fixed censoring creates a problem, and we should not solve the exact constraint, consider the following example.

exampleSuppose $(Y_i^L, Y_i^U) = (0,1)$ or $(0,2)$ with equal probability without predictor information $X$. The oracle prediction set for a miscoverage level $\alpha = 0.5$ (see Equation (ref)) is $C^* = [0,1]$. (Notice that $\inf_{P\in \mathcal{P}_I} P\pqty{Y\in C^*} = 0.5$ and a shorter interval would not achieve the desired miscoverage level.) If we solve the exact empirical analogue of the oracle prediction set with the natural estimator $P_n(a,b) = n^{-1} \sum \indicator{[Y_i^L, Y_i^U] \subset [a,b]}$ as in (ref), we obtain a prediction interval $$\check{C} = \bqty{0, 1 + \indicator{\bar{Y}^U > 1.5}},$$ where $\bar{Y}^U = \frac{1}{n}\sum_i Y_i^U$. Since the lower bound for $(Y_L, Y_U)$ is 0 with probability one, it is easy to see that the lower bracket of $\check{C}$ is $0$. Because $$P_n(0,b) = \begin{cases} 0 &\mathtext{if} 0 < b < 1 \\ n^{-1} \sum \indicator{Y_i^U =1} &\mathtext{if} 1\leq b < 2 \\ 1 &\mathtext{if} b\geq 2 \end{cases},$$ the upper bracket of $\check{C}$ is $2$ if $P_n(0,1) < 0.5$, i.e. there are fewer observations of $Y_i^U =1$ than $Y_i^U = 2$ in the sample, which is equivalent to having the average $\bar{Y}^U > 1.5$. Since $\bar{Y}^U$ has a scaled Binomial distribution $\bar{Y}^U \sim 1 + n^{-1} \text{Binom}(n,0.5)$, and the median for a $\text{Binom}(n,p)$ distribution is between $\lfloor np \rfloor$ and $\lceil np \rceil$, we have that $P(\bar{Y}^U>1.5) \to 0.5$ as $n \to \infty$. This means that $P(\check{C} = [0,2]) \to 0.5$ as $n\to \infty$ and \(\check{C}\) is not consistent. On the other hand, for any $0 < \psi_n < 0.5$, $P_n(0,1) < 0.5-\psi_n$ is equivalent to $\bar{Y}^U > 1.5 + \psi_n$, which has probability approaching $0$ by the central limit theorem if $\sqrt{n}\psi_n \to \infty$ as the sample size $n \to \infty$. Solving the modified optimisation problem in (ref) will then lead to a consistent estimator of the oracle prediction set in this example. Also note that the condition in (ref) is a generalisation of the condition $\sqrt{n}\psi_n \to \infty$.

With (ref) in place of Assumptions (ref) and (ref), we obtain the following result:

theoremUnder Assumptions (ref), and (ref), we have the uniform consistency of $\hat{C}_{I}$ for the oracle prediction interval $C_{I}^{*}$, that is, as the sample size $n\to \infty$, \begin{equation*} \sup_{x\in \mathcal{X}_n} \mu\pqty{\hat{C}_{I}(x)\bigtriangleup C^{*}_{I}(x)} = o_{p}(1). \end{equation*}

We now turn to the more general case where the estimation target is the {oracle prediction set} $C^*\in \mathcal{C}$ defined in (ref).

The next assumption is on the identification of the oracle prediction set, which generalises Assumption (ref). It additionally assumes that the number of disjoint intervals in the oracle prediction set is bounded.

assumption\begin{enumerate} • There exists a unique solution $C^*$ to (ref). For any $\epsilon > 0$, there exists $\delta = \delta(\epsilon) > 0$, such that for $P_{X}$-almost surely $x\in \mathcal{X}$ and any set ${C}'(x) \in \mathcal{C}$ with $\mu({C}'(x)) \leq \mu(C^{*}(x))$ and $\mu(C^{*}(x) \triangle C'(x)) > \epsilon$, we have $P([Y^L, Y^U] \subset {C}'(x)\mid X = x ) < 1 - \alpha-\delta$. • The oracle \(C^*(x) \in \mathcal{C}_{\bar{M}}\), where $\mathcal{C}_{\bar{M}} = \{\sqcup_{m=1}^{M} [a_m,b_m]: a_m < b_m < a_{m+1}, \, M \leq \bar{M} < \infty\}$, for some $\bar{M} < \infty$ and almost surely $x \in \mathcal X$. \end{enumerate}

With the notation $P(C(x);x) = P([Y^{L}, Y^{U}]\subset C(x) \mid X = x)$, and suppose that we have a uniform consistent estimator $P_n(C; x)$. We define the estimator $\hat{C}(x), x\in \mathcal{X}$, for the oracle prediction set $C^*$, allowing for potential fixed censoring.

assumptionThe estimators $\hat{{C}}(x), x\in \mathcal{X}$ are constructed by solving the following problem, for any $x\in \mathcal{X}$, and $M \geq \bar{M}$, \begin{equation*} \min_{C(x)\in \mathcal{C}_M} \mu(C(x)) \mathtext{s.t.} P_{n}(C(x); x) \geq 1 - \alpha - \psi_{n}. \end{equation*} Here $\psi_{n}=o(1)$ is a sequence of positive numbers and $P_{n}$ is a consistent estimator of $P$ such that $$P\pqty{{\sup_{C\in \mathcal{C}_M}}\sup_{x\in \mathcal{X}_n} \abs{P_{n}(C(x);x) - P(C(x);x)} > \psi_{n}}=o(1),$$ as $n\to \infty$.
remark(ref) states conditions under which a kernel smoothing estimator $P_n(C;x)$ is uniformly consistent and provides the corresponding convergence rate. This convergence rate can be used to inform the choice of $\psi_n$. Empirically, varying $\psi_n$ between 0 and small positive values produces negligible changes in the results. The uniformity is restricted to the class $\mathcal{C}_M$ rather than the larger class $\mathcal{C}$ because the proof relies on arguments that are closely tied to $\mathcal{C}_M$ being a VC class of sets. Notably, knowledge of $\bar{M}$ is not required, as it suffices to select an $M$ that is sufficiently large. Even if $M$ is smaller than the true $\bar{M}$, the resulting prediction set can be viewed as an approximation for the oracle prediction set.
theoremUnder Assumptions (ref) and (ref), we have that $\hat{C}$ is a consistent estimator of ${C}^*$ in the sense that, as $n\to \infty$, \begin{equation*} \sup_{x\in \mathcal{X}_n} {\mu\pqty{\hat{C}(x) \bigtriangleup C^*(x)}} = o_{P}(1). \end{equation*}

We have proposed consistent estimators for the oracle prediction set. There has been a recent interest in constructing prediction sets that are finite-sample valid via conformal inference. The next section discusses how to construct a conformal prediction set that is finite-sample valid under interval censoring.

Conformal inference with interval outcomes

Given a consistent estimator for the oracle prediction set \(\hat{C}\), this section considers the problem of constructing a prediction set with finite-sample validity. An estimated prediction set \(\tilde{C}\) with a random sample $\mathcal{D} = \cqty{(X_{i}, Y^{L}_{i}, Y^{U}_{i}), i \in \mathcal{I}}$ is said to be finite-sample marginally valid, if for all \(P\),

equation*[equation* omitted — 77 chars of source]

where \(\Prob\) denotes the joint probability of \(\cqty{\bar{Z}_{i} = (X_{i}, Y_{i}, Y^L_{i}, Y_{i}^{U}): i \in \mathcal{I}\cup\cqty{{n+1}}}\). Notice that this is the finite-sample validity condition typically considered in conformal inference literature, where the probability is the joint distribution for both the observed sample and the individual to be predicted. While the oracle prediction set \(C^*(x)\) depends on the true distribution \(P\), conformal inference constructions lead to a prediction set that is finite-sample valid for all \(P\), making it distribution-free.

In particular, we will achieve finite-sample validity using the method of split conformal inference (see romano2019ConformalizedQuantile and lei2018DistributionfreePredictive). One feasible split conformal inference procedure for the oracle prediction interval $C^*_I$ in our setting can be described as follows. We first split the sample \(\mathcal{D}\) into a training set \(\mathcal{D}_{1} = \cqty{Z_{i}: i\in \mathcal{I}_{1}}\) and a calibration set \(\mathcal{D}_{2}= \cqty{Z_{i}: i \in \mathcal{I}_{2}}\). We calculate the estimator \(\hat{C}_{I}(x) = [\hat{\tau}_{0}(x), \hat{\tau}_{1}(x)]\) for the oracle prediction set with the training set for \(x\in \mathcal{X}\) as in (ref). The following scores \(s_{j} = s(Y^{L}_{j}, Y^{U}_{j}, X_{j}; \hat{\tau}_{0}, \hat{\tau}_{1})\) are then computed for each \(j \in \mathcal{I}_{2}\) in the calibration set, given \(\hat{C}_{I}(x) = [\hat{\tau}_{0}(x), \hat{\tau}_{1}(x)]\), where

equation*[equation* omitted — 117 chars of source]

This score function is a modified version of the score function used in the case of conformalized quantile regression with no interval censoring studied in romano2019ConformalizedQuantile. The score function is negative when the interval \([y^L, y^U]\) is contained in \([t_{0}(x),t_{1}(x)]\) and positive otherwise. By allowing for negative scores, the conformal prediction set \(\tilde{C}_{I}(x)\) penalises both overcoverage and undercoverage.

We can then compute the quantile $$\vartheta_{1 - \alpha} = q\pqty{\cqty{s_{j}: j \in \mathcal{I}_{2}}; (1 - \alpha)(1+ 1/\abs{\mathcal{I}_2})},$$ where \(q\pqty{\mathcal{S}; \zeta}\) denotes the \(\zeta\)-quantile of the set \(\mathcal{S}\) of real numbers. The conformal prediction set is then

equation[equation omitted — 208 chars of source]

Intuitively, when our estimated prediction set \(\hat{C}_{I}\) is wider than necessary, the intervals \(\{(Y_{j}^{L}, Y_{j}^{U}), j \in \mathcal{I}_{2}\}\) tend to be contained within \(\hat{C}_{I}\), leading to a higher proportion of negative scores in the calibration set. As a result, the chosen quantile \(\vartheta_{1 - \alpha}\) is negative, and the conformal prediction interval adjusts for the over-coverage of \(\hat{C}_I\). The correction for under-coverage follows the same logic.

When the prediction set $\hat{C}$ is a union of multiple intervals, we can construct a finite-sample valid conformal prediction set as follows. First we split the samples into a training set $\mathcal{I}_{1}$ and a calibration set $\mathcal{I}_{2}$. We estimate $\hat{C}\in \mathcal{C}$ using data in the training set, and then compute the following conformity score$s_{j}= s(Y^{L}_{j}, Y^{U}_{j}, X_{j}; \hat{C})$, for each $j\in \mathcal{I}_{2}$. In order to allow for $\hat{C}$, which is a union of intervals, we compute the conformity score as follows. For $C(x) = \sqcup_{m=1}^{M} [t_{0m}(x), t_{1m}(x)] \in \mathcal{C}$, where $t_{0m} < t_{1m} < t_{0,m+1}$,

equation[equation omitted — 159 chars of source]

This score can be computed explicitly as follows. With the notation $t_{10} = -\infty, t_{0,M + 1} = \infty$, if there exists $j: 1\leq j \leq M$ such that $t_{0,j} \leq y^L \leq y^U \leq t_{1, j}$, then $s(y^L, y^U; C) = s(y^L, y^U; t_{0j}, t_{1j})$, otherwise, for $i = \max\{j: t_{1,j - 1}\leq y^L\}$, and $k = \min\{j : t_{0,j + 1} \geq y^u\}$,

align*[align* omitted — 273 chars of source]

Notice that $s_{j} \leq 0$ if and only if $[Y^{L}_{j}, Y^{U}_{j}] \subset [t_{0m}, t_{1m}]$ for some $m\leq M$. We can then compute $\vartheta_{1 - \alpha} = q(\{s_{j}: j \in \mathcal{I}_{2}\}; (1 - \alpha)(1+ 1/\abs{\mathcal{I}_2}))$. The conformal prediction set is then given by \footnote{Here, the endpoints of the estimated $\hat{C}$ are adjusted uniformly using the same $\vartheta_{1-\alpha}$. A more sophisticated adjustment method is discussed in (ref), where we propose a minimal Euclidean distance adjustment. We conjecture that this approach could improve the convergence rate of the conformal prediction set. The optimal choice of conformity score remains an interesting question, see for example, thurin2025OptimalTransportbased.}

equation[equation omitted — 165 chars of source]

The next theorem states the finite-sample marginal validity of the conformal prediction set \(\tilde{C}\) defined in (ref) above. It implies the validity of the conformal prediction interval $\tilde{C}_I$.

theoremUnder the assumption that \(\{(X_{i}, Y_{i}, Y_{i}^{L}, Y_{i}^{U}): i\in \mathcal{I}\cup \{n+1\}\}\) are i.i.d.\footnote{In order to perform conformal inference, it is sufficient that the sample satisfies the weaker condition of exchangeability. Here, independence is assumed for simplicity and for establishing the asymptotics for the nonparametric estimators we will use. barber2023ConformalPrediction studies conformal inference beyond exchangeability and is applicable for time series data. Our results can be readily extended to the time series settings combined with their procedure. }, the conformal prediction set $\tilde{C}$ in (ref) is finite-sample marginally valid, that is, for all \(P\), \begin{equation*} \Prob\pqty{Y_{n+1}\in \tilde{C}(X_{n+1})} \geq 1 - \alpha. \end{equation*} where $\Prob$ denotes the joint probability of \(\cqty{\bar{Z}_{i}: i \in \mathcal{I}\cup\cqty{{n+1}}}\).

The following theorem shows that the conformal prediction set \(\tilde{C}(x)\) preserves the consistency of the estimator \(\hat{C}(x)\) for the oracle prediction set. Recall the definition of $s(y^{L}, y^{U}, x; C)$ in Equation (ref), let \(s_{j} = s(Y^{L}_{j}, Y^{U}_{j}, X_{j}; \hat{C})\), and \(s_{j}^{*} = s\pqty{Y_{j}^{L}, Y_{j}^{U}, X_{j}; C^*}\).

theoremGiven an i.i.d. sample \(\mathcal{D} = \{(X_{i}, Y_{i}^{L}, Y^{U}_{i} ): i\in \mathcal{I}\}\) which is split into training \(\mathcal{I}_{1}\) and calibration sets \(\mathcal{I}_{2}\). Let $\tilde{C}(x)$ be the conformal prediction set constructed by the preceding procedure in (ref), if $s(Y^L, Y^U, X; C^*)$ has positive density in a neighbourhood around $0$, then under the assumptions of (ref), \begin{equation*} \sup_{x} \mu\pqty{{\tilde{C}(x) \bigtriangleup C^*(x)}} = o_{p}(1) \end{equation*} as the sample sizes \(n_{1}, n_{2}\to \infty\), here \(n_{1}= \abs{\mathcal{I}_{1}}\) and \(n_{2}= \abs{\mathcal{I}_{2}}\).

It is tempting to look for a finite-sample conditionally valid prediction set \(\tilde{C}(x)\) such that for all \(P\),

equation*[equation* omitted — 88 chars of source]

However, it is impossible to construct a prediction set that is finite-sample conditionally valid and consistent for the oracle prediction set lei2014DistributionfreePrediction, barber2021LimitsDistributionfree, gibbs2024ConformalPrediction, and hence we consider the weaker notion of local validity. Let \(\mathcal{A}= \cqty{A_{k}: k\leq K}\) be a partition of \(\mathcal{X}\), a prediction interval \({C}_{\text{loc}}(x)\) is finite-sample locally valid with respect to \(\mathcal{A}\), if for all \(P\),

equation*[equation* omitted — 117 chars of source]

Let \(\mathcal{I}_{2,k} = \cqty{i\in \mathcal{I}_{2}: X_{i} \in A_{k}}\), \(n_{2,k} = \abs{\mathcal{I}_{2,k}}\), and $$\vartheta_{1-\alpha,k} = q \pqty{\cqty{s_{j}: j \in \mathcal{I}_{2,k}}; (1 - \alpha)\pqty{1+\frac{1}{n_{2,k}}}},$$ the local conformal prediction set is

equation[equation omitted — 269 chars of source]

The finite-sample local validity is a result of exchangeability, similar to lei2014DistributionfreePrediction.

theoremSuppose \(\{(X_{i}, Y_{i}^{L}, Y_{i}^{U}): i\in \mathcal{I}\}\) are i.i.d. Let \(\mathcal{A}\) be a partition of \(\mathcal{X}\). Then the local conformal prediction set \(\tilde{C}_{\text{loc}}(x)\) is finite-sample locally valid with respect to \(\mathcal{A}\).

In the next section, we demonstrate the performance of the conformal and local conformal prediction sets in a set of numerical experiments.

Numerical experiments

We report here the results of a series of numerical experiments to evaluate the performance of the conformal prediction sets. In all the experiments, we randomly draw a sample of predictors \(\cqty{X_{i}, i= 1,2,\dots, n}\), with sample size \(n=2,500\) and each \(X_{i}\sim \text{Unif}[-1.5, 1.5]\). The miscoverage level is fixed at \(\alpha = 0.1\). The unobserved outcomes \(Y_{i}\) and the observed intervals \((Y_{i}^{L}, Y_{i}^{U})\) are then generated based on one of the models described below.

itemize• Model A. The unobserved outcomes \(Y_i\) are generated according to a mixture model with heteroscedasticity, similar to the simulations in lei2014DistributionfreePrediction. \begin{equation} (Y_{i} \mid X_{i} = x)\sim B_{i} N(f(x)+ g(x), \sigma^{2}(x)) + (1-B_{i}) N(f(x)- g(x), \sigma^{2}(x)) \end{equation} Here, the auxiliary random variable \(B_i\) is independently drawn from a Bernoulli distribution with probability \(0.5\), and the following functions are used. \begin{align*} f(x) = 2(x-1)^{2}(x+1);\, g(x) = 4 \sqrt{(x + 0.5) \indicator{x\geq -0.5}};\, \sigma^{2}(x) = \frac{1}{4} + \abs{x}, \end{align*} The observed intervals \((Y_i^L, Y_i^U)\) are obtained by adding independent noise to the unobserved outcomes. \begin{equation} Y_{i}^{L} = Y_{i} - e_{1,i}, \quad Y_{i}^{U} = Y_{i} + e_{2,i}, \quad e_{1,i}, e_{2,i} \sim_{i.i.d.} |N(0,1)|. \end{equation} In this design, the conditional distribution \(P_{Y\mid X}\) is unimodal for smaller values of \(X\) and bi-modal for larger values of \(X\). Regression and quantile regression-based prediction sets will have difficulty capturing the bifurcation shape of the oracle prediction set. • Model B. The unobserved outcomes are generated in the same way as in (ref) above. We randomly select 20% of the unobserved outcomes, and they are interval valued, and in this case and \((Y^{L}_{i}, Y_{i}^{U})\) have a discrete distribution, while the remaining 80% are observed without censoring. Let \(\iota_{i}\) be a random draw from a Bernoulli distribution with probability 0.2, \begin{equation} \begin{cases} Y_{i}^{L} = Y_{i}^{U} = Y_{i} &\mathtext{if} \iota_{i} = 0;\\ Y_{i}^{L} = \lfloor Y_{i} \rfloor, Y_{i}^{U} = Y_{i}^{L} +1 &\mathtext{if} \iota_{i} = 1. \end{cases} \end{equation} This model is designed to evaluate the performance of the conformalized prediction set under fixed censoring. We choose to censor 20% of the outcomes, as approximately 20% of the income data is interval valued in the dataset used for our empirical studies. • Model C. The unobserved outcomes are generated according to the following model, \begin{equation*} Y_{i} = f(X_{i}) + \epsilon_{i}, \quad \epsilon_{i} \sim_{i.i.d.} \chi^{2}_{1.5}, \end{equation*} and \((Y_{i}^{L}, Y_{i}^{U})\) are generated according to (ref). In this case, the conditional distribution is unimodal and skewed.

The observations \(\cqty{(X_{i}, Y_{i}^{L}, Y_{i}^{U}), i = 1,\dots,n}\), are randomly split into a training set $\mathcal{I}_{1}$ and a calibration set $\mathcal{I}_{2}$, with the training set containing 75% of the observations. We estimate the prediction set \(\hat{C}(x)\) based on the kernel smoothing estimator \(P_{n}(C;x)\) with the Epanechnikov kernel for the conditional distribution $P(C; x)$, as described in Equation ((ref)), using the observations in the training set. The conformal prediction sets \(\tilde{C}\) are then constructed according to Equation ((ref)) using the calibration set. Local conformal prediction sets \(\tilde{C}_{\text{loc}}\) are obtained using the local conformalisation procedure, where the support \(\mathcal{X}\) is divided into 5 equal-sized bins.

We then compare the average coverage and volume of the conformal prediction set \(\tilde{C}\) and the local conformal prediction set \(\tilde{C}_{\text{loc}}\) over 100 repetitions with two alternative conformal prediction intervals: \(\tilde{C}_{q2}\), based on a quantile regression assuming a quadratic model, and \(\tilde{C}_{q3}\), based on a quantile regression assuming a cubic model. See romano2019ConformalizedQuantile and (ref) for more details about the construction of a prediction interval based on quantile methods. The coverage is computed by randomly generating 5,000 samples of \((Y^{L}, Y^{U})\) and calculating the proportion of intervals covered by the prediction sets. The volume is calculated as the integrated Lebesgue measure of the prediction sets over the support \(\mathcal{X}\).

Figure (ref) shows the comparison of the prediction sets from the numerical experiments. The prediction sets \(\tilde{C}\) and \(\tilde{C}_{\text{loc}}\) capture the bifurcation region of the conditional distribution in Model A and Model B and have a smaller volume compared to the quantile regression-based prediction sets. Even when the conditional distribution is unimodal, and the oracle prediction set is a single interval for each $x\in \mathcal X$, our proposed estimators still have smaller volume as demonstrated in Model C. In terms of coverage, \(\tilde{C}, \tilde{C}_{\text{loc}}\) and \(\tilde{C}_{q3}\) appear to have coverage close to the nominal level $1-\alpha = 0.9$. Notably, \(\tilde{C}_{\text{loc}}\) appears to have better coverage than \(\tilde{C}\) near the boundaries of the support \(\mathcal{X}\).

figure[figure omitted — 2,504 chars of source]

Empirical study

UK Job Advert Data

We have also applied our method to the UK job advert dataset from Adzuna, an online job advert aggregator in the UK. Adzuna data have been used widely as an indicator of economic activity in the UK, including by the Office for National Statistics as experimental statistics ons2021UsingAdzuna. It has also been used to develop a measure of labour market opportunities for heterogeneous types of workers dias2021WorkerMobility.

The original sample consists of job adverts posted in the United Kingdom between January and March of 2022, which comprise 121,746 observations. Each advertisement contains information on job location, job category, and salary, when available. Salary information is missing for 37.8% (46,077) of the postings. Among the advertisements that report salaries, the majority are listed as salary intervals, accounting for approximately 62% (46,606) of all postings with salary information, making this dataset well-suited for our interval-valued outcome framework.

The outcome of interest is the annual salary associated with each job advert, with $(Y_i^L, Y_i^U)$ denoting the reported lower and upper salary bounds for the $i$-th advert, and we work with the logarithm of the salary bounds. In the analysis below, we focus on the IT industry, which is the largest job category by number of postings (see (ref)). For the 19,385 job advertisements in the IT industry, 40% (7,742 postings) report salary information, of which 64% (4,981 postings) provide salaries in intervals; in our analysis, we restrict attention to the subset of postings that report salary information. Job location is represented by geographic coordinates (latitude and longitude), which we treat as continuous covariates.

For a given location (latitude and longitude), we estimate the 90% oracle prediction interval (see (ref)) for the salary distribution in the IT category conditional on location. We employ the kernel smoothing estimator of the conditional probability using a rectangular kernel applied to geographic distance and solve for the empirical analogue of the feasible optimization problem in Equation (ref), yielding prediction interval $\hat{C} = (\hat{\tau}_0(x), \hat{\tau}_1(x))$. Then we construct the local conformal prediction set $\tilde{C}_{\text{local}}$ according to (ref). To visualise spatial variation in predicted salary distributions, we evaluate the predicted upper and lower bounds on a grid of locations covering the UK and the resulting prediction bounds are then interpolated. Importantly, our approach accommodates the interval-valued nature of the salary data without imposing parametric assumptions or relying on imputation of missing salaries.

Figure (ref) plots the estimated lower and upper bounds of the 90% local conformal prediction interval for annual salaries across geographic locations in the UK. For reference, the median annual pay for full-time employees was £33,000 for the tax year ending on 5 April 2022.\footnote{Source: \url{https://www.ons.gov.uk/employmentandlabourmarket/peopleinwork/earningsandworkinghours/bulletins/annualsurveyofhoursandearnings/2022}. Employee earnings in the UK: 2022, Measures of employee earnings using data from the Annual Survey of Hours and Earnings (ASHE), Office for National Statistics.} We can observe a few patterns in the spatial distribution of the predicted salary bounds. First, the estimated lower bounds exhibit moderate geographic heterogeneity and are relatively even across space and are typically below £30,000. Lower bounds tend to be higher around major cities, while rural areas display lower predicted lower bounds or no job postings. Second, the upper bounds display a much thicker right tail and substantially stronger spatial concentration, with predicted values ranging roughly from £50,000 to £230,000. Areas surrounding major metropolitan centres, such as London, Manchester, Bristol, Birmingham and Edinburgh, stand out with markedly higher upper bounds than the rest of the country, indicating pronounced geographic concentration of high-paying jobs. As a result, the predicted salary interval is widest in London and other large cities, reflecting greater wage dispersion and job heterogeneity in these labour markets.

Table (ref) reports the estimated 90% local conformal prediction intervals for annual IT salaries at selected major UK cities, evaluated at their geographic coordinates. The results illustrate substantial cross-city heterogeneity, particularly in the upper bounds of the predicted salary distributions. London stands out with both relatively higher predicted lower bound and upper bound, leading to the widest prediction interval among the listed cities, consistent with high wage dispersion and job heterogeneity in the capital. Other large metropolitan areas such as Manchester, Edinburgh, and Glasgow also exhibit wide intervals driven primarily by elevated upper bounds, while Birmingham displays a comparatively narrower interval, reflecting more moderate dispersion in predicted salaries. Across all cities, predicted lower bounds are relatively similar and generally below \pounds25,000, whereas differences in upper bounds account for most of the variation in interval widths. We also note that, although the salary intervals reported in individual job postings are typically narrow, their levels vary substantially across postings even within the IT industry. This dispersion leads to the relatively large estimated interval widths.

Overall, these findings highlight the importance of accounting for interval-valued salary information when studying spatial variation in labour market opportunities. In particular, these patterns reinforce the importance of analysing the full conditional salary distribution rather than relying on point estimates, as the upper tail of the distribution varies markedly across locations. Our approach captures spatial variation in both the predicted lower and upper bounds, thereby providing richer insights than simple point estimates such as conditional means or medians, without imposing parametric assumptions or relying on imputation of missing salaries.

table[table omitted — 686 chars of source]
figure[figure omitted — 175 chars of source]
figure[figure omitted — 651 chars of source]

US CPS data

We apply our method to the 2023 US Current Population Survey (CPS) dataset and construct a prediction interval for their income given years of education and experience. Following heckman2008EarningsFunctions, we construct the following variables from the dataset for individuals who report being American citizens and employed:

itemize• Income: We include all individuals who reported their income or income range in the survey. Since 2014, when a person refuses to report their income, the interviewer will ask for a range, and the person's income will be imputed by matching with similar people who have income in the same range. In the 2019 ASEC Updates\footnote{See \href{https://cps.ipums.org/cps/asec_2019_changes.shtml}{https://cps.ipums.org/cps/asec_2019_changes.shtml}}, it is stated that “if respondents did not give a value for income received from a given source, the interviewer now follows up with a question about income ranges. Respondents who did not give a specific value were asked if they received over \$60,000, between \$45,000 and \$60,000, or less than \$45,000 from this income source. Those who replied that they earned less than \$45,000 from this source were further asked if they earned more than \$30,000, between \$15,000 and \$30,000, or less than \$15,000 from this income source.” For people who reported their income as larger than $ \$60,000$, we set the upper bound of the income as $\$400,000$, which is the topcode of the income variable in the dataset. The individuals who did not report their income or their income range are excluded. • Age: We include individuals who are between 18 and 65 years old. • Education: Following heckman2008EarningsFunctions, we construct the years of education for each individual by converting the highest qualification they have obtained into years of education. • Experience: We construct the years of experience for each individual as their age minus the years of education minus 6 (which is presumed to be the entry age into education).

The CPS dataset consists of 146,133 observations. There are 53,171 individuals who satisfy the demographic restrictions above. Among the individuals that satisfy the demographic criteria, 39,333 reported their income either exactly or as a band and the response rate is thus around 74%. The issue of nonresponse in surveys has been studied in the literature; see for example, lillard1986WhatWe and bollinger2019TroubleTails, where they document an increasing proportion of nonresponse of earnings in survey data and the potential bias it introduces. We will not address the nonresponse issue in this paper, but it is an important consideration in practice. Our proposed predictions could assist in the imputation of missing values in a robust manner; see lei2021ConformalInference. In fact, allowing respondents to report their income in bands could help mitigate the non-response problem. For the 39,333 individuals who reported their incomes and satisfy the demographic restrictions, 8,819 of them reported either their gross or net income as a band, which thus amounts to roughly 22.4%. of the respondents. Figure (ref) shows the proportion of individuals who reported their income in a range by years of education.

figure[figure omitted — 247 chars of source]

In our empirical study, we first randomly reserve 20% of the data as a hold-out dataset. The remaining 80% is then further divided into training and calibration sets, with 75% of the data used for training to estimate the prediction interval \(\hat{C}_{I}\) and 25% used for constructing the conformalized prediction interval \(\tilde{C}_{I}\). Table (ref) presents the average estimated prediction intervals across the 100 iterations. The method column indicates the way we estimate the prediction interval. “P” stands for prediction interval $\hat{C}_{I}$ in ((ref)) and “CP” stands for conformal prediction interval $\tilde{C}_{I}$. If “M” is appended, it indicates the prediction intervals are estimated with imputed income for individuals who reported their income in an interval. It can be seen from the table that we would in general, predict a higher income for individuals with more years of education.

The lower bounds of the prediction intervals increase with years of education, as well as the high-density region of the income distribution, which is shown in the columns with $\alpha=0.9$ and a corresponding coverage of 10%. This is consistent with what has been documented in the literature. Having more years of experience tends to increase the predicted income as well. By varying the choice of $\alpha$, we can gain useful insights into the predicted income distribution.

The impact of including interval-valued income data on the prediction interval is more pronounced at smaller $\alpha$ values, as evidenced by the results for $\alpha=0.1$ in the table. This is because, with a 90% coverage target, the tails of the income distribution become more critical. Given that 20% of the observations are reported in ranges and many fall into the highest category, incorporating interval-valued data shifts the right tail of the income distribution. Consequently, the prediction intervals constructed with interval-valued data are significantly different from the ones constructed with imputed income data.

(ref) shows the coverage of the interval valued incomes in the hold-out dataset for the conformal prediction interval $\tilde C$ that are constructed using both interval-valued and exactly reported income data (shown in the “Censored” columns), as well as for the conformal prediction intervals constructed when the interval-valued incomes are imputed (shown in the “Imputed” columns). The significant difference in prediction intervals constructed with interval-valued data and with imputed data may be because there are only five categories for individuals who report their income in a range. This seems somewhat coarse in terms of informing about the true income distribution. Notice that the conformal coverage using the censored model has coverage close to the nominal rate, while the Imputed Coverage is almost always lower. Note that in Tables (ref) and (ref), when we display different levels of coverage \(\alpha\), this provides an empirical estimate of the \((1-\alpha)\)% prediction sets, which can be informative about the shape of the data distribution. Notice that standard reporting of confidence intervals on parameters is meant to summarise sampling uncertainty, while here not only do these conformalized sets take into account sampling uncertainty, but they also estimate prediction sets of the underlying distributions of interest.

{

table[table omitted — 4,861 chars of source]

} {

table[table omitted — 3,169 chars of source]

}

Conclusion

The object of interest in this paper is a prediction set for an outcome given a set of covariates when this outcome is censored. Censoring is a common issue in economic and other data. We first characterise the oracle prediction set, which is the smallest set that maintains a given miscoverage level under censoring. The characterisation leads to a feasible estimation strategy based on the observed sample. We then provide consistent estimators for this oracle prediction set based on nonparametric estimation of the conditional distribution of the observed lower and upper brackets given the conditioning variables. We allow for the prediction set to consist of multiple intervals and consider both random and fixed censoring. Furthermore, we use recent results from conformal inference to obtain a conformal prediction set that maintains a finite sample miscoverage property using a set-based score function. We find that our procedures perform well with simulated data. We also apply our procedures to data from the US Census. There is increasing interest in prediction set estimation, and several intriguing research directions are related to this paper. One promising avenue is to consider prediction sets under partial identification beyond interval censoring. While conformal inference offers finite-sample guarantees on coverage, it remains a research challenge to explore the more traditional statistical inference on the estimated prediction set. Additionally, the design of conformity scores in the context of partial identification presents an interesting and important topic for further investigation.

Declaration of generative AI and AI-assisted technologies in the manuscript preparation process: In preparation for the manuscript, ChatGPT was used to format the figures and tables in the LaTeX document and to provide spelling and grammar checks. GitHub Copilot and Codex, provided code suggestions and auto-completions for the python codes implementing the simulation and empirical study. Refine.ink is also used for proofreading. After using these tools, the authors reviewed and edited the content as needed and take full responsibility for the content of the published article.