Beta regression, Conformal prediction, Logit-Normal, Non-conformity measure, Residuals.

1 Introduction↩︎

In many practical applications, the outcomes of interest take values on a bounded interval \((a, b)\). Examples include rates or proportions of events [1], [2], as well as loss given default, that is, the proportion lost when borrowers default [3]. A common approach to modelling such outcomes with covariates is to first rescale them to the interval \((0,1)\) and then map them to the real line through a function such as the logit [4]. Once the transformed outcomes take values on the real line, they can be modelled using methods from generalized linear models, generalized additive models, distributional regression, or machine learning techniques such as random forests or neural networks. In contrast, beta regression models the outcome directly as beta-distributed on the interval \((0,1)\), incorporating covariates through link functions applied to the mean parameter \(\mu\) and, optionally, the precision parameter \(\phi\), thereby accounting for the bounded and heteroscedastic nature of the data. This is in line with distributional regression [5], where the aim is to model the parameters of a baseline distribution using linear or additive predictors. Since its introduction in the early 2000s [1], [6], beta regression has received considerable attention, leading to efficient software implementations [7], [8], extensions to machine learning contexts [9], [10], and inflated models that accommodate the presence of zeros and/or ones [3], [11]. Other approaches for modelling continuous bounded outcomes include Kumaraswamy regression [12], [13] and continuous binomial regression models [14].

Beyond analyzing the structural and inferential properties of transformation or beta regression models, a question of practical interest is how to construct point and interval predictions for new covariate values [15], a task that, despite its relevance, has received comparatively little attention. To the best of our knowledge, prediction for bounded outcomes has relied primarily on traditional inferential approaches, including bootstrap methods [16]. These methods are typically grounded in asymptotic theory and may require large sample sizes to ensure accurate coverage. Within this prediction framework, conformal prediction [17]–[19] offers a machine learning approach for constructing prediction intervals with finite-sample marginal validity and asymptotic conditional validity. Its core component is a non-conformity measure, a function that quantifies how unusual, or non-conforming, a data point is relative to a model or algorithm trained on the remaining data. The main assumption underlying conformal prediction is exchangeability of the observations, which is often reasonable and allows the method to be coupled with a wide range of statistical models and machine learning algorithms. This assumption can, however, be violated in many regression settings. Extensions of conformal prediction to general regression models have been studied extensively [20], together with several proposals for constructing non-conformity measures in this context [21]. A particular challenge in regression is heteroscedasticity, that is, when the variance (or dispersion) of the response changes with the covariates. In beta regression, the mean and variance are intrinsically linked, so that modelling the mean inherently induces heteroscedasticity. Addressing heteroscedasticity in conformal prediction is itself challenging, since its nature varies across models and outcome types. Several adaptations have been proposed to account for it, including conformalized quantile regression [22], which combines conformal prediction with quantile regression, as well as normalized conformal prediction [18], [21], [23] and Mondrian conformal prediction [24], which both account for heteroscedasticity through the non-conformity measure. A model-based extension of conformal prediction, tailored to distributional regression, was proposed by [25]. Their approach uses the probability integral transform to define non-conformity measures based on an estimated conditional distribution of the response, leading to valid conformal prediction intervals in heteroscedastic settings under some regularity conditions.

Despite the generality of existing conformal prediction frameworks, bounded outcomes present specific challenges that are not adequately addressed by standard theory. First, non-conformity scores based on raw or normalized residuals, standard in heteroscedastic regression [21], [23], produce symmetric prediction intervals that can extend outside the natural support \((0,1)\); truncating to \((0,1)\) restores validity trivially but produces uninformative intervals near the boundaries. Second, heteroscedasticity in bounded outcome regression is not merely an additional modelling assumption but a structural consequence of the support: since \(\text{Var}(Y) \leq \mu(1-\mu)\) for any \(Y \in (0,1)\), where \(\mu = {\mathbb{E}}(Y)\), the variance is constrained to vanish near the boundaries, inducing systematic skewness in the score distribution that generic normalization strategies do not account for. Although the distributional conformal prediction framework of [25] is stated in sufficient generality to accommodate bounded outcomes in principle, it remains generic and does not instantiate regression models tailored to outcomes on \((0,1)\), which requires specialization to a model class with support-constrained heteroscedasticity and verification of the corresponding regularity conditions.

In this work, we propose the use of conformal prediction for regression models with continuous bounded outcomes within the framework of transformation models. This framework includes models of practical interest such as beta regression models and logit-normal regression, which are obtained for different types of transformations of the bounded outcomes. To address heteroscedasticity in both classes of models, we construct non-conformity scores based on model-derived residuals [26]. For comparison, we also consider raw and Pearson residuals, as these are standard choices in practice. This formulation allows for bridging normalized conformal prediction and distributional conformal prediction. Our main contribution is the identification of a quantile-residual-based score that is particularly well suited to bounded outcomes. This score simultaneously accommodates boundedness, heteroscedasticity, and the asymmetry induced near the boundaries of the response space. The novelty of our approach therefore lies not in the residual definitions themselves, but in their principled use for constructing conformal prediction intervals, together with the corresponding validity guarantees. Our development proceeds at three levels of generality. First, we present the split and full conformal prediction algorithms in full generality, so that they apply to any suitable model coupled with a suitable non-conformity score. Second, the theory narrows this scope: we characterize the conditions on the underlying model and score under which the resulting intervals attain finite-sample marginal validity and asymptotic conditional validity. Third, we show that beta regression and logit-normal regression models, coupled with the residual-based scores we propose, satisfy these conditions, so that the general guarantees apply concretely to both classes.

The remainder of the paper is organized as follows. Section 2 describes the modelling approach based within the framework of transformation parametric regression models, and details the main models of interest in this work: beta regression and logit-normal regression. Section 3 introduces the non-conformity measures tailored to these models and presents the split and full conformal prediction algorithms in their general form. Section 4 establishes the conditions for marginal validity and asymptotic conditional validity, and verifies them for the two model classes. Section 5 reports an extensive simulation study of the coverage properties of split and full conformal prediction. Section 6 presents a real-data application and compares conformal prediction with bootstrap-based prediction. Finally, Section 7 concludes with a general discussion and possible extensions. R code and data are available in the GitHub repository https://github.com/ZWU-001/CPBounded.

2 Regression models for continuous bounded outcomes↩︎

In following sections, we focus, without loss of generality, on observations indexed by \(i=1,\dots,n\). For each observation, let \(W_i\in (0,1)\) denote the response random variables, potentially after a linear transformation from a bounded interval \((a,b)\), where \(-\infty < a < b < \infty\). We write \(w_i\) for a realization of \(W_i\). Let \({\boldsymbol{X}}_i \in {\mathbb{R}}^p\) denote the corresponding covariate vector, and let \({\boldsymbol{x}}_i\) be its observed value. To describe the different modelling approaches within a common framework, let \(T:(0,1)\to\mathcal{Y}\subseteq\mathbb{R}\) be a transformation mapping the bounded response space to a subset of the real line. We define the transformed-scale response by \(Y_i=T(W_i)\) and \(y_i=T(w_i)\). There are two cases of practical interest:

• the case where \(T(w) = w\), the identity function, and a parametric model is assumed for such bounded outcomes directly, and

• the case where \(T:(0,1)\to {\mathbb{R}}\) is a continuous map, which allows for using any parametric regression model for unbounded outcomes.

We consider the use of parametric regression models \(Y_i \mid {\boldsymbol{X}}_i = {\boldsymbol{x}}_i \sim F(\cdot ; \boldsymbol{\zeta}, {\boldsymbol{x}}_i)\), where \(F\) is a cumulative distribution function (CDF) supported on the corresponding working scale, and \(\boldsymbol{\zeta}\) denotes the model parameter vector. The conformal prediction algorithms for the new response \(Y_{\text{new}}\) given a new covariate \({\boldsymbol{X}}_{\text{new}}\) in Section 3 will be presented for general choices of \(F\), while Section 4 will characterize the conditions of such models under which theoretical guarantees are obtained. In case (i), the regression model is specified directly on the original bounded response scale. The main example considered here is beta regression, described in Section 2.1; other bounded-response models, such as Kumaraswamy regression and continuous binomial regression, could also be incorporated within the same framework. In case (ii), the response is firstly mapped to an unbounded working scale, allowing the use of parametric, semiparametric, or nonparametric methods, such as a linear predictor, splines, Gaussian processes, or random forests, to name but a few. In this work, we focus on logit-normal regression models with possible heteroscedasticity, as described in Section 2.2. These two models are presented for clarity and practical relevance, but other transformations and regression models can be incorporated into the same framework.

2.1 Beta regression↩︎

We now describe beta regression, which corresponds to case (i) and is of particular interest in practice for modelling continuous responses on the unit interval. Let \(Y\in (0,1)\) be a random variable with beta distribution and shape parameters \(\alpha>0\), and \(\beta>0\). The probability density function (PDF) of \(Y\) is

\[f_{\text{Beta}}(y; \alpha, \beta) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)\Gamma(\beta)} y^{\alpha - 1} (1 - y)^{\beta - 1}, \quad 0 < y < 1.\]

The parameters \((\alpha, \beta)\) can be re-parametrized into \((\mu,\phi)\), where \(\mu = \frac{\alpha}{\alpha+\beta}\in (0,1)\) is the mean and \(\phi = \alpha+\beta > 0\) represents the precision or dispersion parameter [7]. The corresponding pdf is [7]

\[f_{\text{Beta}}(y; \mu, \phi) = \frac{\Gamma(\phi)}{\Gamma(\mu\phi)\Gamma((1 - \mu)\phi)} y^{\mu\phi - 1} (1 - y)^{(1 - \mu)\phi - 1}. \label{eq:betapdf}\tag{1}\] Under this parametrization, \(E(Y) = \mu\), and \(\mathop{\mathrm{Var}}(Y)=\frac{\mu(1-\mu)}{1+\phi}\). We denote the corresponding beta CDF by \(F_{\text{Beta}}(y;\mu,\phi)\). Beta regression (to the mean) is formulated by modelling the mean of each individual, appropriately transformed into \((0,1)\), using a linear predictor [1], [6]. Specifically, the model takes the form \(\varphi(\mu_i) = \theta_0 + {\boldsymbol{x}}_i^{\top}\boldsymbol{\theta}\), where \(\theta_0\) is the intercept, \(\boldsymbol{\theta}\in {\mathbb{R}}^p\) are the regression coefficients, and \(\varphi:(0,1)\to {\mathbb{R}}\) is a continuous monotonic link function (typically, the \(\mathop{\mathrm{logit}}\) function). A natural consequence of modelling the mean in beta regression is that it implicitly defines a model for the variance as well. As a result, beta regression is inherently heteroscedastic. The maximum likelihood estimates (MLEs) of \({\theta_0}\), \({\boldsymbol{\theta}}\) and \(\phi\) can be obtained by maximizing the corresponding likelihood function. The MLEs are not available in closed-form, however, efficient numerical methods and implementations are available in R package betareg and other programming languages [7]. Once the MLEs \(\widehat{\theta_0}\), \(\widehat{\boldsymbol{\theta}}\) and \(\widehat{\phi}\) have been obtained, predictions for a new covariate value \({\boldsymbol{x}}_{\text{new}}\), can be obtained as the predicted mean response \(\widehat{y}_{\text{new}} = \widehat{\mu}_{\text{new}} = \mathop{\mathrm{expit}}(\widehat{\theta_0} + {\boldsymbol{x}}_{\text{new}}^{\top}\widehat{\boldsymbol{\theta}})\), where \(\mathop{\mathrm{expit}}(\cdot)\) denotes the inverse of the \(\mathop{\mathrm{logit}}\) function. A natural extension of the beta regression model to the mean consists of modelling \(\phi\) using covariates. Typically, covariates are included in this parameter through the log-linear model \(\phi_i = \exp\left\{ \eta_0 + {\boldsymbol{x}}_i^{\top}\boldsymbol{\eta}\right\}\), where \(\eta_0\) is the intercept and \(\boldsymbol{\eta}\in {\mathbb{R}}^p\) are the regression coefficients.

2.2 Logit-Normal regression↩︎

We now consider case (ii), focusing on logit-normal regression as a particular example of a transformation regression model. Transformation regression models are obtained by mapping the values \(w_i \in (0,1)\), onto the real line through a continuous, strictly increasing, bijection \(T:(0,1)\to{\mathbb{R}}\). We consider models of the form \[\label{eq:tmodel} y_i=T(w_i) = \psi({\boldsymbol{x}}_{i}) + \epsilon_{i},\tag{2}\] where \(\mu_i=\psi({\boldsymbol{x}}_i) = E[T(w_i)\mid {\boldsymbol{x}}_i]\) is the regression function on the transformed scale and the errors satisfy \(E[\epsilon_i \mid {\boldsymbol{x}}_i] = 0\). In this work, we focus on the parametric specification \(\psi({\boldsymbol{x}}_{i}) = \theta_0 + {\boldsymbol{x}}_{i}^{\top}\boldsymbol{\theta}\) and \(\epsilon_{i} \stackrel{ind.}{\sim} {N}(0, \sigma_i^2)\), where \(\theta_0\) is the intercept and \(\boldsymbol{\theta}\in {\mathbb{R}}^p\) are the regression coefficients in the model for the conditional mean. To allow for heteroscedasticity, the conditional standard deviation is modelled through a separate log-linear predictor, \(\sigma_i = \exp\{\eta_0 + {\boldsymbol{x}}_{i}^{\top}\boldsymbol{\eta}\}\), where \(\eta_0\) and \(\boldsymbol{\eta}\) are corresponding intercept and regression coefficients. The homoscedastic case is recovered by setting \(\boldsymbol{\eta}=\boldsymbol{0}\), so that \(\sigma_1 = \dots = \sigma_n = \sigma\). An implementation of this model can be found in the R package gamlss [27].

Of particular interest is the case where \(T(w) = \mathop{\mathrm{logit}}(w) = \log\left(\dfrac{w}{1-w}\right)\), which leads to classical logit-normal regression. That is, this formulation corresponds to a linear regression model on the logit scale. Consequently, this strategy can be coupled with any tool available for classical linear regression. After fitting this model, the predicted values for a new covariate value \({\boldsymbol{x}}_{\text{new}}\) can be obtained as \(\widehat{y}_{\text{new}} = \widehat{\mu}_{\text{new}}=\widehat{\theta}_0+{\boldsymbol{x}}_{\text{new}}^{\top}\widehat{\boldsymbol{\theta}}\) and accordingly, \(\widehat{w}_{\text{new}} = \mathop{\mathrm{expit}}(\widehat{y}_{\text{new}})=\mathop{\mathrm{expit}}(\widehat{\theta}_0+{\boldsymbol{x}}_{\text{new}}^{\top}\widehat{\boldsymbol{\theta}})\), where \(\widehat{\theta}_0\) and \(\widehat{\boldsymbol{\theta}}\) denote the estimated intercept and regression coefficient vector.

3 Conformal prediction for transformation and regression models↩︎

In this section, we present the proposed strategy for constructing conformal prediction intervals for continuous bounded responses using the regression models introduced in Section 2. Throughout this section, let \(\mathcal{D}_n=\{(w_i,{\boldsymbol{x}}_i)\}_{i=1}^n\) denote the observed data, where \(w_i\in(0,1)\), and let \(y_i=T(w_i)\) be the corresponding response on the working scale. We denote \(W_{\text{new}}\in(0,1)\) as the new bounded response associated with the covariate vector \({\boldsymbol{X}}_{\text{new}}={\boldsymbol{x}}_{\text{new}}\), and let \(Y_{\text{new}}=T(W_{\text{new}})\) and \(y_{\text{new}}\) be the corresponding response on the working scale and its realization. Let \(F(\cdot; \widehat{\boldsymbol{\zeta}},{\boldsymbol{x}}_i)\) denote the fitted CDF of the working-scale response \(Y_i\mid{\boldsymbol{X}}_i={\boldsymbol{x}}_i\), where \(\widehat{\boldsymbol{\zeta}}\) denotes the estimated parameters based on the \(\mathcal{D}_n\). For both regression models, as introduced in Section 2, \(\widehat{\boldsymbol{\zeta}}=(\widehat{\theta}_0, \widehat{\boldsymbol{\theta}},\widehat{\eta}_0,\widehat{\boldsymbol{\eta}})\) considering heteroscedasticity. When \(T\) is the identity map, so that \(y_i=w_i\), this is the fitted CDF of the original bounded response, as in beta regression. When \(T\) maps \((0,1)\) to the real line \(\mathbb{R}\), this is the fitted CDF of the transformed response, as in logit-normal regression. The conformal construction is carried out on the working scale; the resulting interval is then mapped back to the original bounded response scale whenever \(T\) is not the identity. For notational convenience, we use the working-scale pairs \(\{(y_i,{\boldsymbol{x}}_i)\}_{i=1}^n\) in the remainder of this section. We begin by introducing residual-based non-conformity measures, which are derived from model-specific residuals in general regression settings [26]. These measures are tailored to each model and account for heteroscedasticity where appropriate. We then describe the general split conformal and full conformal algorithms and focus on the models and the non-conformity measures presented in Section 3.1. Without loss of generality, we set \(\widehat{y}_1, \dots, \widehat{y}_n\) to be the point prediction for working-scale observations \(1, \dots, n\) generated by the fitted regression model. Specifically, for the beta regression model, we have \(\widehat{y}_i=\mathop{\mathrm{expit}}(\widehat{\theta}_0+{\boldsymbol{x}}_{i}^{\top}\widehat{\boldsymbol{\theta}})\), where \(\widehat{\theta}_0\) and \(\widehat{\boldsymbol{\theta}}\) denote the estimated intercept and regression coefficient vector, respectively, obtained from the fitted regression model. For the logit-normal regression model, we have \(\widehat{y}_i=\widehat{\theta}_0+{\boldsymbol{x}}_{i}^{\top}\widehat{\boldsymbol{\theta}}\), and accordingly, the prediction in the original scale would be \(\widehat{w}_i=\mathop{\mathrm{expit}}(\widehat{\theta}_0+{\boldsymbol{x}}_{i}^{\top}\widehat{\boldsymbol{\theta}})\).

3.1 Non-conformity measures↩︎

The three scores introduced below are presented in increased order of complexity in terms of how they incorporate information about the fitted model’s structure. Raw residuals 3 target only the conditional location; Pearson residuals 4 additionally standardize by the conditional scale and thus adapt to heteroscedasticity; and quantile residuals 5 use the entire fitted conditional distribution through the probability integral transform, thus also accommodating asymmetry and tail behaviour.

• Raw residuals: A raw-residuals based non-conformity measure is defined as: \[\label{eq:raw95res} r_i = \vert y_i - \widehat{y}_i\vert.\tag{3}\]

• Pearson residuals: A Pearson-residuals based non-conformity measure is defined as: \[\label{eq:pearson95res} r^{P}_{i} = \frac{|y_i - \widehat{y}_i|}{\sqrt{\widehat{\text{Var}}(Y_i \mid {\boldsymbol{x}}_i)}},\tag{4}\] where for beta regression, \(\widehat{\text{Var}}(Y_i \mid {\boldsymbol{x}}_i) = \frac{\widehat{\mu}_i(1-\widehat{\mu}_i)}{1+\widehat{\phi}_i}\); \(\widehat{\mu}_i=\widehat y_i\) is the fitted mean and \(\widehat{\phi}_i\) is the fitted precision parameter; for the logit-normal regression, \(\widehat{\operatorname{Var}}(Y_i\mid {\boldsymbol{x}}_i)=\widehat{\sigma}_i^2\), where \(\widehat{\sigma}_i\) is the fitted conditional standard deviation. This type of residual adjusts the difference between the observed value and the fitted mean prediction by the estimated standard deviation, assuming it exists, thus accounting for heteroscedasticity through localized weighting. The resulting non-conformity measure resembles those used in normalized conformal prediction, where absolute raw residuals are scaled by general functions \(\tau_i({\boldsymbol{x}}_i)\). Various choices for \(\tau_i\) have been explored in the context of regression models [21], [24], with the goal of capturing different forms of heteroscedasticity.

• Quantile residuals: A quantile-residuals based non-conformity measure is defined as: \[\label{eq:res95quantile} r_i^{Q} = \left \vert \Phi^{-1} \{ F(y_i; \widehat{\boldsymbol{\zeta}}, {\boldsymbol{x}}_i) \} \right \vert.\tag{5}\] The choice of \(\Phi^{-1}\) is convenient because it makes the quantile residual reduce to the standardized (Pearson) residual in the special case of a heteroscedastic normal model: there, the probability integral transform gives \(F(y_i;\widehat{\boldsymbol{\zeta}}_i,{\boldsymbol{x}}_i)=\Phi\{(y_i-\widehat{y}_i)/\widehat{\sigma}_i\}\), so that Equations 4 and 5 coincide, \(r^{P}_{i}=r^{Q}_{i}=|y_i-\widehat{y}_i|/\widehat{\sigma}_i\). This identity is specific to the Gaussian case; for non-Gaussian models such as beta regression the two scores differ, with the quantile residual additionally capturing the asymmetry and boundary behaviour of the conditional distribution.

These scores form a coherent family rather than an generic list. This is illustrated by the heteroscedastic normal case, in which the Pearson and quantile scores coincide: the quantile residual generalizes the Pearson residual to non-Gaussian, bounded conditional distributions, and recovers the distributional conformal prediction score of [25] as an interpretable special case (using the uniform distribution instead of Gaussian). In this sense, the non-conformity measure is not chosen generically but aligned with the model used to fit the data, a residual-based framework that links the conformal score to the model’s distributional structure and bridges normalized and distributional conformal prediction.

The interplay between heteroscedasticity and the choice of non-conformity score is central here. In transformation regression, heteroscedasticity is introduced explicitly through a covariate-dependent scale \(\sigma_i\) (Section 2.2), whereas in beta regression it is intrinsic, since modelling the mean \(\mu_i\) already shapes the variance, with further flexibility when \(\phi_i\) is covariate-dependent (Section 2.1). For bounded outcomes, however, the difficulty is not only one of changing variance: near the boundaries of \((0,1)\) the conditional distribution becomes asymmetric and heavier-tailed, so the symmetric Pearson-type score 4 may reflect boundary-induced skewness rather than local scale. Quantile residuals 5 instead adapt to the full conditional distribution, accounting for scale, asymmetry, and tail behaviour in a single principled measure, which is why we recommend them for (transformed) bounded outcomes.

Our framework is organized by a single principle: the non-conformity score is aligned with the residual structure of the fitted model. The quantile-residual score 5 is the member of this family obtained from the probability integral transform, which is also why it relates to the baseline distributional conformal score of [25]. This coincidence, however, is a property of that particular residual and not of the framework: other residuals for beta regression [26], [28], [29] lead to non-conformity measures that are valid within our construction yet do not reduce to a rank-based distributional conformal score. The two viewpoints therefore intersect at the quantile residual rather than one containing the other. Among these residuals, the members differ in how directly they accommodate heteroscedasticity; for bounded outcomes we recommend the quantile residual because it adapts to the full conditional distribution, capturing scale, asymmetry, and boundary behaviour in a single measure, whereas several alternatives capture only some of these. Normalized non-conformity measures developed for heteroscedastic regression [21] can also be embedded in the framework, though they are not generally tailored to the support constraint of bounded outcomes.

We note that these quantile residuals [26] are distinct from the quantile scores used in conformalized quantile regression [22]. The latter are built from conditional quantile functions of the response, typically estimated by quantile regression, and measure the signed distance of the response from an estimated quantile band. By contrast, the score in 5 is obtained from the fitted conditional distribution \(F(\cdot; \widehat{\boldsymbol{\zeta}}, {\boldsymbol{x}}_i)\) through the probability integral transform, so that conditional quantiles enter only at interval construction, when \(F\) is inverted; in this sense the score is closer in spirit to the distributional conformal prediction framework [25] rather than the quantile-regression family. Both approaches address heteroscedasticity, but through different mechanisms: covariate-dependent quantile bands in the former, and the covariate-dependent conditional distribution in the latter, which for bounded outcomes additionally captures boundary-induced asymmetry.

3.2 Split and full conformal prediction algorithms↩︎

We now present split and full conformal prediction algorithms for bounded continuous outcomes. The algorithms are stated in relatively general terms. These algorithms will later be used within the two modelling frameworks described in Sections 2.1-2.2 (beta regression and logit-normal regression), together with the residual-based scores in Section 3.1, with particular emphasis on the quantile-residual score. The corresponding model-specific interval constructions are detailed in Sections 3.3 and 3.4.

The algorithms are stated in general terms through a regression model, the non-conformity score \(S\), and the fixed working-scale mapping \(T\) used before model fitting. In our framework, as is discussed earlier, for beta regression no mapping is applied, so \(y=T(w)=w\), \(\mathcal{Y}=(0,1)\), and the model is fitted directly to the bounded response. For logit-normal regression, \(y=T(w)=\mathop{\mathrm{logit}}(w)\), \(\mathcal{Y}=\mathbb{R}\), \(T^{-1}\) is the \(\mathop{\mathrm{expit}}\) map, and the model is fitted to the transformed responses. In both cases, the conformal construction is carried out on the scale used for model fitting, while the final prediction interval is reported on the original bounded response scale. For clarity, we write \(S(y_i,{\boldsymbol{x}}_i;\widehat{\mathcal{M}}(\widehat{\boldsymbol{\zeta}}))\) in Algorithm 1 for the non-conformity score assigned to a candidate pair \((y_i,{\boldsymbol{x}}_i)\), where \(y_i \in\mathcal{Y}\) is a response value on the working scale, \({\boldsymbol{x}}_i \in\mathbb{R}^p\) is the associated covariate vector, and \(\widehat{\mathcal{M}}(\widehat{\boldsymbol{\zeta}})\) denotes the fitted regression model with estimated parameter vector \(\widehat{\boldsymbol{\zeta}}\). The fitted model provides the quantities needed to evaluate the score, such as the fitted conditional location, conditional variance, or fitted conditional CDF. \(\widehat{\mathcal{M}}^{(y)}(\widehat{\boldsymbol{\zeta}})\) in Algorithm 2 represents the fitted model with the candidate value \(y\).

The split conformal method is computationally efficient, as it fits the regression model only once, making it well suited for large datasets. However, it relies on a random split of the data into training and calibration sets, which may reduce statistical efficiency and introduce variability in the resulting prediction intervals. In contrast, the full conformal method uses the entire dataset for both model fitting and calibration by refitting the model for each candidate response value. This typically produces narrower intervals but incurs substantially higher computational cost due to repeated model training across the prediction grid. These can be observed from empirical results in Section 5–6.

3.3 Split conformal prediction intervals construction↩︎

The Algorithm 1 is generic as they can be used with various regression models and score functions. Depending on different models and score types, the resulting split conformal prediction intervals may take different forms. In all scenarios, the prediction interval is constructed by solving the inequality

\[\left\{ y: S(y,{\boldsymbol{x}}_{\text{new}};\widehat{\mathcal{M}}(\widehat{\boldsymbol{\zeta}}))\le q_{1-\alpha} \right\}. \label{eq:scores}\tag{6}\]

With the new covariate \({\boldsymbol{x}}_{\text{new}}\) and the fitted model \(\widehat{\mathcal{M}}(\widehat{\boldsymbol{\zeta}})\), we can obtain the new point prediction \(\widehat{y}_{\text{new}}\). The explicit form of this inequality depends on the chosen non-conformity score. For example, the raw residual score depends only on the fitted conditional location \(\widehat y_{\text{new}}\), the Pearson residual score additionally depends on the fitted conditional variance or scale at \({\boldsymbol{x}}_{\text{new}}\), and the quantile residual score depends on the full fitted conditional CDF \(F(\cdot;\widehat{\boldsymbol{\zeta}},{\boldsymbol{x}}_{\text{new}})\). We now specify the resulting interval forms for the two modelling frameworks considered in this paper.

• Beta regression model: For this model, no transformation is applied to the response. Two variants of the Beta regression model are considered: one in which only the mean parameter \(\mu\) depends on covariates, and another one where both the mean \(\mu\) and the precision parameter \(\phi\) are modelled as functions of covariates. For each scenario, two types of nonconformity scores are employed: the Pearson residual in 4 and the quantile residual in 5 . Under the Pearson residual, 6 becomes

  \[\left\lbrace y: \frac{\left \vert y-\widehat{y}_{\text{new}} \right \vert}{\sqrt{\widehat{\text{Var}}(Y_{\text{new}} \mid {\boldsymbol{x}}_{\text{new}})}} \leq q_{1-\alpha} \right\rbrace.\]

  This results in the interval

  \[\left(\max\left(0,\widehat{y}_{\text{new}}-q_{1-\alpha}\sqrt{\widehat{\text{Var}}(Y_\text{new} \mid {\boldsymbol{x}}_\text{new})}),\min(1, \widehat{y}_{\text{new}}+q_{1-\alpha}\sqrt{\widehat{\text{Var}}(Y_\text{new} \mid {\boldsymbol{x}}_\text{new})}\right)\right),\]

  where \(\sqrt{\widehat{\text{Var}}(Y_\text{new} \mid {\boldsymbol{x}}_\text{new})}=\sqrt{\widehat{y}_{\text{new}}(1 - \widehat{y}_{\text{new}})/(1 + \widehat{\phi}_{\text{new}})}\) and \(\widehat{\phi}_{\text{new}}\) denotes the estimated dispersion of the new observation (including both homoscedastic and heteroscedastic scenarios). In this setting, the resulting interval may extend beyond the unit interval \((0, 1)\). This limitation arises because the interval construction accounts only for heteroscedasticity through the variance term \(\sqrt{\widehat{\text{Var}}(Y_\text{new}\mid {\boldsymbol{x}}_\text{new})}\), while ignoring the skewness of the beta distribution near the boundaries. When \(\widehat{y}_{\text{new}} \to 0^+\), residuals exhibit positive skewness; when \(\widehat{y}_{\text{new}} \to 1^-\), negative skewness prevails. To ensure that the prediction interval respects the bounded response space, we truncate it to \((0,1)\). This truncation is not part of the theoretical construction in Section 4 and does not affect the validity guarantees; since \(Y_{\text{new}} \in (0,1)\) almost surely for beta models, the coverage event is unchanged. Under the quantile residual, 6 becomes

  \[\left\{y: \left \vert \Phi^{-1} \{ F_{\text{Beta}}(y; \widehat{\mu}_{\text{new}}, \widehat{\phi}_{\text{new}}) \} \right \vert \leq q_{1-\alpha}\right\}.\]

  This results in the interval:

  \[(F^{-1}_{\text{Beta}}(\Phi(-q_{1-\alpha}); \widehat{\mu}_{\text{new}},\widehat{\phi}_{\text{new}}), F^{-1}_{\text{Beta}}(\Phi(q_{1-\alpha}); \widehat{\mu}_{\text{new}},\widehat{\phi}_{\text{new}})),\]

  where \(F^{-1}_{\text{Beta}}\) is the inverse beta CDF. This inverse transformation ensures the interval remains within \((0, 1)\), naturally accounts for skewness, and eliminates the need for truncation.

• Logit-normal regression model: Under the homoscedastic setting (\(\sigma_1 = \dots = \sigma_n = \sigma\)), we use the absolute raw residual from 3 as the score function. Here \(\widehat{y}_{\text{new}}\) denotes the fitted conditional mean on the transformed scale. 6 then becomes

  \[\{y: \left \vert y-\widehat{y}_{\text{new}} \right \vert \leq q_{1-\alpha}\}.\]

  This results in the prediction interval:

  \[(\mathop{\mathrm{expit}}(\widehat{y}_{\text{new}}-q_{1-\alpha}),\mathop{\mathrm{expit}}(\widehat{y}_{\text{new}}+q_{1-\alpha})).\]

  The function \(\mathop{\mathrm{expit}}\) is used to map the prediction interval in \(\mathbb{R}\) back to the \((0,1)\) scale. Under the heteroscedastic setting (\(\sigma_i\)’s are covariate-dependent), we use the absolute Pearson residual from 4 . In this case, inequality 6 becomes

  \[\left\{y: \frac{\left \vert y-\widehat{y}_{\text{new}} \right \vert}{\widehat{\sigma}_{\text{new}}} \leq q_{1-\alpha}\right\}.\]

  This results in the interval form:

  \[\left(\mathop{\mathrm{expit}}(\widehat{y}_{\text{new}}-q_{1-\alpha}\widehat{\sigma}_{\text{new}}), \mathop{\mathrm{expit}}(\widehat{y}_{\text{new}}+q_{1-\alpha}\widehat{\sigma}_{\text{new}})\right),\]

  where \(\widehat{\sigma}_{\text{new}}\) is the estimated standard deviation of the new observation.

3.4 Full conformal prediction intervals construction↩︎

A straightforward implementation of full CP (Algorithm 2) evaluates the non-conformity score at every candidate response value \(y\) on a predefined grid \(\mathcal{Y}\), augmenting and refitting the model for each candidate, and then identifies the smallest and largest values satisfying the conformal inclusion criterion. Although conceptually simple, this exhaustive search can be computationally demanding when the required precision level is high or the fitted model is complex. We therefore employ an adaptive interval-finding algorithm that locates the conformal interval boundaries without evaluating every grid point. In the case of the beta regression models, we numerically solve for the interval within the unit interval \((0,1)\) by applying an interval-finding algorithm with a user-specified tolerance \(\varepsilon>0\). Based on Algorithm 2, our aim is to identify the set:

\[\mathcal{C}^{\text{full}}_{1-\alpha}({\boldsymbol{x}}_{\text{new}}) = \{ y \in (0,1) : S(y,{\boldsymbol{x}}_{\text{new}};\widehat{\mathcal{M}}^{(y)}(\widehat{\boldsymbol{\zeta}})) \leq q_{1-\alpha} \}.\]

To check whether a candidate value belongs to \(\mathcal{C}^{\text{full}}_{1-\alpha}\), we define an indicator-type function:

\[I: (0, 1) \rightarrow \{\texttt{TRUE}, \texttt{FALSE}\},\]

where \(I(\xi) = \texttt{TRUE}\) if the augmented nonconformity score for candidate \(y\) satisfies the inclusion criterion, and FALSE otherwise. Our objective is to find the values \(L\) and \(U\), i.e., lower bound and upper bound for the prediction interval, such that \(I(\xi) = \texttt{TRUE}\) for all \(\xi \in [L,U] \subset (0,1)\). The core idea here is to locate a point inside the prediction interval, i.e., satisfying the inclusion criterion. Our first step is to decide an initial guess \(\xi_0\in(0,1)\). The algorithm distinguishes two cases depending on whether \(\xi_0\) lies inside the interval.

3.4.0.1 Case 1: The initial guess lies within the inclusion region (\(I(\xi_0) = \texttt{TRUE}\))

In this case, a valid point inside the conformal prediction set is known. The algorithm proceeds as follows:

• To find the lower endpoint \(L\), it performs a binary search on the interval \((0, \xi_0)\), repeatedly halving subintervals to locate the point where \(I(\xi)\) transitions from FALSE to TRUE within tolerance \(\varepsilon\).

• To find the upper endpoint \(U\), it performs a similar binary search on \((\xi_0, 1)\), locating the transition from TRUE to FALSE.

3.4.0.2 Case 2: The initial guess lies outside the inclusion region (\(I(\xi_0) = \texttt{FALSE}\))

In this case, we need to perform an exploratory search to find a point inside the interval.

• Initialization: A coarse grid is set up over \((0,1)\), typically starting with points \(0\), \(0.5\), and \(1\). The endpoints \(0\) and \(1\) are assigned with FALSE values.

• Grid Expansion: If no TRUE values are detected on the initial grid, i.e., \(I(0.5)=\texttt{FALSE}\), the grid is refined iteratively by inserting midpoints between all adjacent pairs. This process continues until a point \(\xi' \in (0,1)\) is found such that \(I(\xi') = \texttt{TRUE}\).

• Boundary Localization: Once a valid inclusion point \(\xi'\) is located, its immediate FALSE neighbors on both sides are identified. These provide bounds within which the transitions occur.

• Binary Search: The lower and upper boundaries of the inclusion set are then refined within tolerance \(\varepsilon\) using the binary search method from Case 1.

For logit-normal regression models under the full conformal setting, the candidate response values need to be searched over a grid on the transformed real line. One choice to define the candidate grid \(\mathcal{Y}\) is using an expanded interval around, e.g., the classical prediction bounds in normal linear regression:

\[\left[l-\rho_{\text{grid}}(u-l), u+\rho_{\text{grid}}(u-l)\right],\]

where \(l\) and \(u\) denote the lower and upper bounds of the classical prediction interval on the transformed scale. The parameter \(\rho_{\text{grid}}\) is only a numerical grid-expansion parameter for this specific full conformal implementation. It is not a conformal tuning parameter and does not enter the theoretical validity result. Any sufficiently large choices would lead to the same conformal interval up to the grid resolution, provided that the grid contains the conformal inclusion region. In practice, this can be checked by verifying that the conformal inclusion rule is false at both endpoints; otherwise, the grid needs to be enlarged. To construct the prediction interval, similarly, we firstly set an initial guess \(\xi_0\in\mathbb{R}\). If \(I(\xi_0) = \texttt{TRUE}\), we perform the binary search on the interval \((l-\rho_{\text{grid}}(u-l), \xi_0)\) and the interval \((\xi_0,u+\rho_{\text{grid}}(u-l))\) similar to Case 1. If \(I(\xi_0) = \texttt{FALSE}\), we adopt a similar strategy outlined in Case 2. The only difference is that instead of inserting midpoints of candidate values, we insert the values corresponding to the floor of the midpoint indices.

To reduce computational cost, evaluations of \(I(\xi)\) are stored and reused (i.e., memorized), thereby avoiding redundant computation during grid refinement and binary search phases. The theoretical results in Section 4 concern the exact conformal set \(\mathcal{C}^{\mathrm{full}}_{1-\alpha}\). The procedure described above is a numerical implementation that approximates the endpoints of this set up to the user-specified tolerance \(\varepsilon\), provided that the initial search domain contains the conformal inclusion region. Accordingly, neither \(\varepsilon\) nor \(\rho_{\mathrm{grid}}\) enters the theoretical validity results.

4 Theoretical Guarantees↩︎

We now establish theoretical guarantees for conformal prediction intervals for bounded outcomes based on Algorithms 1–2. Proposition [prop:marginal] establishes finite-sample marginal validity for both split and full conformal prediction under exchangeability. Although exchangeability-based marginal validity for generic conformal prediction is well established [19], these guarantees are not immediate in our setting. In particular, they require verifying that the non-conformity scores (S) are symmetric and almost surely free of ties under the bounded-outcome models considered. Since this property does not automatically hold for arbitrary model-based residuals, it is established in Lemma 1 for the quantile residuals in 5 . Establishing asymptotic conditional validity is even more challenging. The bounded support of the response induces heteroscedasticity and boundary-driven asymmetry in the conditional distribution, so conditional validity requires regularity conditions tailored to the distribution of the fitted scores rather than a direct application of existing results. Proposition [prop:conditional] establishes asymptotic conditional validity under regularity conditions.

We first outline the framework and notation through the following assumptions.

• For \(i=1,\dots,n\), let \(W_i\in(0,1)\) denote the original bounded response and let \(Y_i=T(W_i)\) denote the corresponding working-scale response, where \(T:(0,1)\to\mathcal{Y}\subseteq\mathbb{R}\) is a continuous and strictly increasing mapping. Let \({\boldsymbol{X}}_i\in\mathbb{R}^p\) denote the covariate vector, and define \({\boldsymbol{Z}}_i=(Y_i,{\boldsymbol{X}}_i)\) and the testing point \({\boldsymbol{Z}}_{\text{new}}=(Y_{\text{new}},{\boldsymbol{X}}_{\text{new}})\). The corresponding realizations are denoted by \({\boldsymbol{z}}_i=(y_i,{\boldsymbol{x}}_i)\) and \(y_i=T(w_i)\), where \({\boldsymbol{z}}_i\) and \({\boldsymbol{x}}_i\) denote the corresponding realizations.

• We consider a parametric regression model for the working-scale response, \(Y_i\mid {\boldsymbol{X}}_i={\boldsymbol{x}}_i \sim F(\cdot;\boldsymbol{\zeta},{\boldsymbol{x}}_i)\), where \(F(\cdot;\boldsymbol{\zeta},{\boldsymbol{x}}_i)\) is a continuous CDF supported on \(\mathcal{Y}\), and \(\boldsymbol{\zeta}\) denotes the model parameter vector. Let \(\widehat{\boldsymbol{\zeta}}\) be the estimator of \(\boldsymbol{\zeta}\), computed from the observed sample by a fitting rule that is invariant to permutations of its inputs, such as maximum likelihood estimation. Under possible model misspecification, we denote the probability limit of \(\widehat{\boldsymbol{\zeta}}\) by \(\boldsymbol{\zeta}^\star\), which is assumed to exist and to be an interior point of the parameter space.

• Suppose that the true conditional distribution of \(y_i \mid {\boldsymbol{x}}_i\) is denoted by \(\Psi(y_i \mid {\boldsymbol{x}}_i)\). Under the correct model specification, \(\Psi(y_i \mid {\boldsymbol{x}}_i)=F(y_i; \boldsymbol{\zeta}^*,{\boldsymbol{x}}_i)\), whereas under the misspecified model, \(\Psi(y_i \mid {\boldsymbol{x}}_i) \neq F(y_i; \boldsymbol{\zeta}^*, {\boldsymbol{x}}_i)\).

The above assumptions cover both model frameworks considered in this paper. For beta regression, \(T\) is the identity map. Hence, \(Y_i=W_i\), \(\mathcal{Y}=(0,1)\), and \(F(\cdot;\boldsymbol{\zeta},{\boldsymbol{x}}_i)\) is the conditional CDF of the bounded response itself, which corresponds to \(F(y_i;\boldsymbol{\zeta},{\boldsymbol{x}}_i)=F_{\text{Beta}}(y_i;\mu_i,\phi_i)\), where \(\mu_i\) and possibly \(\phi_i\) are modelled through covariates. For the logit-normal regression model, \(T\) is a strictly increasing map (\(\mathop{\mathrm{logit}}\)) from \((0,1)\) to an unbounded working scale \(\mathbb{R}\). In this case, \(Y_i\) is the transformed response and \(F(\cdot;\boldsymbol{\zeta},{\boldsymbol{x}}_i)\) denotes the conditional CDF of \(Y_i\mid {\boldsymbol{X}}_i={\boldsymbol{x}}_i\) on the transformed scale, which corresponds to the Gaussian CDF \(F(y_i;\boldsymbol{\zeta},{\boldsymbol{x}}_i)=\Phi\left(\frac{y_i-\mu_i}{\sigma_i}\right)\) associated with the fitted regression model on the \(\mathop{\mathrm{logit}}\) scale. These mild assumptions suffice to establish finite-sample marginal validity for the conformal prediction intervals produced by Algorithm 1 (split conformal) and Algorithm 2 (full conformal), holding in the presence of heteroscedasticity.

As shown in Lemma 1 in the Appendix 9, the quantile residual score \(r_i^Q\) in 5 satisfies the no-ties condition almost surely under both the beta and logit-normal regression models, thus representing examples where the marginal validity results in Proposition [prop:marginal] hold. In both frameworks, the exchangeability requirement of Proposition [prop:marginal] holds under the standard sampling assumption underlying the prediction problem. Specifically, if the observed pairs \((W_i,{\boldsymbol{X}}_i)\), together with \((W_{\text{new}},{\boldsymbol{X}}_{\text{new}})\), are drawn exchangeably from a common joint distribution \(P\), as is routinely assumed in both beta and logit-normal regression, then exchangeability is preserved on the working scale.

To extend our results to asymptotic conditional validity, consider now the following conditions.

• The estimated conditional CDF converges uniformly to its probability limit: \[\sup_{Y_i, {\boldsymbol{X}}_i} \left|F(Y_i; \boldsymbol{\zeta}^*,{\boldsymbol{X}}_i) - F(Y_i; \widehat{\boldsymbol{\zeta}}, {\boldsymbol{X}}_i)\right| = o_p(1), \quad \text{for }i = 1,\dots,n, \text{new},\] where “\(\text{new}\)” denotes the index for a new observation.

• Suppose that the fitted regression model \(\widehat{\mathcal{M}}\) converges in probability to a limit model \({\mathcal{M}}^{(*)}\). Let \(V_i := S(Y_i,{\boldsymbol{X}}_i;{\mathcal{M}}^{(*)})\) denote the oracle non-conformity score corresponding to \({\mathcal{M}}^{(*)}\). The conditional score distribution \(G_\star(v) := \Pr(V_{\mathrm{new}} \le v \mid \mathbf{X}_{\mathrm{new}} = \mathbf{x}_{\mathrm{new}})\) admits a density \(g_\star\) that is positive on a neighbourhood \(\mathcal{N}\) of the \((1-\alpha)\)-quantile \(q^{*}_{1-\alpha} = G_\star^{-1}(1-\alpha)\).

• For \(i = 1,\dots,n\), let \(G_i(v) := \Pr(V_i \le v)\) denote the distribution function for \(V_i\), and set \(m_n(v) := \frac{1}{n}\sum_{i=1}^{n} G_i(v)\). There exists a neighbourhood \(\mathcal{N}\) of \(q_{1-\alpha}^*\) such that \[m_n(v) \;\longrightarrow\; G_\star(v), \quad\text{for every } v \in \mathcal{N}.\]

The following result shows that these conditions are sufficient to guarantee the asymptotic conditional validity of split and full conformal prediction under the use of the quantile residual score \(r^Q_i\).

Condition C1 requires the estimated conditional CDF to converge uniformly to its probability limit, in the sense that the maximum discrepancy between \(F(Y_i; \widehat{\boldsymbol{\zeta}}, {\boldsymbol{X}}_i)\) and \(F(Y_i; \boldsymbol{\zeta}^*,{\boldsymbol{X}}_i)\) vanishes in probability. For the beta regression model, sufficient conditions for C1 include a compact parameter space \(\Delta \subset \mathbb{R}^{2(p+1)}\) and the existence of constants \(M > 0\) and \(\varepsilon > 0\) such that \(\|\mathbf{X}_i\| \leq M\) and \(Y_i \in [\varepsilon, 1-\varepsilon]\) for all \(i\). For the logit-normal regression model, sufficient conditions are a compact parameter space \(\Delta \subset \mathbb{R}^{2(p+1)}\) and \(\|\mathbf{X}_i\| \leq M\), since \(Y_i\in\mathbb{R}\) and the Gaussian CDF is uniformly continuous in its parameters. For beta regression, Condition C1 is more general than these sufficient conditions, as it permits outcomes arbitrarily close to the boundaries of \((0,1)\). Condition C2 is a mild regularity requirement ensuring that \(G_{\star}\) is strictly increasing at the relevant quantile; it holds for beta distributions whenever the \((1-\alpha)\)-quantile of the score lies in the interior of the support. Condition C3 requires that the conditional score distribution be covariate-invariant in a neighbourhood of the target quantile. That is, that the oracle score \(V=\left|\Phi^{-1}(F(Y;\boldsymbol{\zeta}^*,{\boldsymbol{x}}))\right|\) be asymptotically pivotal there (fixed distribution), and that the averaged score distribution \(m_n(v)=\frac{1}{n}\sum_{i=1}^n G_i(v)\) converge to this common law \(G_\star\). Although C3 is in fact a strong homogeneity condition, and the reason is structural: the empirical quantile \(\widehat{q}_{1-\alpha}\) is computed from the calibration scores alone and therefore does not depend on \({\boldsymbol{x}}_{new}\), so it admits a single probability limit. A covariate-free limit can coincide with the conditional target \(q_{1-\alpha}^* = G_\star^{-1}(1-\alpha)\) at every test point only if that target is itself covariate-free near the quantile; C3 is precisely this requirement. Under correct specification it holds automatically for both model frameworks: the probability integral transform \(F(Y;\boldsymbol{\zeta}^*; {\boldsymbol{x}})\) is \(\mathrm{Uniform}(0,1)\) at every \({\boldsymbol{X}}={\boldsymbol{x}}\), so \(V\sim\mathrm{HalfNormal}(0,1)\) irrespective of the covariates, \(G_\star\) does not depend on \({\boldsymbol{x}}_{new}\). Under misspecification the pivot breaks, the conditional score law varies with \({\boldsymbol{x}}_{new}\), and C3 becomes a genuine and restrictive assumption, which is a well-known limitation that non-localized conformal prediction attains conditional, as opposed to marginal, validity only under score homogeneity [30].

5 Simulation study↩︎

In this section, we conduct a simulation study to evaluate the empirical marginal coverage (at the 90% nominal level) and the average width of prediction intervals constructed using split and full conformal prediction methods across four modelling frameworks:

1. logit-normal regression,

2. heteroscedastic logit-normal regression,

3. beta regression with covariates affecting only the mean (\(\mu\)) parameter,

4. beta regression with covariates affecting both the mean (\(\mu\)) and precision (\(\phi\)) parameters.

The experiment is designed as it follows.

• Covariate generation: The covariate vectors \({\boldsymbol{x}}_i = (x_{i1},x_{i2},x_{i3})^{\top}\), \(i=1,\dots,n\), are generated from a multivariate normal distribution: \({\boldsymbol{x}}_i \stackrel{i.i.d.}{\sim} N_3 (\boldsymbol{0},\boldsymbol{\Sigma})\), where \(\boldsymbol{\Sigma}\) is a compound symmetric matrix with \(\Sigma_{jk} = 1\) if \(j=k\), and \(\Sigma_{jk} = 0.5\) if \(j \neq k\).

• Sample sizes and replications: Sample sizes considered are \(n \in \{50, 100, 500, 1000\}\). For each simulation setting, the procedure is repeated across multiple Monte Carlo replications to evaluate empirical coverage. Specifically, 10,000 replications are used for \(n=50\) to improve stability at small sample sizes, while 1,000 replications are used for the other values of \(n\). In each replication, a new test point is generated to assess whether the constructed prediction interval contains the true response value and the interval width.

• Response generation: We consider four data-generating mechanisms (S1-S4) corresponding to the four models under evaluation:

  1. logit-normal regression model (S1): To ensure comparability with the beta regression model (S3), the standard deviation \(\sigma\) in the transformation model is calibrated to approximately match the marginal mean and variance of a beta distribution with covariate-dependent mean \(\mu\). This is achieved by minimizing a loss function that balances squared differences in empirical means and variances between simulated logit-normal and beta responses. The resulting calibration yields approximate pairings of \(\sigma = 1.50,0.90,0.63,0.45\) corresponding to \(\phi = 2,5,10,20\), respectively. For parameters \(\sigma \in \{0.90,0.63,0.45\}\), we use regression coefficients \((\theta_0, \boldsymbol{\theta}^{\top})^{\top} = (0.5,0.4,-0.3,0.3)^{\top}\), while for the higher-variance case \(\sigma = 1.50\), we adopt the more conservative setting \((0.4,0.25,-0.2,0.2)^{\top}\). Responses are simulated as \[w_i = \text{expit}(\theta_0 + \boldsymbol{x}_i^{\top} \boldsymbol{\theta} + \epsilon_i), \quad \epsilon_i \overset{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma^2).\]

  2. Heteroscedastic logit-normal regression model (S2): To ensure that the standard deviation \(\sigma_i\) remains within the range of 0.45 and 1,5, we model it as a log-linear function of the covariates: \[\log(\sigma_i) = \eta_0 + {\boldsymbol{x}}_i^{\top}\boldsymbol{\eta},\] where \((\eta_0, \boldsymbol{\eta}^{\top})^{\top} = (-0.2, 0.06, 0.06, 0.06)^{\top}\) is a fixed coefficient vector. Given that \({\boldsymbol{x}}_i \sim \mathcal{N}_3(\mathbf{0}, \boldsymbol{\Sigma})\), it follows that \(\log(\sigma_i) \sim \mathcal{N}(\eta_0, \boldsymbol{\eta}^{\top} \boldsymbol{\Sigma} \boldsymbol{\eta})\). Applying the empirical rule for the normal distribution, we approximate \[\log(\sigma_i) \in [\mu - 3\sigma, \mu + 3\sigma] \approx [-0.65,\;0.25],\] which implies \(\sigma_i \in [\exp(-0.65),\;\exp(0.25)] \approx [0.52,\;1.28] \subset [0.45,1.5]\). For the mean (\(\mu\)), we use regression coefficients \((\theta_0, \boldsymbol{\theta}^{\top})^{\top} = (0.4, 0.25, -0.2, 0.2)^{\top}\). The responses are then generated as \[w_i = \text{expit}(\theta_0 + {\boldsymbol{x}}_i^{\top} \boldsymbol{\theta}+ \epsilon_i),\] where \(\epsilon_i \overset{\text{i.i.d.}}{\sim} \mathcal{N}\left(0,\;\exp\{2(\eta_0 + {\boldsymbol{x}}_i^{\top} \boldsymbol{\eta})\} \right)\).

  3. Beta regression model with covariate-dependent mean (S3): To ensure that the conditional mean \(\mu_i\) remains well within the unit interval and avoids values near the boundaries (0 and 1), the regression coefficients are chosen to control the scale of the linear predictor. For precision parameters \(\phi \in \{5,\;10,\;20\}\), we set \((\theta_0, \boldsymbol{\theta}^{\top})^{\top} = (0.5,\;0.4,\;-0.3,\;0.3)^{\top}\). For \(\phi = 2\), a more conservative choice \((\theta_0, \boldsymbol{\theta}^{\top})^{\top} = (0.4,\;0.25,\;-0.2,\;0.2)^{\top}\) is used to mitigate the higher variability associated with lower precision. The responses are simulated from a beta distribution with parameters \((\mu_i \phi,\;(1 - \mu_i) \phi)\), where \[\mu_i = \text{expit}(\theta_0 + \boldsymbol{x}_i^{\top} \boldsymbol{\theta}).\]

  4. Beta regression model with covariate-dependent mean and dispersion (S4): In this scenario, the precision parameter \(\phi_i\) is modelled as a log-linear function of the covariates to ensure it remains within the range of 2 and 20: \[\log(\phi_i) = \eta_0 + {\boldsymbol{x}}_i^{\top} \boldsymbol{\eta},\] with fixed coefficients \((\eta_0, \boldsymbol{\eta}^{\top})^{\top} = (1.85, 0.15, 0.15, 0.15)^{\top}\). Given that \({\boldsymbol{x}}_i \sim \mathcal{N}_3(\mathbf{0}, \boldsymbol{\Sigma})\), it follows that \(\log(\phi_i) \sim \mathcal{N}(\eta_0, \boldsymbol{\eta}^{\top} \boldsymbol{\Sigma} \boldsymbol{\eta})\). Using the empirical rule for the normal distribution, we approximate \[\log(\phi_i) \in [\mu - 3\sigma,\;\mu + 3\sigma] \approx [0.74,\;2.96],\] which implies \(\phi_i \in [\exp(0.74),\;\exp(2.96)] \approx [2.1,\;19.3] \subset [2,20]\).

     The mean is modelled as \[\mu_i = \text{expit}(\theta_0 + {\boldsymbol{x}}_i^{\top} \boldsymbol{\theta}),\] with coefficients \((\theta_0, \boldsymbol{\theta}^{\top})^{\top} = (0.4, 0.25, -0.2, 0.2)^{\top}\). The responses are then drawn from a beta distribution with parameters \((\mu_i \phi_i,\;(1 - \mu_i) \phi_i)\).

  5. Model frameworks: The construction of conformal prediction intervals follows the models and scores introduced in Sections 2 and 3. For the logit-normal model (M1), the absolute raw residual is used as the non-conformity score, while for the heteroscedastic logit-normal model (M2), the absolute Pearson residual is employed. In the beta regression framework, two modelling strategies are evaluated: modelling with covariate-dependent \(\mu\) only under the absolute Pearson residual (M3(P)) and the absolute quantile residual (M3(Q)), and modelling with both covariate-dependent \(\mu\) and \(\phi\) under the absolute Pearson residual (M4(P)) and the absolute quantile residual (M4(Q)).

  6. Tolerance and parameter choice for the full conformal: In this simulation study, we set \(\varepsilon=10^{-4}\). To achieve this, for models M1 and M2, prediction intervals are constructed over the real line using a fixed grid with step size \(10^{-4}\). For models M3 and M4, the prediction interval is inside \((0,1)\). Therefore, we solve for the interval numerically with a tolerance of \(10^{-4}\), matching with the cases of M1 and M2, rather than constructing a initial grid. \(\rho_{\text{grid}}=3\) is set as the grid expansion parameter for the full conformal under logit-normal regression model, as is introduced in Section 3.4.

As a benchmark, we construct bootstrap-based prediction intervals approach from the work of [16], which relies on standardized weighted residuals derived from the score function of the beta regression model. In this technique, we consider the beta regression model based on covariate-dependent \(\mu\) and \(\phi\) and use the percentile quantiles to construct the interval. Again, to improve stability, 10,000 replications are used for \(n=50\) and 1,000 replications for other values of \(n\).

Simulation results↩︎

The results for scenarios S1 and S3, under \(\sigma = 0.63\) and \(\phi = 10\), respectively, together with results for S2 and S4, are presented in Table 1. It also includes the empirical results under corresponding settings. The rest of the simulation results are presented in Table 4 of Appendix 10. Table 1 reveals several key patterns regarding coverage accuracy and interval efficiency across different simulation scenarios and modelling frameworks. First, both split and full conformal prediction (CP) methods achieve empirical coverage rates close to the nominal 90% level in all scenarios, regardless of whether the model is correctly specified or misspecified. Notably, full conformal demonstrates superior calibration in small samples (\(n=50\)), more closely approximating the nominal level (e.g., 90.4% vs. 92.4% for M1 under S1). This aligns with theoretical expectations, as the full CP avoids the additional randomness introduced by sample splitting. However, as the sample size increases, the coverage differences between the split and full CP methods decrease. In some cases, the split CP performs comparably, suggesting that both methods are asymptotically well-calibrated under regular conditions.

The full CP produces substantially narrower intervals across all settings. This efficiency advantage is most pronounced in smaller samples (e.g., width 0.4946 vs. 0.6425 for M4(P) under S1, \(n=50\)) but persists in larger samples, without compromising coverage. This gain in efficiency arises because the full CP uses the entire dataset for both model fitting and residual calculation, avoiding the step of data splitting that reduces effective sample size in the split CP. This advantage is especially valuable in small-sample regimes, where the calibration set in split CP may be too small to reliably estimate quantiles, leading to over-conservative intervals. Even as the sample size increases, the full CP retains a moderate efficiency benefit due to its more stable and concentrated residual distribution.

Models incorporating heteroscedasticity (M2 and M4) generally produce wider prediction intervals than their simpler counterparts M1 and M3 when the additional variance or dispersion structure is unnecessary or misspecified (e.g., M2 presents wider intervals than M1 almost throughout under S1). This supports the idea that unnecessary heteroscedastic modelling can inflate interval width. When the heteroscedastic models are correctly specified, this width penalty is still visible in small samples (e.g., S4 at \(n=50\): M4 is wider than M3). This behaviour is consistent with the parsimony principle: although M2 and M4 are more flexible, they require estimation of additional variance or dispersion components, which can inflate the conformal score quantiles in finite samples. As the sample size increases, this estimation penalty diminishes, and the benefit of modelling heteroscedasticity becomes more apparent (e.g., S2 at \(n=1000\): M2 is narrower than M1).

A clear reduction in both split and full CP interval widths is observed across most scenarios as the sample size increases from \(n = 50\) to \(n = 1000\). This progressive narrowing performance reflects the improvement in parameter estimation precision by larger samples. For instance, under the scenario S3 using full conformal prediction, the interval width for model M3(P) decreases from 0.4912 at \(n = 50\) to 0.4615 at \(n = 1000\), while M4(P) exhibits an even more pronounced reduction from 0.5189 to 0.4626. The greatest relative decrease in width tends to occur between \(n = 50\) and \(n = 100\), highlighting the sensitivity of conformal methods to sample size in small-sample regimes. Importantly, these gains in interval efficiency are achieved without compromising marginal coverage, which remains close to the nominal level across all cases. This indicates that increased sample size enhances the stability and concentration of residual-based quantiles, therefore tightening prediction intervals while preserving their statistical validity.

The bootstrap method achieves empirical coverage close to the nominal 90% level in most scenarios and sample sizes. However, it shows mild under-coverage when \(n=50\), in contrast to the conformal methods, which remain close to or above the nominal level. This may reflect the reliance of the bootstrap procedure on model-based resampling approximations, which can underestimate predictive uncertainty in small samples. Since the bootstrap method is based on a beta regression model with covariate-dependent mean and precision, it is naturally aligned with S3 and S4. Under these beta-generated scenarios, it is competitive in terms of interval width; for example, when \(n=500\) under S4, its width is 0.5912, compared with 0.5934 for M4(Q). In contrast, under S1 and S2, where the data are generated from logit-normal models, the bootstrap method is applied under model misspecification. In these cases, its intervals are generally slightly wider than those from the proposed conformal methods; for instance, when \(n=1000\) under S1, its width is 0.4372, compared with 0.4358 for M3(Q).

Table 3 in Appendix 10 reports the Monte Carlo errors under 1,000 iterations corresponding to the empirical coverage and width in Table 1. The Monte Carlo standard errors for coverage are uniformly small across all frameworks, scenarios, and sample sizes, ranging from approximately 0.007 to 0.010, confirming that the empirical coverage estimates in Table 1 are stable and reliable. The standard errors for interval width are similarly small and decrease with increasing \(n\), reflecting the improved concentration in larger samples. Notably, the standard errors for width are larger under heteroscedastic models (M2, M4), especially at \(n=50\), consistent with the greater variability in interval construction when dispersion is estimated from limited data. For the bootstrap, errors are of similar magnitude to those of the conformal methods, indicating that the simulation estimates are comparable. Other Monte Carlo errors are reported in Table 5 of the Appendix.

6 Applications↩︎

In this section, we analyze a dataset containing information on college-aged women (18–25 years old) at Brigham Young University, originally collected by [31] and later analyzed by [32]. One observation with an implausible age entry (age = 1) was removed from the original dataset of \(184\) individuals, resulting in a final sample size of \(n=183\). This dataset includes the body fat percentage measured by “Siri’s equation” [33] defined as \[\text{Body Fat Percentage} = \left( \frac{495}{\text{Body Density}} \right) - 450.\] While this equation can theoretically yield values at 0 or 1 (after scaling), such extremes are not biologically plausible. In this dataset, body fat percentages range from 0.0747 to 0.3849, well within the open interval \((0,1)\), which supports the use of beta and transformation regression models. In addition to BMI (Body Mass Index: calculated as weight in kilograms divided by the square of height in meters), weight, height, and age, the dataset contains 14 standardized circumference measurement (converted to centimeters) for specific body regions. Given the modest sample size, we select a subset of eight covariates to enhance model stability and interpretability. For each covariate, we report a brief description along with summary statistics in the form \((\min,\, \operatorname{mean},\, \max)\): BMI - \((15.80,\, 21.60,\, 31.14)\), Neck (the minimal circumference of the neck, measured perpendicular to its long axis) - \((26.00,\, 31.48,\, 36.00)\), Chest (circumference measured in a horizontal plane at the level of the sixth rib, taken at the end of a normal expiration) - \((43.00,\, 85.02 ,\, 100.00)\), Hips (the maximum horizontal circumference around the buttocks) - \((82.50,\, 96.93,\, 115.00)\), Waist (the minimal horizontal measurement, taken at the end of a normal expiration) - \((58.00,\, 69.46,\, 96.00)\), Forearm (the maximum circumference measured perpendicular to the long axis of the forearm) - \((20.00,\, 23.49,\, 29.00)\), PThigh (Proximal Thigh: circumference immediately below the gluteal fold) - \((45.00,\, 57.38,\,69.50)\), and Wrist (circumference measured perpendicular to the long axis of the forearm) - \((13.50,\, 15.65,\, 19.00)\). We focus on this subset of covariates for the following reasons: BMI, Hips, and Waist have been chosen as preferable variables for body fat in multivariable linear models [32]. Chest and Neck measurements provide information about upper body size and can enhance model accuracy, particularly when fat is not primarily distributed around the Hips and Waist. Notably, recent studies have shown that Neck circumference is a robust marker of body adiposity, especially in women [34]. Forearm and Wrist measurements reflect peripheral limb size and offer insights into body frame and muscle development. Finally, PThigh contributes additional information about the lower body and helps balance the representation of upper and lower body characteristics.

To implement conformal prediction, the dataset is randomly split into 90% for interval construction and 10% reserved as a test set (20 samples) for final evaluation. Within the 90% construction set, split conformal methods further divide the data equally into training and calibration subsets. Table 2 summarizes the empirical coverage and average interval width obtained under the modelling frameworks considered. Figure 3 presents a visualization of prediction intervals among all model frameworks. Each panel corresponds to a different model under either the split or full conformal prediction setting, with test points plotted on the x-axis and predicted intervals shown vertically. Points where the true body fat percentage falls outside the predicted interval are highlighted. Those green or red points are predictions from the transformation regression model or the beta regression model based on the training set (split CP) or the training set plus the calibration set (full CP). Across all settings, conformal prediction intervals consistently achieve the nominal 90% coverage. As observed in both the simulation study and the real data application, full conformal methods yield slightly narrower intervals than their split counterparts. Introducing heteroscedasticity into either the transformation or beta regression models increases the width of the prediction intervals, as expected. Notably, the beta regression models outperform the transformation models by achieving narrower intervals while maintaining the same level of coverage. Among all configurations, the beta regression model using quantile residuals produces the narrowest intervals.

4pt

Notably, the missed observations are concentrated around test points 12 or 18 across the different model frameworks. This pattern is robust to both the model specification (transformation vs. beta regression) and the CP methodology (split vs. full), suggesting these test samples may represent difficult-to-predict individuals rather than isolated model deficiencies. These observations are not located near the extremes of the body fat distribution, indicating that their prediction difficulty is not due to marginal outlyingness in the response. Instead, the missed observations may be conditionally unusual given their covariate profiles. For example, atypical body shape profiles, such as high BMI but low limb circumference, or vice versa, which could result in covariate combinations not well-represented in the dataset used for model construction. This also suggests a connection with multivariable modelling: the non-conformity of these observations may arise from unusual dependence patterns among the body-measurement covariates, rather than from extreme body-fat values alone.

While the bootstrap provides reasonable efficiency under correct model specification, it depends critically on the accuracy of both the mean and dispersion submodels, and relies on asymptotic approximations. In contrast, the conformal procedures considered here provide finite-sample marginal coverage guarantees under exchangeability, without requiring correct specification of the working regression model. Notably, the full conformal methods applied to the beta regression framework often yield narrower intervals than the bootstrap, demonstrating both statistical efficiency and robustness. These findings underscore the advantages of full CP, particularly when using quantile-based residuals. To assess the consistency of interval estimates across CP models, we perform a sensitivity analysis based on the union and intersection of the prediction intervals. For each test observation, the union interval is obtained by taking the smallest lower endpoint and the largest upper endpoint among all fitted CP intervals. It therefore gives a conservative summary of the range supported by at least one model. The intersection interval is obtained by taking the largest lower endpoint and the smallest upper endpoint among all fitted CP intervals. It represents the region jointly supported by all models, and is empty if the largest lower endpoint exceeds the smallest upper endpoint. The absence of empty intersections confirms a strong level of agreement among the models considered. We report their empirical coverage and average width as a sensitivity analysis of the stability of the prediction intervals across model frameworks and nonconformity scores.

7 Discussion↩︎

We study the construction of conformal prediction intervals for continuous bounded outcomes using model-informed non-conformity measures that account for the inherent heteroscedasticity of these data. Because the scores are derived from the fitted model through the probability integral transform, they remain interpretable and, unlike symmetric Pearson-type residuals, they respect the support of the response and account for the boundary-induced asymmetry of the conditional distribution, so that the resulting intervals stay within the admissible range. We establish finite-sample marginal validity for both split and full conformal prediction under exchangeability and a mild no-ties condition, together with asymptotic conditional validity under additional regularity conditions. Our simulation study supports these results, with all methods attaining coverage at or near the nominal level, including under model misspecification, and remaining competitive with, or more stable than, a bootstrap benchmark, in particular avoiding the mild under-coverage of the bootstrap in small samples.

A second practical consideration is the trade-off between the two algorithms. Split conformal prediction fits the model once and is essentially instantaneous, making it well suited to large datasets and real-time use, but the data split widens the intervals, especially in small samples. Full conformal prediction uses the entire sample for both fitting and calibration and consistently produces narrower intervals, at the cost of repeated refitting and noticeably higher and more variable runtimes. The adaptive interval-finding routine we employ mitigates this cost, and developing faster approximations to full conformal in this setting is a worthwhile direction. In practice, we recommend the quantile-residual score as the default, given its boundary-respecting behaviour and theoretical interpretability, paired with full conformal when computational resources allow and split conformal when speed or scale is the priority.

We focus on logit-normal regression models and beta regression because of their wide applicability, interpretability, and well-established software implementations. Alternative models and algorithms [10], [12], [14] could equally be employed and integrated with the methodology introduced here, and different quantile definitions or normalised conformal prediction methods could be used to further accommodate heteroscedasticity or other features of the data. Other potential extensions include the prediction of multivariate bounded outcomes, conformal methods for zero- and/or one-inflated beta regression models, which pose an additional challenge due to the mixed discrete–continuous nature of the response, and the extension to hierarchical data structures such as spatial, longitudinal, or dynamic settings.

In our real-data application, we also performed a sensitivity analysis comparing the union and intersection of the conformal prediction intervals across methodological frameworks. A natural question follows: since we use four models for bounded responses together with various non-conformity measures, could their predictions be combined, rather than relying on a single model or on the union or intersection of their intervals? One option is to stack the models into an ensemble and aggregate the corresponding prediction sets (see Section 10.4 of [35] for a review of stacking methods, and [36] for a recent proposal), although this requires further analysis to guarantee the validity of the resulting intervals. The real-data application also points to a connection with multivariate modelling. The two non-conforming points we detected did not have marginally extreme responses; rather, they corresponded to atypical joint covariate profiles (for example, high BMI together with low limb circumference). This suggests that non-conformity can be driven by the dependence structure of the data rather than by the marginal distribution of the response alone, which is in line with findings that atypical observations are often best identified through explicit modelling of dependence rather than marginal behaviour [37].

8 Acknowledgements↩︎

We thank Prof. Vladimir Vovk for helpful discussions about this project that helped us add precision to the results presented here. We thank three reviewers and the Action Editor for their constructive comments provided.

9 Proof of technical results↩︎

9.1 Lemmas↩︎

We now prove that the quantile-based non-conformity scores 5 are symmetric and almost surely free of ties.

The following lemma shows that the empirical quantile convergence is by C1-C3 for continuous regression models, covering beta regression and logit-normal regression models coupled with quantile residuals \(r_i^Q\) in 5 .

9.2 Proof of Proposition [prop:marginal]↩︎

Let \(\alpha\in(0,1)\) be the fixed significance level, and define rank indices \(k := \lceil (1-\alpha)(n+1)\rceil\) for full CP and \(k_{\mathrm{cal}} := \lceil (1-\alpha)(n_{cal}+1)\rceil\) for split CP.

9.2.0.1 Full CP

Since \(T:(0,1)\to \mathcal{Y}\) is a strictly increasing mapping by Assumption A1, it has a well-defined inverse \(T^{-1}:\mathcal{Y}\to(0,1)\). By construction, Algorithm 2 searched over \(y\in\mathcal{Y}\), so \(\widetilde{\mathcal{C}}^{\text{full}}_{1-\alpha}({\boldsymbol{X}}_{\text{new}})\subseteq\mathcal{Y}\). Since \(Y_{\text{new}}=T(W_{\text{new}})\) and \(T^{-1}\) is strictly increasing, \(\{W_{\text{new}}\in\mathcal{C}^{\text{full}}_{1-\alpha}({\boldsymbol{X}}_{\text{new}})\}\) if and only if \(\{Y_{\text{new}}\in\widetilde{\mathcal{C}}^{\text{full}}_{1-\alpha}({\boldsymbol{X}}_{\text{new}})\}\).

Fix the candidate value at the true response, \(y=Y_{new}\). We show that under assumptions A1-A3, the proof then follows by standard CP theory [35]. For notational simplicity, in the remainder of the proof we use the index \(n+1\) interchangeably with “new” to denote the test observation. Let \(\varrho := \sum_{j=1}^{n+1}\mathbf{1}\{s_j^{(y)} \le s_{new}^{(y)}\}\) be the rank of the test score among the augmented scores. The scores \((s_1^{(y)},\dots,s_{new}^{(y)})\) are exchangeable and almost surely tie-free from Lemma 1, so \(\varrho\) is uniformly distributed on \(\{1,\dots,n+1\}\). By construction of Algorithm 2, the threshold \(q_{1-\alpha}^{(y)}\) is the \(k\)-th smallest of \(\{s_1^{(y)},\dots,s_n^{(y)}\}\) (with the convention \(q_{1-\alpha}^{(y)}=+\infty\) when \(k>n\)), and \(Y_{new}\in\widetilde{\mathcal{C}}^{\text{full}}_{1-\alpha}({\boldsymbol{X}}_{new})\) if and only if \(s_{new}^{(y)}\le q_{1-\alpha}^{(y)}\). In the absence of ties, \(\{s_{new}^{(y)}\le q_{1-\alpha}^{(y)}\} = \{\varrho \le k\}\), so \[\Pr\!\left(Y_{new}\in\widetilde{\mathcal{C}}^{\text{full}}_{1-\alpha}({\boldsymbol{X}}_{new})\right) \;=\; \Pr(\varrho \le k) \;=\; \frac{\min(k,\,n+1)}{n+1} \;=\; \frac{\lceil (1-\alpha)(n+1)\rceil}{n+1}.\] The elementary bounds \((1-\alpha)(n+1)\le \lceil (1-\alpha)(n+1)\rceil < (1-\alpha)(n+1)+1\) lead to \[1-\alpha \;\le\; \Pr\!\left(Y_{new}\in\widetilde{\mathcal{C}}^{\text{full}}_{1-\alpha}({\boldsymbol{X}}_{new})\right) \;\le\; 1-\alpha+\frac{1}{n+1}.\]

9.2.0.2 Split CP

By construction in Algorithm 1, the model is fitted once, on \(\{{\boldsymbol{Z}}_i:i\in\mathcal{I}_{\text{train}}\}\), so the fitted statistic \(\widehat{\boldsymbol{\zeta}}\), and hence the score function \(z\mapsto s(z;\widehat{\mathcal{M}}(\widehat{\boldsymbol{\zeta}}))\), is a fixed measurable function once \(\mathcal{I}_{\text{train}}\) is given. Conditional on \(\mathcal{I}_{\text{train}}\), the calibration responses and the test point, \(\{{\boldsymbol{Z}}_i:i\in\mathcal{I}_{\text{cal}}\}\cup\{{\boldsymbol{Z}}_{new}\}\), remain exchangeable and are independent of the training fold; applying the common fixed score function gives exchangeable scores \(\{s_i:i\in\mathcal{I}_{\text{cal}}\}\cup\{s_{new}\}\), almost surely tie-free by the assumptions of the Proposition. Repeating the rank argument aforementioned with \(n_{cal}\) in place of \(n\), and \(k_{cal}\) in place of \(k\) gives, conditionally on \(\mathcal{I}_{\text{train}}\), \[1-\alpha \;\le\; \Pr\!\left(Y_{new}\in\widetilde{\mathcal{C}}^{\text{split}}_{1-\alpha}({\boldsymbol{X}}_{new})\,\middle|\,\mathcal{I}_{\text{train}}\right) \;\le\; 1-\alpha+\frac{1}{n_{cal}+1}.\]

Taking expectations over \(\mathcal{I}_{\text{train}}\) and applying the law of total expectation yields the unconditional bound \(1-\alpha \;\le\; \Pr\!\left(Y_{new}\in\widetilde{\mathcal{C}}^{\text{split}}_{1-\alpha}({\boldsymbol{X}}_{new})\right) \;\le\; 1-\alpha+\frac{1}{n_{cal}+1}\), which completes the proof since the equivalence \(\{W_{\text{new}}\in\mathcal{C}^{\text{split}}_{1-\alpha}({\boldsymbol{X}}_{\text{new}})\} \iff \{Y_{\text{new}}\in\widetilde{\mathcal{C}}^{\text{split}}_{1-\alpha}({\boldsymbol{X}}_{\text{new}})\}\) holds as before.

9.3 Proof of Proposition [prop:conditional]↩︎

We prove the result for full CP; the argument for split CP is analogous, with the calibration set \(\mathcal{J}_{\text{cal}}\) replacing \(\{1,\ldots,n\}\). As is verified in the Proof of Proposition 1, the equivalence \(\{W_{\text{new}}\in\mathcal{C}^{\text{full}}_{1-\alpha}({\boldsymbol{X}}_{\text{new}})\} \iff \{Y_{\text{new}}\in\widetilde{\mathcal{C}}^{\text{full}}_{1-\alpha}({\boldsymbol{X}}_{\text{new}})\}\) holds.

Define the conditional coverage gap \[\Delta_n({\boldsymbol{x}}_{new}) := \Pr\!\left(Y_{new} \in \widetilde{\mathcal{C}}^{\text{full}}_{1-\alpha}({\boldsymbol{X}}_{new}) \mid {\boldsymbol{X}}_{new} = {\boldsymbol{x}}_{new}\right) - (1-\alpha).\] We show \(\Delta_n({\boldsymbol{x}}_{new}) = o_p(1)\). Define \(s_{\text{new}} := S(Y_{\text{new}};\widehat{\boldsymbol{\zeta}},{\boldsymbol{X}}_{\text{new}})\) as the non-conformity score for the new observation based on the fitted model \(\widehat{\mathcal{M}}\), and \(V_{\text{new}}:=S(Y_{\text{new}};\boldsymbol{\zeta}^*,{\boldsymbol{X}}_{\text{new}})\) as the oracle non-conformity score for the new observation on the limit model \(\mathcal{M}^{(*)}\). Then, define \(R_{new} := s_{new} - V_{new}\). By Condition C1, \(R_{new} = o_p(1)\). Denote \(G_\star(v) := \Pr(V_{new} \leq v \mid {\boldsymbol{X}}_{new} = {\boldsymbol{x}}_{new})\), and define \(q_{1-\alpha}^* := G_\star^{-1}(1-\alpha)\). By construction of Algorithm 2, \(Y_{new} \in \widetilde{\mathcal{C}}^{\text{full}}_{1-\alpha}({\boldsymbol{X}}_{new})\) if and only if \(s_{new} \leq \widehat{q}_{1-\alpha}\), where \(\widehat{q}_{1-\alpha}\) is the empirical \((1-\alpha)\)-quantile of \(\{s_i : i = 1,\ldots,n\}\). Adding and subtracting \(\Pr(s_{new} \leq q_{1-\alpha}^* \mid {\boldsymbol{x}}_{new})\), \[\begin{align} \Delta_n({\boldsymbol{x}}_{new}) &= \underbrace{\Pr(s_{new} \leq \widehat{q}_{1-\alpha} \mid {\boldsymbol{x}}_{new}) - \Pr(s_{new} \leq q_{1-\alpha}^* \mid {\boldsymbol{x}}_{new})}_{(A)} + \underbrace{\Pr(s_{new} \leq q_{1-\alpha}^* \mid {\boldsymbol{x}}_{new}) - (1-\alpha)}_{(B)}. \label{eq:decomp} \end{align}\tag{7}\]

9.3.0.1 Term (B) vanishes.

Write \(s_{new} = V_{new} + R_{new}\) with \(R_{new} = o_p(1)\). For any \(\delta > 0\) we bound \(\Pr(s_{new} \leq q_{1-\alpha}^* \mid {\boldsymbol{x}}_{new})\) from above and below.

For the upper bound, partitioning on \(\{|R_{new}| > \delta\}\) and \(\{|R_{new}| \leq \delta\}\) and noting that on \(\{|R_{new}| \leq \delta\}\) the inequality \(s_{new} \leq q_{1-\alpha}^*\) together with \(V_{new} = s_{new} - R_{new} \leq s_{new} + \delta\) gives \(V_{new} \leq q_{1-\alpha}^* + \delta\),

\[\{s_{new} \leq q_{1-\alpha}^*\} \cap \{|R_{new}| \leq \delta\} \subseteq \{V_{new} \leq q_{1-\alpha}^* + \delta\}. \label{eq:subsets}\tag{8}\]

Using \(\Pr(B \cap C) \geq \Pr(B) - \Pr(C^c)\) with \(B = \{s_{new} \leq q_{1-\alpha}^*\}\) and \(C = \{|R_{new}| \leq \delta\}\) together with 8 ,

\[\begin{align} \Pr\!\left(s_{new} \leq q_{1-\alpha}^* \mid {\boldsymbol{x}}_{new}\right) &\leq \Pr\!\left(V_{new} \leq q_{1-\alpha}^* + \delta \mid {\boldsymbol{x}}_{new}\right) + \Pr(|R_{new}| > \delta) \\ &= G_\star(q_{1-\alpha}^* + \delta) + \Pr(|R_{new}| > \delta). \end{align}\]

Similarly, for the lower bound, on \(\{|R_{new}| \leq \delta\}\) the inequality \(V_{new} \leq q_{1-\alpha}^* - \delta\) implies \(s_{new} = V_{new} + R_{new} \leq q_{1-\alpha}^*\), so that \[\{V_{new} \leq q_{1-\alpha}^* - \delta\} \cap \{|R_{new}| \leq \delta\} \subseteq \{s_{new} \leq q_{1-\alpha}^*\}.\] Now, set \(B = \{V_{new} \leq q_{1-\alpha}^* - \delta\}\) and \(C = \{|R_{new}| \leq \delta\}\), \[\begin{align} \Pr\!\left(s_{new} \leq q_{1-\alpha}^* \mid {\boldsymbol{x}}_{new}\right) &\geq \Pr\!\left(V_{new} \leq q_{1-\alpha}^* - \delta \mid {\boldsymbol{x}}_{new}\right) - \Pr(|R_{new}| > \delta) \\ &= G_\star(q_{1-\alpha}^* - \delta) - \Pr(|R_{new}| > \delta). \end{align}\]

Combining the two bounds, \[G_\star(q_{1-\alpha}^* - \delta) - \Pr(|R_{new}| > \delta) \leq \Pr\!\left(s_{new} \leq q_{1-\alpha}^* \mid {\boldsymbol{x}}_{new}\right) \leq G_\star(q_{1-\alpha}^* + \delta) + \Pr(|R_{new}| > \delta).\] Fix \(\eta > 0\). By continuity of \(G_\star\) at \(q_{1-\alpha}^*\) (Condition C2), choose \(\delta > 0\) such that \(G_\star(q_{1-\alpha}^* + \delta) - G_\star(q_{1-\alpha}^* - \delta) \leq \eta/2\), and then \(N\) such that \(\Pr(|R_{new}| > \delta) \leq \eta/2\) for all \(n > N\), which is possible since \(R_{new} = o_p(1)\). For \(n > N\), both bounds lie within \(\eta\) of \(G_\star(q_{1-\alpha}^*) = 1 - \alpha\), so that \[\left|\Pr\!\left(s_{new} \leq q_{1-\alpha}^* \mid {\boldsymbol{x}}_{new}\right) - (1-\alpha)\right| \leq \eta.\] As \(\eta > 0\) was arbitrary, \(\Pr(s_{new} \leq q_{1-\alpha}^* \mid {\boldsymbol{x}}_{new}) = 1 - \alpha + o(1)\), so that \((B) = o(1)\).

9.3.0.2 Term (A) vanishes.

By the indicator identity \(|\mathbf{1}\{t \leq a\} - \mathbf{1}\{t \leq b\}| = \mathbf{1}\{\min(a,b) < t \leq \max(a,b)\}\), \[|(A)| \leq \Pr\!\left(\min(\widehat{q}_{1-\alpha}, q_{1-\alpha}^*) < s_{new} \leq \max(\widehat{q}_{1-\alpha}, q_{1-\alpha}^*) \mid {\boldsymbol{x}}_{new}\right).\] Fix \(\epsilon > 0\). On \(\{|\widehat{q}_{1-\alpha} - q_{1-\alpha}^*| \leq \epsilon\}\) this event is contained in \(\{q_{1-\alpha}^* - \epsilon < s_{new} \leq q_{1-\alpha}^* + \epsilon\}\), so \[|(A)| \leq \Pr\!\left(q_{1-\alpha}^* - \epsilon < s_{new} \leq q_{1-\alpha}^* + \epsilon \mid {\boldsymbol{x}}_{new}\right) + \Pr\!\left(|\widehat{q}_{1-\alpha} - q_{1-\alpha}^*| > \epsilon\right).\] Fix \(\eta > 0\). By continuity of \(G_\star\), choose \(\epsilon, \delta > 0\) with \(G_\star(q_{1-\alpha}^* + \epsilon + \delta) - G_\star(q_{1-\alpha}^* - \epsilon - \delta) \leq \eta/3\); then, arguing via \(\{|R_{new}| \leq \delta\}\) as in Term (B), \[\Pr\!\left(q_{1-\alpha}^* - \epsilon < s_{new} \leq q_{1-\alpha}^* + \epsilon \mid {\boldsymbol{x}}_{new}\right) \leq \frac{\eta}{3} + \Pr(|R_{new}| > \delta).\] Since \(R_{new} = o_p(1)\) and \(\widehat{q}_{1-\alpha} \xrightarrow{\Pr} q_{1-\alpha}^*\) (Lemma 2), there is \(N\) with \(\Pr(|R_{new}| > \delta) \leq \eta/3\) and \(\Pr(|\widehat{q}_{1-\alpha} - q_{1-\alpha}^*| > \epsilon) \leq \eta/3\) for \(n > N\). Hence \(|(A)| \leq \eta\) for \(n > N\), and as \(\eta\) was arbitrary, \(|(A)| = o(1)\).

Combining the bounds on \((A)\) and \((B)\) gives \(\Delta_n({\boldsymbol{x}}_{new}) = o_p(1)\), which establishes asymptotic conditional validity for full CP.

10 Additional simulation results↩︎

During the simulation process for \(n = 50\) with Monte Carlo iteration \(N = 10{,}000\) under the split conformal prediction (CP) setting, we encountered occasional failures in model fitting due to numerical instability, particularly in heteroscedastic models. When such cases occurred, the corresponding random seed was discarded and replaced with a new one to ensure the total number of valid replications remained at \(N = 10,000\). The affected random seeds, which caused model-fitting failures and were excluded from the simulation, are listed below:

• S1-M2: \(\sigma = 1.5\): 8608, \(\sigma = 0.9\): 8608, \(\sigma = 0.63\): 8608, \(\sigma = 0.45\): 8608.

• S3-M2: \(\phi = 2\): 3403, \(\phi = 5\): 2807.

• S2-M2: 3311.

For the bootstrap-based prediction intervals, occasional numerical failures were also observed when \(n=50\). The skipped bootstrap seeds were as follows: S1-Bootstrap: 6952; S2-Bootstrap: 9437; S3-Bootstrap: 3052, 6087; and S4-Bootstrap: 3635, 8247, 9180, 9644. These issues were limited to a very small number of seeds and did not materially affect the overall validity or stability of the simulation results.

All non-bootstrap simulations were conducted in R version 4.3.1 (2023-06-16) on a machine running macOS (Darwin), equipped with an Apple M2 Pro processor and 16 GB of RAM, using a single CPU core. The bootstrap-based simulations, which required substantially greater computational resources, were carried out on the CREATE high-performance computing cluster at King’s College London using Slurm batch jobs in R version 4.5.1. For these bootstrap jobs, we requested 1 node, 1 task, 8 CPU cores, and 32 GB of memory from the CPU partition.

Table 6 presents the corresponding average CPU time (in seconds) and its standard deviation under 100 iterations (for non-bootstrap simulations) and 1,000 iterations (for bootstrap simulations) across the same simulation scenarios and modelling frameworks as Table 1. We can observe that the split CP exhibits near-instantaneous computation times across all scenarios, averaging between 0.001 and 0.018 seconds regardless of sample size or model complexity. This efficiency stems from its working regime: a single model fit on the training subset followed by straightforward quantile calculations on the calibration set. Standard deviations remain negligible (\(\leq 0.007\) seconds), which indicates consistent runtime stability. This makes the split method ideal for real-time applications or large-scale iterative analyses where computational resources are constrained. In contrast, the full CP exhibits substantially higher computational costs due to its need to refit models for each candidate point. Additionally, the full CP presents significant runtime variability (with a standard deviation up to 0.119 seconds). This volatility seems to stem from algorithmic sensitivity to random covariate effects and optimization convergence in dispersion parameters. Despite these computational demands, the full CP’s more efficient intervals remain an advantage. Practitioners must therefore weigh the critical trade-off: split conformal for computational efficiency when chasing for the speed versus full conformal for statistical precision when resources permit. We also observe that generally, as \(\phi\) increases (or \(\sigma\) decreases), the computational time for the full CP decreases, which may be because the model is easier to fit under smaller variability. The computational costs for other scenarios are presented in Table 7.

11 Empirical evidence for the application↩︎

As further empirical evidence motivating the introduction of heteroscedasticity into CP models, we present the residual-versus-fitted plots using the training set plus the calibration set for both the logit-normal model and the beta regression model in Figure 4. These plots reveal the non-constant variance patterns that can support the introduction of heteroscedasticity in the construction of conformal and bootstrap prediction intervals. For the raw residual plot in the logit-normal regression model, the residuals are approximately centered around zero, but there is a funnel-shaped pattern suggesting that the spread of residuals may vary with the level of the fitted values. For the Pearson residual in the beta regression model, the pattern is more pronounced here, with increased residual dispersion observed at higher fitted values. These diagnostic plots therefore provide visual signal for the need of heteroscedastic conformal prediction models in our analysis.

References↩︎