Abstract
The pintervals package aims to provide a unified framework for constructing prediction intervals and calibrating predictions in a model-agnostic setting using set-aside calibration data. It comprises routines to construct inductive conformal as well as parametric and bootstrapped prediction intervals from any model that outputs point predictions. Several R packages and functions already exist for constructing prediction intervals, but they often focus on specific modeling frameworks or types of predictions, or require manual customization for different models or applications. By providing a consistent interface for a variety of prediction interval construction approaches (all model-agnostic), pintervals allows researchers to apply and compare them across different modeling frameworks and applications.Predictive models and forecasting have become increasingly important across a wide range of scientific disciplines, industry, and policy-making. Models ranging from simple linear regression to complex machine learning algorithms are now routinely used to generate predictions for a variety of outcomes, from economic indicators and public health metrics to social and political phenomena, climate modeling, and energy use (see for instance Hegre et al. 2017; Petropoulos et al. 2022; Lyu et al. 2021; Klyuev et al. 2022; Kochkov et al. 2024; Gong et al. 2022). In the social sciences, in particular, this reflects a growing recognition that forecasting can complement traditional causal inference by providing actionable, forward-looking insights to policymakers, civil society, and researchers alike (Nanlohy et al. 2017).
Most such models, however, provide only point predictions, which do not convey their inherent uncertainty. Yet, in many applications, understanding the uncertainty around point predictions is crucial for effective decision-making, as it provides a more complete picture of the risks and potential outcomes associated with different scenarios (Raftery 2016). For instance, in conflict forecasting, it may be more important for stakeholders, in order to effectively prioritize their resources, to know the likelihood that a conflict exceeds a certain threshold of fatalities, rather than just the expected number of fatalities (e.g., Hegre et al. 2024). Similarly, the width of a prediction interval informs decision-makers about the level of confidence they may place in a particular point prediction, which is essential for risk assessment and planning.
This limitation has led to a growing consensus that point predictions alone are not sufficient for effective communication with stakeholders or for providing actionable support to decision-makers, especially in high-stakes domains such as conflict forecasting. When forecasts inform real-world interventions, it is not enough to know what the most likely outcome is; decision-makers must also understand how uncertain the outcome is, and how much confidence they can place in the prediction (Raftery 2016; Hegre et al. 2024; Brandt et al. 2014; Chadefaux 2017). For example, a forecast suggesting a moderate risk of violence may warrant very different policy responses depending on whether the associated uncertainty is low or high. Transparent communication of predictive uncertainty can therefore help ensure that responses are proportionate to the level of risk, improving both reliability and accountability.
To quantify the uncertainty around predictions, researchers commonly turn to prediction intervals. A variety of techniques exist for constructing such intervals, including parametric methods (i.e., assuming a specific distribution for the prediction errors), bootstrapping (resampling errors from a calibration set), quantile regression, and Bayesian credible intervals. Each approach has strengths and limitations. Parametric intervals, when available analytically or approximately, may be easy to compute but rely on strong distributional or asymptotic assumptions. Bootstrap methods are flexible, but they assume that the distribution of errors is constant across the data, and they tend to underestimate variance (due to sampling from an empirical distribution function rather than the true cumulative distribution function). Quantile regression provides asymmetric intervals but requires careful specification and is model-specific. Bayesian methods allow for full posterior distributions but are also model-specific and require careful prior specification (for a review of these methods, see Tian et al. 2022; Zhang et al. 2020; Hesterberg 2015). Furthermore, Bayesian uncertainty is epistemic rather than aleatory, so credible sets cannot in general be interpreted as frequentist confidence sets. Bayesian posterior inference is also subject to the false confidence theorem (Balch et al. 2019; Martin 2019; Carmichael and Williams 2018): there always exist sets that do not contain the true value-to-be-predicted, yet are assigned arbitrarily large posterior probability, arbitrarily often over repeated sampling of data sets.
The limitations of these standard methods may result in unreliable/misleading uncertainty quantification, particularly in complex or heterogeneous data settings common in social science applications. There has thus been a demand for more-flexible, model-agnostic approaches to uncertainty quantification that offer verifiable reliability guarantees, and that can be applied across a wide range of predictive models and data types. To meet this demand, conformal prediction (CP) is proving to be a viable approach to uncertainty quantification.
Under the assumption of exchangeable data observations1, CP produces prediction sets with guaranteed coverage for any fixed sample size, at any chosen confidence level. A key advantage of CP is its model-agnostic nature: it can be applied on top of any prediction model, from simple linear regression models to complex machine learning algorithms, without requiring model-specific modifications (Shafer and Vovk 2008; Vovk et al. 2022). Recent extensions, including clustered CP [CCP; Ding et al. (2023); Hjort et al. (2024)], distance-weighted CP [DWCP; Hjort et al. (2025)], and bin-conditional CP [BCCP; Randahl et al. (2025)], have further enhanced its applicability by addressing local calibration challenges in settings with skewed or imbalanced outcome distributions.
This article serves as an introduction to the pintervals package (Randahl and Williams 2025), and as an expository on the theoretical and computational foundations of the supported methods for uncertainty quantification. We also showcase the package using panel data on county election votes by party in the United States presidential election of 2024. The package is designed to be flexible and extensible, allowing users to apply it to a wide range of predictive modeling tasks in the social sciences and beyond.
Several R packages and functions already exist for constructing
prediction intervals. For example, the forecast package
(Hyndman et al. 2025) provides methods for
time series forecasting, while the tidymodels and
mlr3 packages (Kuhn and Wickham
2020; Lang et al. 2019) offer tools for machine learning
workflows. Parametric prediction intervals can be constructed using the
stats package [for linear and generalized linear
regression models; R Core Team (2025)], or
by manually computing intervals (when analytically available) based on
standard distribution functions (e.g., qt for
t-distribution, qnorm for normal distribution, etc.).
Bootstrapped prediction intervals can be constructed using the
boot package [namely, it provides functions for
resampling; Angelo Canty and B. D. Ripley
(2024)], or by manually implementing bootstrapping procedures
using the stats package. Additionally, bias-corrected
and accelerated (BCa) bootstrap confidence intervals can be obtained
using the bcaboot package (Efron
and Narasimhan 2020), which provides a second-order accurate
alternative to standard bootstrap intervals. CP intervals can also be
obtained through a range of packages, such as the
conformalbayes package (McCartan
2025) for Bayesian CP, the conformalForecast
package (Wang and Hyndman 2025) for time
series forecasting, the LCP package (Guan 2021) for localized CP, and the
predictset package (Coverdale
2026) for constructing CP sets directly from a model-based
workflow.
Beyond these packages, there are also various tutorials on how to obtain CP intervals manually using functions available in packages such as tidymodels (Kuhn and Wickham 2020), and marginaleffects (Arel-Bundock et al. 2025; Arel-Bundock 2025).
Beyond R, tools for constructing prediction intervals and CP sets are available across several statistical computing environments. In Python, the MAPIE library (Cordier et al. 2023) is a widely used, scikit-learn-compatible package for model-agnostic uncertainty quantification that implements a range of CP methods for regression, classification, and time series, and is compatible with models from scikit-learn, TensorFlow, and PyTorch. The crepes package (Boström 2022, 2024) is another Python package for CP that implements standard, normalized, and Mondrian conformal regressors and predictive systems, and is fully model-agnostic. In Julia, the ConformalPrediction.jl package (Altmeyer 2023) implements CP for both regression and classification tasks and is designed to work with any supervised model trained in the MLJ.jl machine learning framework (Blaom et al. 2020).
While existing packages offer valuable tools for constructing prediction intervals, they generally lack a unified, model-agnostic interface. Most are tied to specific modeling frameworks or require extensive manual adjustments to accommodate different models or applications, making it challenging for researchers to apply and compare them consistently across datasets and prediction tasks. The pintervals package addresses this limitation by providing a coherent and model-agnostic framework for generating prediction intervals. It implements several widely used modern and classical approaches, including inductive CP (with subcategories such as weighted, clustered, and bin-conditional), bootstrap-based intervals, and parametric methods, and is compatible with any model that produces point predictions. The package is designed to be user-friendly and extensible, allowing researchers to easily apply it to a wide range of predictive modeling tasks. The package is available from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=pintervals, and its development version, issue tracker, and contribution guidelines are hosted on GitHub at https://github.com/Doktorandahl/pintervals.
Prediction intervals quantify the uncertainty around model predictions by providing a range within which the true outcome is stated to be contained at a given confidence level, e.g., a 90% prediction interval. Unlike confidence intervals, which characterize uncertainty about fixed but unknown model parameters, prediction intervals capture the uncertainty in subsequent data realizations, making them an essential tool for evaluating and communicating the reliability of forecasting models. In applied research, particularly in the social sciences, this distinction is crucial, as stakeholders and policymakers are often more interested in the uncertainty of predicted outcomes than in parameter inference.
The pintervals package implements three classes of model-agnostic prediction intervals: conformal, bootstrapped, and parametric. These methods differ in their underlying assumptions and estimation strategies but share a common interface and workflow. Both the conformal and bootstrap approaches rely on a separate calibration set containing observed outcome values, used to calibrate prediction intervals based on the empirical distribution of prediction errors. Parametric intervals, in contrast, assume a specified distribution for the prediction errors of which the parameters can either be estimated from calibration data or supplied directly. Together, the software package provides a flexible and coherent framework for quantifying predictive uncertainty across a wide range of models and applications.
CP is a general, distribution-free method that produces prediction intervals with guaranteed marginal coverage, with minimal assumptions on the underlying prediction model. The method relies on the construction of a non-conformity score to quantify a degree of dissimilarity of a given datum point versus a collection of data examples. Under the assumption that non-conformity scores are exchangeable, their empirical distribution can be used to construct a p-value for testing the null hypothesis that the non-conformity score of a test point is exchangeable with the observed data scores, to formally test whether to rule out a candidate outcome value. The conformal p-value is a usual frequentist p-value, in that it is stochastically no smaller than a uniform(0,1) random variable, it gives finite-sample control of type 1 errors, and it corresponds to a prediction interval. Namely, a CP interval at level \(1-\alpha\) contains all outcome values for which their conformal p-value exceeds \(\alpha\). For a comprehensive background on CP, see Vovk et al. (2022), and for its connections to other statistical paradigms, see Williams (2026) and Williams and Liu (2024).
The (full/transductive) CP algorithm constructs a prediction interval for \(y_{N+1}\) given a corresponding feature \(x_{N+1}\), by comparing to a training set of feature-response pairs \((x_1,y_1), \dots, (x_N, y_N)\). The non-conformity score \(s_i\) may measure the dissimilarity between a feature-response pair \((x_i, y_i)\) and the remaining data points \((x_1, y_1), \dots, (x_{i-1}, y_{i-1}), (x_{i+1}, y_{i+1}), \ldots, (x_{N+1}, y_{N+1})\). Although the distribution of the data and of the non-conformity scores are both unknown, assuming exchangeability implies that \(s_{N+1}\) is equally likely to take any rank among the \(N+1\) scores, i.e., \(\text{rank}(s_{N+1}) \sim \text{uniform}\{1,\dots,N+1\}\). Consequently, if \(\hat{q}_{1-\alpha}\) denotes the \(1-\alpha\) empirical quantile of the \(N+1\) scores, then, assuming exchangeability, \(\mathsf{P}(s_{N+1} \leq \hat{q}_{1-\alpha}) \geq 1-\alpha\). In other words,
\[ \text{rank}(s_{N+1})/(N+1) := \frac{1}{N+1}\sum_{i=1}^{N+1} \mathbf{1}\{s_{i} \geq s_{N+1}\} \]
is a p-value of the null hypothesis that \(s_{N+1}\) is indeed exchangeable with the \(s_1, \dots, s_{N}\). A CP set, at any level \(1-\alpha \in [0,1]\), can be assembled by considering all values that \(y_{N+1}\) could take and only including those for which \(\text{rank}(s_{N+1})/(N+1) > \alpha\).
In a supervised setting, we typically create a prediction model \(\hat{f}\) such that \(\hat{f}(x_{N+1})\) represents the model’s best guess of \(y_{N+1}\). In a regression setting, the canonical non-conformity score is perhaps the absolute residual, \(s_i := \lvert y_i - \hat{f}(x_i)\rvert\), while in classification a standard choice is one minus the softmax probability assigned to the true class. In either case, the CP procedure yields a marginal prediction set
\[ C_{1-\alpha}(x_{N+1}) := \big\{ y \in \mathcal{Y} \, : \, s_{N+1}(y) \leq \hat{q}_{1-\alpha}\big\}, \]
where \(s_{N+1}(y)\) denotes the non-conformity score corresponding to \(y_{N+1} = y\). This construction ensures the following finite-sample marginal coverage guarantee,
\[ \mathsf{P}\{Y_{N+1} \in C_{1-\alpha}(X_{N+1})\} \geq 1-\alpha, \]
which holds with virtually no assumptions on the underlying prediction model or choice of non-conformity score.
There are two main variants of the CP algorithm: the transductive (or full CP) algorithm (Vovk et al. 1999; Saunders et al. 1999; Lei et al. 2013; Vovk 2015) and inductive (or split CP) algorithm (Papadopoulos et al. 2002; Papadopoulos 2008). The transductive utilizes the entire training data, making it maximally statistically efficient, but (typically) requires re-training of the prediction model for each test value \(y\), making it computationally costly. The inductive version trades off some share of statistical efficiency for a gain in computational efficiency, by setting aside a portion of the data as a calibration set, untouched during training, but requiring only a single training of the prediction model. The resulting prediction intervals constructed are typically evaluated based on empirical coverage on a test set, with reference to the nominal \(1-\alpha\) confidence level, for each \(\alpha \in [0,1]\); if the coverage gap exceeds usual Monte Carlo simulation variation, then the implication is that the non-conformity scores are not exchangeable. In the latter case, alternative formulations of the non-conformity measure and/or prediction model may better account for structural components of the data, to yield exchangeable scores. In practical applications, the size of the prediction intervals are also of interest, as “tighter” intervals are more informative to the user.
While standard CP (SCP) provides valid coverage guarantees marginally, it does not account for potential heterogeneity in the data, which can lead to over- or under-coverage in certain regions of the data. To address this, Mondrian CP (MCP) extends the conformal framework by partitioning the calibration set into groups based on covariates or the outcome label (Vovk et al. 2003). The disjoint groups enable the construction of group-conditional prediction intervals, i.e., resulting in a CP set for each group. Supposing the feature space \(\mathcal{X}\) can be partitioned into disjoint subsets \(\mathcal{X}_1, \dots, \mathcal{X}_K\), MCP guarantees, for each \(k \in \{1, \dots, K\}\) and \(\alpha \in [0,1]\), that
\[ \mathsf{P}\{Y_{N+1} \in C_{1-\alpha}(x_{N+1}) \mid x_{N+1} \in \mathcal{X}_k\} \geq 1-\alpha. \]
In MCP, the calibration step is performed separately per group rather than across the entire calibration data, involving the computation (for each \(\alpha\)) of \(K\), \(1-\alpha\) level quantiles \(\hat{q}_{1-\alpha}^{(1)}, \dots, \hat{q}_{1-\alpha}^{(K)}\). While group-conditional validity may hold in MCP, the limiting factor for efficiency (e.g., narrow intervals) is that a sufficient number of calibration points must be available in each group. Finally, it is known (Lei and Wasserman 2014) that exact conditional coverage guarantees, i.e., achieving valid coverage conditional on the exact value \(x_{N+1}\) rather than on a Mondrian group membership, are impossible in the context of CP.
CCP is an extension of MCP for situations where the number of groups \(K\) is large and/or some groups contain too few calibration points to efficiently estimate the required quantile. The approach, first proposed in Ding et al. (2023), is to cluster the \(K\) initial groups into \(M \ll K\) clusters, and perform calibration at the cluster level rather than at the group level. By pooling information across similar groups, CCP increases the statistical efficiency and may also help to avoid numerical issues arising from too few examples per group.
The key idea of CCP is to cluster together groups whose empirical distribution of non-conformity scores are similar. As shown in Ding et al. (2023), the expected coverage gap for any given group depends directly on the similarity between the groups clustered together (note that CCP guarantees validity at the cluster level). When the empirical cumulative distribution functions clustered together are similar, the coverage gap is small. In practice, clustering can be done, for example, through k-means-clustering or based on the Kolmogorov-Smirnov distance between empirical cumulative distribution functions (Ding et al. 2023; Hjort et al. 2024).
It is known that the SCP methods often under-cover in lower density regions of the outcome space \(\mathcal{Y}\), while over-covering in higher density regions of \(\mathcal{Y}\) (Guan 2023). While MCP and CCP are useful for addressing under- and over-coverage for different groups in the data, these methods are not easily applicable if the outcome space has no natural or well-defined grouping structure. One example is data arising from skewed distributions. BCCP is a special case of the MCP that addresses this issue by partitioning the calibration data into user-defined bins of the outcome space, such as low, medium, and high values. BCCP can thus be conceptualized as an MCP method for continuous outcomes, where the bins are analogous to the labels in label-conditional MCP. The BCCP method is described in detail in Randahl et al. (2025).
For BCCP, non-conformity scores are calculated on a calibration data set within each user-defined bin. As the true bin is not known for new test points, for a new test point, the BCCP algorithm calculates its bin-specific non-conformity score within each bin. A CP set is then constructed to include all bins for which the bin-specific conformal p-value of the test point exceeds \(\alpha\). For univariate outcomes, the final prediction set is formed by either the union of disjoint intervals corresponding to each bin in the CP set, or by contiguizing the intervals using the left- and right-most endpoints across the bins in the CP set. Both formations achieve at least the desired coverage level, but the contiguized intervals may suffer over-coverage.
BCCP is particularly useful in cases where the outcome space is skewed or imbalanced, as it allows for more flexible and adaptive CP sets that can better reflect the distribution of the outcome variable versus SCP. It also allows the user to define the bins based on their knowledge of the data and the specific application, providing a more tailored approach to uncertainty quantification. Similar to MCP and CCP, BCCP requires carefully balancing the number and specification of bins with the quantity of data within each bin, to maximize statistical efficiency.
Another modification of the CP procedure is a class of methods referred to as weighted CP. Rather than calculating the quantile of interest once (or once per group or cluster), weighted CP calculates \(\hat{q}_{1-\alpha}\) as a weighted quantile, where the scores \(s_{1}, \dots, s_{N}\) are weighted according to some weights \(w_{1}, \dots, w_{N}\). By calculating a weighted quantile, the method places more emphasis on non-conformity scores that are determined closer to exchangeable with the test instance, e.g., with reference to time or space. For example, a weighted CP framework could take weights of the form \(w_i \propto \exp\{-\text{dist}(x_i, x_{N+1})\}\), where \(\text{dist}(x_i, x_{N+1})\) is an application-specific distance function quantifying the distance between calibration and test features. The prediction interval for \(y_{N+1}\) is then calibrated using the weighted \(1-\alpha\) quantile of the scores.
The literature contains several variations of this, including data-driven weights to account for covariate shift (Tibshirani et al. 2019), fixed weights to handle non-exchangeable data (Foygel Barber et al. 2023), spatial weights (Mao et al. 2023), and feature-distance weights (Guan 2023). The coverage guarantees of the method depend on the choice of weighting function; see Barber and Tibshirani (2025) for a review of the statistical properties of weighted CP methods.
The performance of the distance-weighted CP (DWCP) methods relies heavily on the choice of distance metric and kernel. The weighting methods might suffer from the “curse of dimensionality” in high-dimensional data sets, in particular if the features are not normalized to a common scale. Choosing an overly narrow kernel reduces the effective sample size used for the estimate of the calibration quantile (Foygel Barber et al. 2023), which can lead to arbitrarily poor local coverage (Guan 2023). To mitigate this, Guan (2023) propose a lower-dimensional weighting function (referred to as a localizer) tailored to the CP application, but demonstrate its usefulness only for low-dimensional data. Furthermore, the choice of distance metric should be tailored to the specific type of drift present in the data, as Euclidean distance or Mahalanobis distance is not necessarily the right choice for every application. Lee et al. (2025), for instance, develop a weight function for weighted CP under temporal drift, while the choice of weighting method in spatial settings is discussed in Lou et al. (2025) and Mao et al. (2023). Given this sensitivity to the choice of distance metric, kernel, and bandwidth, we recommend that practitioners report results under multiple weighting specifications as a sensitivity check.
Outside the realm of CP, another widely used approach to constructing prediction intervals is based on bootstrapping. The bootstrap is a resampling technique that approximates the sampling distribution of any quantity of interest (e.g., a point prediction) by repeatedly drawing samples with replacement from the observed data. This approach is particularly appealing when the distribution of prediction errors is unknown or complex, as it provides an empirical estimate of predictive uncertainty without requiring parametric assumptions (Beran 1990; Tian et al. 2022).
Bootstrapped prediction intervals are expected to achieve nominal coverage under the assumption that the calibration set is representative of the data-generating process, and that the prediction errors are approximately independent and identically distributed. The algorithm proceeds by first computing the prediction errors on the calibration set and then repeatedly resampling these errors with replacement to form an empirical error distribution (Tian et al. 2022; Hesterberg 2015). For each new test observation, a large number of resampled errors are added to the model’s point prediction, generating a simulated distribution of possible outcomes. The lower and upper bounds of the prediction interval are then obtained by taking \(\hat{q}_{\alpha/2}\) and \(\hat{q}_{1-\alpha/2}\).
Several variants of the bootstrap exist, including the residual bootstrap, wild bootstrap, and parametric bootstrap. Variants differ both in how the resampling is performed, and in the assumptions they require about the underlying model or error structure.
Parametric prediction intervals are constructed under the assumption that the prediction errors follow a specific probability distribution. This approach models uncertainty analytically rather than empirically, allowing the interval bounds to be derived directly from the assumed distribution and its estimated parameters. Common choices include the normal, log-normal, Poisson, or negative binomial distributions, depending on the nature of the outcome variable and the model’s residuals.
For example, if the prediction errors are assumed to follow a centered normal distribution, their standard deviation can be estimated from a calibration set, and the prediction interval for a new observation is obtained by adding and subtracting the appropriate quantiles of the normal distribution around the point prediction. More generally, any distribution for which quantile functions are available can be used to construct prediction intervals, including user-specified or mixture distributions.
The main advantage of parametric prediction intervals lies in their computational efficiency and smooth, easily interpretable behavior. However, their validity depends critically on the appropriateness of the assumed error distribution: when this assumption is violated the resulting intervals may exhibit biased or unreliable coverage.
The methods described above differ both in the assumptions they require and in the type of coverage they guarantee. The table below summarizes these differences. The CP-based methods provide finite-sample guarantees under exchangeability, either marginally (SCP) or conditionally on feature-defined groups (MCP), clusters of groups (CCP), or bins of the outcome space (BCCP). For DWCP, the guarantee depends on how well the weights capture departures from exchangeability. Bootstrapped and parametric intervals, in contrast, offer no general finite-sample guarantee. Their validity rests on asymptotic arguments or on correct specification of the error distribution. In practice, the choice between methods therefore involves a trade-off between the strength of the coverage guarantee, the level at which coverage is required (marginal, group, or outcome range), and the amount of calibration data available within each group or bin.
| Method | Main assumption | Coverage guarantee |
|---|---|---|
| SCP | Exchangeability of \((X_i, Y_i)_{i=1}^{N+1}\) | Marginal: \(\mathsf{P}\{Y_{N+1} \in C(X_{N+1})\} \geq 1-\alpha\) |
| MCP | Exchangeability within each subset \(\mathcal{X}_k\) of the feature space | Group-conditional: \(\mathsf{P}\{Y_{N+1} \in C(X_{N+1}) \mid X_{N+1} \in \mathcal{X}_k\} \geq 1-\alpha\) for all \(k\) |
| CCP | Exchangeability within each cluster \(\mathcal{C}_m\) of subsets with similar score distributions, with clusters formed on separate data | Cluster-conditional: \(1-\alpha\) within each \(\mathcal{C}_m\); for each subset \(\mathcal{X}_k \subseteq \mathcal{C}_m\), \(\mathsf{P}\{Y_{N+1} \in C(X_{N+1}) \mid X_{N+1} \in \mathcal{X}_k\} \geq 1-\alpha-\varepsilon_m\), where \(\varepsilon_m\) depends on the dissimilarity of the subsets pooled in \(\mathcal{C}_m\) |
| BCCP | Exchangeability within pre-specified bins \(\mathcal{Y}_b\) of the outcome space | Bin-conditional: \(\mathsf{P}\{Y_{N+1} \in C(X_{N+1}) \mid Y_{N+1} \in \mathcal{Y}_b\} \geq 1-\alpha\) for all \(b\) |
| DWCP | Exchangeability holds after reweighting (weights correctly capture the exchangeability violation) | Depends on the weights: marginal if well-specified, otherwise approximate or local |
| Bootstrap | Calibration errors are i.i.d. and representative of the test errors (error distribution does not depend on \(X\)) | None in finite samples; approximately \(1-\alpha\) marginal coverage as the calibration set grows, with a tendency to under-cover in small samples |
| Parametric | Errors follow the specified distribution \(F_\theta\), with \(\theta\) known or consistently estimated | Exact marginal coverage if \(F_\theta\) is correct and \(\theta\) known; asymptotic if \(\theta\) is estimated; none under misspecification |
The pintervals package provides a unified and model-agnostic interface for constructing prediction intervals using the methods described above. All functions in the package share a consistent syntax and workflow, allowing users to easily switch between different approaches and apply them to a wide range of predictive models. Each function takes as input a vector of predicted values for the test set, as well as predicted and true outcomes from a calibration set2, either as a two-column tibble or as separate vectors, and returns a tibble with the lower and upper bounds of the prediction intervals, alongside the corresponding point predictions.
All CP-based methods are implemented using the inductive/split CP algorithm. The transductive/full CP algorithm was only discussed conceptually above. For methods that rely on partitioning the data, such as MCP, CCP, or BCCP, an additional grouping variable must be specified. This can be provided either as a separate vector or as a column in the calibration tibble.
Prediction intervals are constructed using the
pinterval_*() function family, which includes the following
main functions:
pinterval_conformal() for standard conformal prediction
intervals;pinterval_mondrian() for Mondrian conformal prediction
intervals;pinterval_ccp() for clustered conformal prediction
intervals;pinterval_bccp() for bin-conditional conformal
prediction intervals;pinterval_bootstrap() for bootstrapped prediction
intervals;pinterval_parametric() for parametric prediction
intervals.All functions follow the same general structure and support additional arguments to control key aspects of interval construction—such as the desired coverage level, the number of bootstrap replications, or the distributional parameters for parametric intervals. This unified design facilitates flexible experimentation and direct comparison across different interval estimation methods within a common analytical framework. The output of each function is a tibble containing the point predictions along with the lower and upper bounds of the prediction intervals and any additional relevant information (e.g., group or cluster assignments for MCP and CCP).
The pinterval_conformal() function implements SCP
intervals using the inductive CP algorithm. It takes the predicted
values from the test set and the calibration data as input, and returns
a tibble with the lower and upper bounds of the prediction intervals.
The calibration data should be a two-column tibble of predicted (first
column) and true (second column) values from the calibration set;
alternatively the predicted values can be passed as a vector for the
argument calib and the true values as a vector for the
argument calib_truth.
As additional arguments, the user can specify the desired coverage level and the non-conformity score function to be used. Absolute errors are used by default, with relative, zero-adjusted relative, heterogeneous, and raw errors as the alternatives. The user can also specify a custom real-valued function that takes the predicted and true values as input and returns a vector of non-conformity scores. This allows for flexibility in defining the non-conformity score, based on the specific application or data characteristics.
The functions pinterval_mondrian() and
pinterval_ccp() implement MCP and CCP, which allows for
partitioning the calibration set into disjoint groups based on a
grouping variable.3
The syntax for MCP and CCP builds on
pinterval_conformal(), but requires an additional column in
the calibration tibble for the grouping variable (third column). The
grouping variable can also be passed as a vector for the argument
calib_class.
For CCP, the user can also define the number of clusters to be used, or choose to optimize the number of clusters by minimizing the Caliński-Harabasz index of the proposed clustering structure (Caliński and Harabasz 1974). To avoid overfitting of the clusters, the calibration data can be further divided into one partition which is used to cluster the classes, and one partition used to calculate the resulting prediction intervals.4
The function pinterval_bccp() implements BCCP intervals,
which allows for partitioning the calibration set into user-defined bins
of the outcome space. The syntax is similar to
pinterval_mondrian() and pinterval_ccp(), in
that it uses an additional column, bins, in the calibration
tibble to specify the bins, or a vector for the argument
calib_bins. The user can also specify whether to contiguize
the intervals using the argument contiguize, which is set
to FALSE by default. If set to TRUE, the
function will return contiguous intervals by taking the left- and
right-most endpoints across the bins, otherwise it will contiguize
adjacent bin-specific intervals. When contiguize is
FALSE the function outputs a tibble where the prediction
intervals with corresponding lower- and upper bounds are in a
list-column called intervals, otherwise the function
outputs a standard tibble with point predictions and the contiguized
lower- and upper bounds of the intervals. The user can also specify the
desired coverage level and the non-conformity score function to be used,
as in pinterval_conformal().
All CP methods in pintervals also support DWCP, which gives more weight to non-conformity scores from observations in the calibration set with values close to the new test observation.
DWCP is enabled by setting distance_weighted_cp = TRUE.
DWCP additionally requires calibration and test covariates, supplied
through the arguments distance_features_calib and
distance_features_pred. Distances are calculated as either
the Mahalanobis distance (default) or the Euclidean distance between the
rows of the calibration and test covariate data. Optionally, the
distances can be normalized after calculation (preferred for Euclidean
distances), either through min-max normalization or by dividing by the
standard deviation of the distances (z-score normalization). When DWCP
is used, the non-conformity scores from the calibration set are weighted
according to their distance to the new test observation before computing
the quantiles for the prediction interval.
In addition, the user can specify a weight function, which translates
the distances into weights for the non-conformity scores. By default,
the function uses the Gaussian kernel: \(w =
\exp(-d^2)\), where \(d\) is the
Euclidean distance between the predicted value of the new test
observation and the predicted values in the calibration set. Built-in
alternatives are the Cauchy kernel: \(w =
\frac{1}{1+d^2}\), the logistic kernel: \(w = \frac{1}{1+\exp(d)}\), and the
reciprocal linear kernel: \(w =
\frac{1}{1+d}\). The user can also specify a custom
distance-weighting function using the argument
distance_weighting_function. This function should take a
vector of distances as input and return a vector of non-negative
weights, which will be used to weight the non-conformity scores. This
allows for considerable flexibility in defining the distance-weighting
function based on the specific application or data characteristics.
The pinterval_bootstrap() function implements
bootstrapped prediction intervals. It uses the prediction errors from
the calibration set to mimic the distribution of errors via resampling.
The function takes as input the predicted values for the test set and
the calibration data, either as a two-column tibble (predicted values in
the first column and true values in the second), or as separate vectors
supplied to the arguments calib and
calib_truth.
As additional arguments, the user can specify the desired coverage
level (default is 0.9), the number of bootstrap samples (default is
1000), and the type of error to bootstrap from. The argument
error_type can be set to "raw" to use raw,
signed, errors, or "absolute" to use absolute errors with
random signs. The user can also specify whether to use distance-weighted
resampling by setting the argument distance_weighted to
TRUE. This will give more weight to errors from
observations in the calibration set with predicted values close to the
new test observation, which can help improve the accuracy of the
prediction intervals. The distance-weighted bootstrap procedure follows
the same principles as DWCP described above, being controlled by the
arguments distance_features_calib,
distance_features_pred, distance_type,
normalize_distance, and weight_function.
The pinterval_parametric() function implements
parametric prediction intervals based on a user-specified probability
distribution. The user supplies the predicted values for the test set, a
desired distribution (e.g., normal, log-normal, Poisson, negative
binomial), and either a calibration set with predicted and true values
or the parameters of the distribution.
Natively supported distributions include common continuous and count
distributions such as the normal, log-normal, Poisson, and negative
binomial distributions. The user specifies the distribution using the
argument dist. For native distributions, the relevant
parameters of the distribution (e.g., standard deviation from the normal
distribution and dispersion parameter from the negative binomial
distribution) are estimated from the calibration set. Alternatively, the
user can provide the parameters directly by passing a named list with
the parameters through the argument pars. Unlike the other
pinterval_*() functions,
pinterval_parametric() does not require calibration data if
the parameters of the distribution are provided directly. Similarly, for
one-parameter distributions (e.g., the Poisson and Chi-square
distributions), no calibration data are needed as the single parameter
can be estimated directly from the predicted values.
Beyond the natively supported distributions, the user can also use
any arbitrary distribution by supplying a custom quantile function to
the argument dist. This allows for considerable
flexibility, including custom or mixture distributions using for
instance the mistr package (Sablica and Hornik 2023, 2020). When specifying
a custom distribution, the user needs to provide the parameters of the
distribution as a named list to the argument pars. Ideally,
these parameters should be estimated on a calibration data set. The
function will then use the specified quantile function to compute the
lower and upper bounds of the prediction intervals.
Parametric intervals are simple, efficient and interpretable, but rely on correct specification of the underlying distribution. Users are encouraged to verify distributional assumptions using the calibration errors before applying this method.
To illustrate the practical use of the pintervals package, we demonstrate how it can be applied to generate prediction intervals for county-level relative vote share for the Democratic party5 in the 2024 US presidential election. This example provides a concrete, real-world application of the methods introduced above, using a typical regression problem workflow where predictions are obtained as point-predictions from a fitted model and prediction intervals are constructed using the methods available in pintervals to quantify the uncertainty around these predictions.
The data are bundled with the pintervals package as
the county_turnout_panel dataset. It is a county-level
panel covering the five US presidential elections of 2008, 2012, 2016,
2020, and 2024, with one row per county per election year for the 50
states and the District of Columbia. The county presidential vote
returns are taken from the MIT Election Data and Science Lab (MIT Election Data and Science Lab 2025), while
the demographic and socioeconomic covariates are drawn directly from the
US Census Bureau API, with each covariate matched to the American
Community Survey five-year vintage ending nearest the election (U.S. Census Bureau 2024). The dataset includes
variables such as total population, median household income, educational
attainment, and the racial, age, and gender composition of the counties,
along with each county’s US Census region and division and its centroid
coordinates. The outcome variable of interest in this example is the
county-level relative Democratic vote share, defined as the Democratic
party’s share of the presidential vote in the county minus the party’s
vote share in the electorate as a whole in the same election. Unlike the
cross-sectional county_turnout dataset, the panel does not
contain precomputed point predictions; instead, we generate
out-of-sample predictions within the example below by training a random
forest (Wright and Ziegler 2017)
regression model on the earlier elections, using eight county-level
demographic and socioeconomic covariates as predictors.
The goal of the analysis is to construct prediction intervals around the point predictions of the relative vote share of the Democratic party using the different methods implemented in pintervals and evaluate their empirical coverage. We will compare the performance of the different CP methods (standard, Mondrian, clustered, distance-weighted, and bin-conditional), as well as bootstrapped prediction intervals and parametric prediction intervals assuming normal distributions of errors.
We begin the example below by loading the necessary packages, the county turnout panel dataset, and splitting the data into training, calibration and test sets.
As a training dataset, we use all county-years for the presidential elections in 2008–2016, for a total of 9,325 county-year observations. We then use the 3,105 counties for the presidential election in 2020 as calibration data, and the 3,097 counties in 2024 as the test data. This mimics a real-world forecasting exercise, where the training, calibration, and test data sets are sequential.
To generate the point predictions, we then train a random forest model with the relative vote share for the democratic party as the prediction target, and a set of county-level characteristics as the predictors. We use this random forest model to generate out-of-sample predictions for both the calibration and test sets, without retraining the model in between.
library(pintervals)
library(dplyr)
library(tidyr)
library(ranger)
data("county_turnout_panel", package = "pintervals")
set.seed(101010)
train <- county_turnout_panel %>%
filter(election_year < 2020) %>%
mutate(
rel_dem_share = dem_share - dem_share_electorate
) %>%
na.omit()
calib <- county_turnout_panel %>%
filter(election_year == 2020) %>%
mutate(
rel_dem_share = dem_share - dem_share_electorate
) %>%
na.omit()
test <- county_turnout_panel %>%
filter(election_year == 2024) %>%
mutate(
rel_dem_share = dem_share - dem_share_electorate
) %>%
na.omit()
m <- ranger(
rel_dem_share ~ nonwhite_pct + foreignborn_pct + female_pct +
age29andunder_pct + median_hh_inc + clf_unemploy_pct +
lesscollege_pct + rural_pct,
data = train
)
calib$predicted <- predict(m, data = calib)$predictions
test$predicted <- predict(m, data = test)$predictionsAfter obtaining the point predictions for the relative vote share of
the democratic party, we can proceed to construct prediction intervals
using the methods available in pintervals. For each
method, we will use the point predictions from the random forest model
and the true vote share for the Democratic party from the calibration
set to generate the intervals for the test set. We will then evaluate
the empirical coverage using the interval_coverage()
function from the package. We begin with SCP intervals, constructed
using pinterval_conformal(), followed by an evaluation of
their empirical coverage. In Appendix A, we
repeat this procedure evaluating the predictions using the interval
score instead of the coverage.
conformal_intervals <- pinterval_conformal(
pred = test$predicted,
calib = calib$predicted,
calib_truth = calib$rel_dem_share,
alpha = 0.1
)
conformal_coverage <- interval_coverage(
truth = test$rel_dem_share,
lower_bound = conformal_intervals$lower_bound,
upper_bound = conformal_intervals$upper_bound
)
conformal_coverage
#> [1] 0.8898934We can repeat this procedure for the parametric and bootstrapped prediction intervals, adjusting the function calls as needed to accommodate the specific requirements of each method. For the parametric intervals, we assume a normal distribution of the prediction errors.
norm_intervals <- pinterval_parametric(
pred = test$predicted,
calib = calib$predicted,
calib_truth = calib$rel_dem_share,
dist = "norm",
alpha = 0.1
)
norm_coverage <- interval_coverage(
truth = test$rel_dem_share,
lower_bound = norm_intervals$lower_bound,
upper_bound = norm_intervals$upper_bound
)
bootstrap_intervals <- pinterval_bootstrap(
pred = test$predicted,
calib = calib$predicted,
calib_truth = calib$rel_dem_share,
alpha = 0.1,
n_bootstrap = 1000
)
bootstrap_coverage <- interval_coverage(
truth = test$rel_dem_share,
lower_bound = bootstrap_intervals$lower_bound,
upper_bound = bootstrap_intervals$upper_bound
)
c(norm_coverage, bootstrap_coverage)
#> [1] 0.8918308 0.8844043Here we can see that all methods achieve coverage close to the nominal 90% level. In Appendix B, we present a simulation study comparing the performance of these methods across multiple calibration-test splits.
Next, we demonstrate how to use MCP and CCP to improve coverage within subsets of the data. In this example, we will use the US Census region for the counties as the grouping variable for the MCP. For the CCP, we use the US Census divisions as the grouping variable. In addition, we utilize DWCP using geographic (Euclidean) distances between counties.6 We will then evaluate the empirical coverage both in aggregate and within each of the four US Census regions.
For the MCP method the user must ensure that the calibration data are sufficiently large in each group to obtain reliable estimates of the non-conformity score quantiles. If a group has too few calibration points, the resulting prediction intervals may be too wide or too narrow, leading to poor coverage. In this example, the calibration data contains 1,420 counties in the US Census Region South, 1,054 counties in the Midwest, 414 counties in the West, and 217 counties in the Northeast.
For the CCP method, however, the important consideration is not the number of calibration points in each pre-defined group, which in our example range from 67 in the New England Division to 617 in the West North Central Division, but rather the number of calibration points in each proposed cluster. The number of clusters can be optimized using either the Caliński-Harabasz index (default) or a desired minimum cluster size. In this example, we use the Caliński-Harabasz index to optimize the number of clusters and set the maximum number of clusters to 5.
mondrian_intervals <- pinterval_mondrian(
pred = test$predicted,
pred_class = test$region,
calib = calib$predicted,
calib_truth = calib$rel_dem_share,
calib_class = calib$region,
alpha = 0.1
)
mondrian_coverage <- interval_coverage(
truth = test$rel_dem_share,
lower_bound = mondrian_intervals$lower_bound,
upper_bound = mondrian_intervals$upper_bound
)
clustered_conformal_intervals <- pinterval_ccp(
pred = test$predicted,
pred_class = test$division,
calib = calib$predicted,
calib_truth = calib$rel_dem_share,
calib_class = calib$division,
alpha = 0.1,
optimize_n_clusters = TRUE,
max_n_clusters = 5
)
clustered_conformal_coverage <- interval_coverage(
truth = test$rel_dem_share,
lower_bound = clustered_conformal_intervals$lower_bound,
upper_bound = clustered_conformal_intervals$upper_bound
)
distance_weighted_cp_intervals <- pinterval_conformal(
pred = test$predicted,
calib = calib$predicted,
calib_truth = calib$rel_dem_share,
alpha = 0.1,
distance_weighted_cp = TRUE,
distance_features_calib = calib %>%
select(latitude, longitude),
distance_features_pred = test %>%
select(latitude, longitude),
distance_type = "euclidean",
normalize_distance = "sd"
)
distance_weighted_cp_coverage <- interval_coverage(
truth = test$rel_dem_share,
lower_bound = distance_weighted_cp_intervals$lower_bound,
upper_bound = distance_weighted_cp_intervals$upper_bound
)
c(mondrian_coverage, clustered_conformal_coverage,
distance_weighted_cp_coverage)
#> [1] 0.8931224 0.9015176 0.8911850We again see that the methods achieve coverage close to the nominal 90% level in aggregate. To evaluate the group-wise coverage within each US Census region, we first add region, method, and the relative Democratic party vote share to each set of intervals, combine the results from all CP methods, and then compute the empirical coverage within each group.
conformal_intervals <- conformal_intervals %>%
mutate(
region = test$region,
method = "SCP",
rel_dem_share = test$rel_dem_share
)
mondrian_intervals <- mondrian_intervals %>%
mutate(
region = test$region,
method = "MCP",
rel_dem_share = test$rel_dem_share
)
clustered_conformal_intervals <- clustered_conformal_intervals %>%
mutate(
region = test$region,
method = "CCP",
rel_dem_share = test$rel_dem_share
)
distance_weighted_cp_intervals <- distance_weighted_cp_intervals %>%
mutate(
region = test$region,
method = "DWCP",
rel_dem_share = test$rel_dem_share
)
bootstrap_intervals <- bootstrap_intervals %>%
mutate(
region = test$region,
method = "Bootstrap",
rel_dem_share = test$rel_dem_share
)
norm_intervals <- norm_intervals %>%
mutate(
region = test$region,
method = "Normal",
rel_dem_share = test$rel_dem_share
)
all_conformal_intervals <- bind_rows(
conformal_intervals,
mondrian_intervals,
clustered_conformal_intervals,
distance_weighted_cp_intervals,
bootstrap_intervals,
norm_intervals
)
group_wise_coverage <- all_conformal_intervals %>%
group_by(method, region) %>%
summarise(
coverage = interval_coverage(
truth = rel_dem_share,
lower_bound = lower_bound,
upper_bound = upper_bound
),
.groups = "drop"
)
group_wise_coverage %>%
pivot_wider(
names_from = region,
values_from = coverage
)
#> # A tibble: 6 × 5
#> method Midwest Northeast South West
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 Bootstrap 0.924 0.804 0.874 0.860
#> 2 CCP 0.910 0.852 0.911 0.872
#> 3 DWCP 0.910 0.785 0.898 0.875
#> 4 MCP 0.898 0.856 0.903 0.867
#> 5 Normal 0.919 0.789 0.897 0.855
#> 6 SCP 0.918 0.789 0.894 0.855The results show that the Mondrian and CCP algorithms achieve coverage that is more uniform and close to the nominal level across all US Census regions compared to SCP, bootstrapped, and parametric prediction intervals.
Analyzing the results further, we below compute the mean absolute error (MAE) of coverage across regions for each method to quantify the uniformity of coverage. The results show that the Mondrian and clustered CP have MAE values around half of those of the alternative methods. The simulation study in Appendix B further supports these findings, showing that the group-aware CP methods consistently achieve coverage closer to the nominal level across the US Census regions.
mae_coverage_regions <- group_wise_coverage %>%
group_by(method) %>%
summarize(mae_coverage = mean(abs(coverage - 0.9))) %>%
pivot_wider(
names_from = method,
values_from = mae_coverage
)
mae_coverage_regions
#> # A tibble: 1 × 6
#> Bootstrap CCP DWCP MCP Normal SCP
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 0.0466 0.0243 0.0382 0.0203 0.0444 0.0446Finally, we illustrate how BCCP can be used to improve coverage within different ranges of the outcome variable. In this example, we will create bins for relative Democratic party vote shares lower than -30%, between -30% and -20%, between -20% and -10%, and above -10% and generate prediction intervals using BCCP.7 We will then evaluate the empirical coverage both in aggregate and within each bin.
calib <- calib %>%
mutate(
truth_bin = case_when(
rel_dem_share < -0.3 ~ 1,
rel_dem_share < -0.2 ~ 2,
rel_dem_share < -0.1 ~ 3,
TRUE ~ 4
)
)
bccp_contiguized_intervals <- pinterval_bccp(
pred = test$predicted,
calib = calib$predicted,
calib_truth = calib$rel_dem_share,
calib_bins = calib$truth_bin,
alpha = 0.1,
contiguize = TRUE
)
bccp_contiguized_coverage <- interval_coverage(
truth = test$rel_dem_share,
lower_bound = bccp_contiguized_intervals$lower_bound,
upper_bound = bccp_contiguized_intervals$upper_bound
)
bccp_discontigous_intervals <- pinterval_bccp(
pred = test$predicted,
calib = calib$predicted,
calib_truth = calib$rel_dem_share,
calib_bins = calib$truth_bin,
alpha = 0.1,
contiguize = FALSE
)
bccp_discontigous_coverage <- interval_coverage(
truth = test$rel_dem_share,
intervals = bccp_discontigous_intervals$intervals
)
c(bccp_contiguized_coverage, bccp_discontigous_coverage)
#> [1] 0.9053923 0.8898934Both contiguized and non-contiguized BCCP intervals achieve coverage close to the nominal 90% level in aggregate. To evaluate the bin-wise coverage within each bin of relative Democratic vote share, we add the bin, method, and vote share to each set of test-set intervals, combine the results from both BCCP variants and the other methods, and then compute the empirical coverage within each bin.
test <- test %>%
mutate(
truth_bin = case_when(
rel_dem_share < -0.3 ~ 1,
rel_dem_share < -0.2 ~ 2,
rel_dem_share < -0.1 ~ 3,
TRUE ~ 4
)
)
bccp_contiguized_intervals <- bccp_contiguized_intervals %>%
mutate(
truth_bin = test$truth_bin,
method = "BCCP(c)",
rel_dem_share = test$rel_dem_share
)
bccp_discontigous_intervals <- bccp_discontigous_intervals %>%
mutate(
truth_bin = test$truth_bin,
method = "BCCP(d)",
rel_dem_share = test$rel_dem_share
)
conformal_intervals <- conformal_intervals %>%
mutate(
truth_bin = test$truth_bin,
method = "SCP",
rel_dem_share = test$rel_dem_share
)
bootstrap_intervals <- bootstrap_intervals %>%
mutate(
truth_bin = test$truth_bin,
method = "Bootstrap",
rel_dem_share = test$rel_dem_share
)
norm_intervals <- norm_intervals %>%
mutate(
truth_bin = test$truth_bin,
method = "Normal",
rel_dem_share = test$rel_dem_share
)
all_intervals_bins <- bind_rows(
bccp_contiguized_intervals,
bccp_discontigous_intervals,
conformal_intervals,
bootstrap_intervals,
norm_intervals
)
bin_wise_coverage <- all_intervals_bins %>%
group_by(method, truth_bin) %>%
summarise(
coverage = interval_coverage(
truth = rel_dem_share,
lower_bound = lower_bound,
upper_bound = upper_bound,
intervals = intervals
),
.groups = "drop"
)
bin_wise_coverage_wide <- bin_wise_coverage %>%
pivot_wider(
names_from = truth_bin,
values_from = coverage
)
bin_wise_coverage_wide
#> # A tibble: 5 × 5
#> method `1` `2` `3` `4`
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 BCCP(c) 0.929 0.966 0.885 0.842
#> 2 BCCP(d) 0.928 0.929 0.865 0.842
#> 3 Bootstrap 0.761 0.975 0.958 0.826
#> 4 Normal 0.827 0.980 0.955 0.801
#> 5 SCP 0.819 0.980 0.955 0.800As shown above, both contiguized and non-contiguized BCCP achieve coverage closer to the nominal level within each bin compared to SCP, bootstrapped, and parametric prediction intervals. Bins are defined as: (1) relative vote share < -0.3, (2) -0.3 \(\leq\) relative vote share < -0.2, (3) -0.2 \(\leq\) relative vote share < -0.1, and (4) relative vote share \(\geq\) -0.1.
Analyzing the results further, we below compute the MAE of coverage across bins for each method to quantify their performance. The MAE is lowest for the non-contiguized BCCP method, followed closely by the contiguized version, while all other methods show substantially higher MAE. The simulation study in Appendix B further supports these findings, showing that BCCP consistently achieves coverage closer to the nominal level within different ranges of the outcome variable.
mae_coverage_bins <- bin_wise_coverage %>%
group_by(method) %>%
summarize(mae_coverage = mean(abs(coverage - 0.9))) %>%
pivot_wider(
names_from = method,
values_from = mae_coverage
)
mae_coverage_bins
#> # A tibble: 1 × 5
#> `BCCP(c)` `BCCP(d)` Bootstrap Normal SCP
#> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 0.0423 0.0376 0.0865 0.0768 0.0790This paper has introduced the pintervals package for R, a unified and extensible framework for constructing prediction intervals using conformal, bootstrapped, and parametric methods. The package is designed to lower the barrier for applied researchers and practitioners to incorporate formal uncertainty quantification into their predictive analyses. By offering a consistent syntax and flexible interface, pintervals allows users to move seamlessly between different approaches to interval construction while maintaining transparent and reproducible workflows.
Through an empirical example predicting the county-level relative vote share of the Democratic party in the 2024 US presidential election, we illustrated how the package can be applied to real-world data, and how alternative methods differ in coverage performance. The example further demonstrated how extensions such as MCP, CCP, BCCP, and DWCP can address heterogeneity and local calibration challenges across different subgroup dimensions of the data or prediction target. The results confirm that conformal methods achieve empirical coverage close to the nominal level not only overall but also within meaningful subpopulations, an important property in applied research where distributional assumptions may fail or subgroup balance is imperfect. These results are further validated in a simulation study presented in Appendix B, showing the robustness and reliability of the CP methods implemented in pintervals.
Beyond the example presented here, pintervals can be readily integrated into a wide range of predictive modeling workflows, from classical regression to complex machine learning pipelines. The design emphasizes interpretability and extensibility, making it straightforward to evaluate new methods or tailor existing ones to domain-specific needs. The package can be applied to any setting where exchangeability holds, at least approximately or within groups, including temporal and spatial data, as illustrated by the sequential train-calibration-test design and the geographically weighted DWCP intervals in the example above. Future work could extend the package with additional functionality, such as adaptive updating of intervals as new data become available, or joint conditioning on feature-defined groups (as in MCP and CCP) and outcome bins (as in BCCP).8
We believe that pintervals will provide an accessible, principled, and reproducible way for researchers and analysts to quantify predictive uncertainty. By combining simplicity of use with rigorous statistical foundations, the package aims to bridge the gap between methodological innovation and practical implementation, enabling more-trustworthy inference and decision-making in applied predictive modeling.
The results in this vignette were obtained using R 4.5.2 (R Core Team 2025) with the pintervals 1.3.3, ranger 0.17.0 (Wright and Ziegler 2017; Wright 2025), dplyr 1.2.1 (Wickham et al. 2026), and tidyr 1.3.2 (Wickham et al. 2025) packages. R itself and all packages used are available from CRAN at https://CRAN.R-project.org/. The development version of pintervals is available at https://github.com/Doktorandahl/pintervals.
This appendix repeats the worked example of the example section, using the interval score as the metric for comparison rather than the empirical coverage. For a central \(1 - \alpha\) prediction interval \([l, u]\) and an observed outcome \(y\) the interval score is defined as:
\[ S_\alpha(l, u; y) = (u - l) + \frac{2}{\alpha}(l - y)\,\mathbf{1}\{y < l\} + \frac{2}{\alpha}(y - u)\,\mathbf{1}\{y > u\} \]
The interval score is a proper scoring rule that rewards narrow
intervals while penalizing intervals that fail to contain the outcome.
Lower values indicate sharper, better-calibrated intervals. The
calibration and test sets, the estimated random forest model, and every
pinterval_*() call are identical to the main text. The only
change in the replication code is the evaluation step which uses
interval_score() rather than
interval_coverage(). Because the interval score trades the
width of the intervals against miscoverage, we also report the mean
interval width of each method, computed with
interval_width().
We begin with the standard conformal, normal parametric, and bootstrapped intervals.
conformal_score <- interval_score(
truth = test$rel_dem_share,
lower_bound = conformal_intervals$lower_bound,
upper_bound = conformal_intervals$upper_bound,
alpha = 0.1
)
norm_score <- interval_score(
truth = test$rel_dem_share,
lower_bound = norm_intervals$lower_bound,
upper_bound = norm_intervals$upper_bound,
alpha = 0.1
)
bootstrap_score <- interval_score(
truth = test$rel_dem_share,
lower_bound = bootstrap_intervals$lower_bound,
upper_bound = bootstrap_intervals$upper_bound,
alpha = 0.1
)
c(conformal_score, norm_score, bootstrap_score)
#> [1] 0.3883673 0.3882974 0.3890984conformal_width <- interval_width(
lower_bound = conformal_intervals$lower_bound,
upper_bound = conformal_intervals$upper_bound
)
norm_width <- interval_width(
lower_bound = norm_intervals$lower_bound,
upper_bound = norm_intervals$upper_bound
)
bootstrap_width <- interval_width(
lower_bound = bootstrap_intervals$lower_bound,
upper_bound = bootstrap_intervals$upper_bound
)
c(conformal_width, norm_width, bootstrap_width)
#> [1] 0.3024480 0.3035879 0.3007931The three methods produce almost identical aggregate interval scores. This is not surprising, given the small difference in coverage rates noted in the main text. The mean interval widths are similarly close, ranging from 0.301 for the bootstrapped intervals to 0.304 for the normal parametric intervals.
We next score the MCP, CCP, and DWCP intervals.
mondrian_score <- interval_score(
truth = test$rel_dem_share,
lower_bound = mondrian_intervals$lower_bound,
upper_bound = mondrian_intervals$upper_bound,
alpha = 0.1
)
clustered_conformal_score <- interval_score(
truth = test$rel_dem_share,
lower_bound = clustered_conformal_intervals$lower_bound,
upper_bound = clustered_conformal_intervals$upper_bound,
alpha = 0.1
)
distance_weighted_cp_score <- interval_score(
truth = test$rel_dem_share,
lower_bound = distance_weighted_cp_intervals$lower_bound,
upper_bound = distance_weighted_cp_intervals$upper_bound,
alpha = 0.1
)
c(mondrian_score, clustered_conformal_score,
distance_weighted_cp_score)
#> [1] 0.3842136 0.3759636 0.3852191mondrian_width <- interval_width(
lower_bound = mondrian_intervals$lower_bound,
upper_bound = mondrian_intervals$upper_bound
)
clustered_conformal_width <- interval_width(
lower_bound = clustered_conformal_intervals$lower_bound,
upper_bound = clustered_conformal_intervals$upper_bound
)
distance_weighted_cp_width <- interval_width(
lower_bound = distance_weighted_cp_intervals$lower_bound,
upper_bound = distance_weighted_cp_intervals$upper_bound
)
c(mondrian_width, clustered_conformal_width,
distance_weighted_cp_width)
#> [1] 0.3028763 0.3046582 0.3016525Here too we see that the overall interval scores are similar, with a slightly lower (better) score for the CCP algorithm. As the CCP intervals are, on average, marginally wider than the MCP and DWCP intervals, the lower score of the CCP algorithm stems from smaller penalties for miscoverage rather than from sharper intervals.
To obtain the interval score within each US Census region, we reuse
the combined tibble all_conformal_intervals assembled in
the main text and summarize with interval_score() and
interval_width().
group_wise_score <- all_conformal_intervals %>%
group_by(method, region) %>%
summarise(
score = interval_score(
truth = rel_dem_share,
lower_bound = lower_bound,
upper_bound = upper_bound,
alpha = 0.1
),
width = interval_width(
lower_bound = lower_bound,
upper_bound = upper_bound
),
.groups = "drop"
)
group_wise_score %>%
group_by(method) %>%
mutate(
mean_score = mean(score),
mean_width = mean(width)
) %>%
ungroup() %>%
select(-width) %>%
pivot_wider(
names_from = region,
values_from = score
)
#> # A tibble: 6 × 7
#> method mean_score mean_width Midwest Northeast South West
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 Bootstrap 0.415 0.301 0.355 0.498 0.390 0.417
#> 2 CCP 0.387 0.302 0.363 0.397 0.371 0.415
#> 3 DWCP 0.420 0.303 0.365 0.528 0.370 0.419
#> 4 MCP 0.410 0.311 0.363 0.480 0.375 0.421
#> 5 Normal 0.422 0.304 0.365 0.523 0.375 0.426
#> 6 SCP 0.422 0.302 0.364 0.524 0.375 0.426The results clearly show that the CCP and MCP algorithms produce lower mean interval scores when aggregating across the US Census regions. The mean interval widths, in contrast, differ little between the methods, ranging from 0.301 to 0.311. The improvement of the group-aware algorithms therefore does not come from narrower intervals, but from adapting the intervals to each region.
Finally, we score the bin-conditional conformal intervals, both in aggregate and within each range of the outcome.
bccp_contiguized_score <- interval_score(
truth = test$rel_dem_share,
lower_bound = bccp_contiguized_intervals$lower_bound,
upper_bound = bccp_contiguized_intervals$upper_bound,
alpha = 0.1
)
bccp_discontigous_score <- interval_score(
truth = test$rel_dem_share,
intervals = bccp_discontigous_intervals$intervals,
alpha = 0.1
)
c(bccp_contiguized_score, bccp_discontigous_score)
#> [1] 0.4052939 0.3979040bccp_contiguized_width <- interval_width(
lower_bound = bccp_contiguized_intervals$lower_bound,
upper_bound = bccp_contiguized_intervals$upper_bound
)
bccp_discontigous_width <- interval_width(
intervals = bccp_discontigous_intervals$intervals
)
c(bccp_contiguized_width, bccp_discontigous_width)
#> [1] 0.3204990 0.3103355Both BCCP variants produce wider intervals than the other methods, with the discontiguous version being narrower than the contiguized one since it can exclude the gaps between non-adjacent bins from the prediction set.
bin_wise_score <- all_intervals_bins %>%
group_by(method, truth_bin) %>%
summarise(
score = interval_score(
truth = rel_dem_share,
lower_bound = lower_bound,
upper_bound = upper_bound,
intervals = intervals,
alpha = 0.1
),
width = interval_width(
lower_bound = lower_bound,
upper_bound = upper_bound,
intervals = intervals
),
.groups = "drop"
)
bin_wise_score %>%
group_by(method) %>%
mutate(
mean_score = mean(score),
mean_width = mean(width)
) %>%
ungroup() %>%
select(-width) %>%
pivot_wider(
names_from = truth_bin,
values_from = score
)
#> # A tibble: 5 × 7
#> method mean_score mean_width `1` `2` `3` `4`
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 BCCP(c) 0.405 0.321 0.405 0.364 0.401 0.451
#> 2 BCCP(d) 0.397 0.310 0.396 0.358 0.390 0.445
#> 3 Bootstrap 0.390 0.301 0.467 0.320 0.320 0.454
#> 4 Normal 0.387 0.304 0.423 0.318 0.323 0.483
#> 5 SCP 0.387 0.302 0.423 0.317 0.322 0.484The results indicate that the BCCP algorithm does not produce lower interval scores on average compared to the alternative methods. However, looking closely at the results, the BCCP algorithms seem to have an edge over the competing algorithms in both bin 1 and bin 4, and perform worse than the competing algorithms in bins 2 and 3. This is in line with expectations given how the interval score trades the interval width versus the miscoverage. As can be seen in the results for coverage, methods not taking the ranges of the prediction target into account tend to over-cover the middle of the distribution and under-cover the tails (i.e., bins 1 and 4). The BCCP algorithm achieves this by widening the prediction intervals, with mean interval widths of 0.321 and 0.310 for the contiguized and discontiguous versions, compared with 0.301 to 0.304 for the other methods. This leads to a higher interval score in the regions where other algorithms have narrower, over-covering, intervals. Conversely, because the BCCP algorithm achieves near the nominal coverage rate in bins 1 and 4, the interval score for these bins are lower for the BCCP methods. These results are further confirmed by the simulation study in Appendix B below.
This appendix presents a simulation study designed to systematically evaluate the empirical coverage properties of the prediction interval methods implemented in the pintervals package. While the main text illustrates these methods using a single train-calibration-test split, mimicking a real-world forecasting scenario, the simulation study assesses their performance across repeated random splits of the calibration and test data. Crucially, the training data remains constant throughout the simulation, only the calibration and test data are re-assigned by merging these two partitions and then splitting them into 50/50 calibration and test partitions. In total, we run 1000 simulations.
Across all simulation settings, the target coverage level is fixed at \(1 - \alpha = 0.9\). For each simulation, prediction intervals are constructed using the same procedures described in the main text, and performance is summarized using (i) empirical coverage, i.e., the proportion of observations in the test set whose true value is within the corresponding prediction interval, averaged across simulation runs, and (ii) the mean absolute error (MAE) of coverage relative to the target coverage level. The MAE metric is defined as
\[ \text{MAE} = \frac{1}{G} \sum_{g=1}^{G} \left| \hat{c}_g - (1 - \alpha) \right|, \]
where \(\hat{c}_g\) denotes empirical coverage in group \(g\), and \(G\) is the number of groups or bins under consideration. Lower MAE values indicate more uniform and locally valid coverage. The MAE is computed for each simulation iteration and then averaged across all runs to summarize performance.
In addition to the empirical coverage, we also evaluate the methods for generating prediction intervals using the interval score, using the same setup as in Appendix A.
The first set of simulations evaluates aggregate coverage, averaging across all observations in the test set without conditioning on subgroups or outcome ranges. For each simulation iteration, the data are randomly split into calibration and test sets, and prediction intervals are constructed using SCP, parametric prediction intervals, and bootstrap-based intervals. Here, we compare the SCP method against parametric intervals assuming normal distribution of errors, as well as bootstrapped intervals.
The results, in the table below, show that all three methods achieve nominal coverage very close to the nominal level of 0.9, with MAE values below 0.01. Overall, the differences between the methods are small although there is a minimal tendency towards overcoverage for the normal prediction intervals and towards undercoverage for the bootstrapped intervals.
| Method | Mean coverage | MAE coverage |
|---|---|---|
| Normal | 0.902 | 0.005 |
| SCP | 0.899 | 0.006 |
| Bootstrap | 0.896 | 0.007 |
The second set of simulations evaluates coverage within predefined subgroups of the data. These tests are motivated by the concern that methods with correct aggregate coverage may still exhibit systematic under- or over-coverage within substantively meaningful groups. In each simulation run, prediction intervals are constructed using SCP, MCP, CCP, and DWCP, as well as bootstrapped and parametric (normal) intervals.
Group membership is defined using geographic classifications, corresponding to US Census regions and divisions, as described in the main text. Coverage is computed separately within each group, while MAE of coverage is calculated across groups to quantify uniformity.
The results, in the table below, show that while all methods achieve empirical coverage close to the nominal level in aggregate, the group-aware CP methods (MCP, CCP) have group-wise coverage much closer to the nominal level compared to SCP, bootstrapped, and parametric intervals. In addition, these methods exhibit substantially lower MAE of coverage across regions compared to SCP, bootstrapped, and parametric intervals. This indicates that the group-aware methods provide more uniform and reliable coverage within subpopulations, addressing potential heterogeneity in the data.
| Method | Midwest | Northeast | South | West | Mean coverage | MAE coverage |
|---|---|---|---|---|---|---|
| MCP | 0.899 | 0.891 | 0.898 | 0.895 | 0.897 | 0.015 |
| CCP | 0.899 | 0.886 | 0.900 | 0.882 | 0.896 | 0.016 |
| DWCP | 0.916 | 0.818 | 0.890 | 0.896 | 0.895 | 0.031 |
| Normal | 0.927 | 0.826 | 0.901 | 0.879 | 0.902 | 0.032 |
| Bootstrap | 0.924 | 0.825 | 0.891 | 0.877 | 0.896 | 0.033 |
| SCP | 0.925 | 0.823 | 0.898 | 0.877 | 0.899 | 0.033 |
The third set of simulations focuses on coverage within different ranges of the outcome variable. This scenario is particularly relevant when prediction errors vary systematically across the outcome space, such as in the tails of the distribution. In each simulation iteration, the calibration data are partitioned into bins based on the observed outcome, using thresholds chosen to approximately balance the number of observations across four bins. Prediction intervals are then constructed using BCCP, with both contiguized and non-contiguized variants, alongside standard conformal, parametric, and bootstrap-based intervals. Coverage is then computed separately within each outcome bin. As in the grouped simulations, performance is summarized using mean aggregate coverage and the MAE of bin-wise coverage relative to the nominal level.
The results, in the table below, indicate that both contiguized and non-contiguized BCCP achieve empirical coverage closer to the nominal level within each outcome bin compared to SCP, bootstrapped, and parametric prediction intervals. The MAE of coverage across bins is lowest for the non-contiguized BCCP method, followed by the contiguized version. Both of the BCCP methods have substantially lower MAE across the bins, compared to the alternatives. These findings highlight the effectiveness of BCCP in addressing heterogeneity in prediction errors across the outcome space, leading to more reliable uncertainty quantification. As expected, the aggregate coverage for the contiguized BCCP is moderately higher than nominal due to the effect of contiguization.
| Method | Bin 1 | Bin 2 | Bin 3 | Bin 4 | Mean coverage | MAE coverage |
|---|---|---|---|---|---|---|
| BCCP(d) | 0.898 | 0.899 | 0.899 | 0.899 | 0.899 | 0.012 |
| BCCP(c) | 0.898 | 0.967 | 0.928 | 0.899 | 0.924 | 0.030 |
| Normal | 0.817 | 0.986 | 0.964 | 0.842 | 0.902 | 0.073 |
| SCP | 0.812 | 0.985 | 0.963 | 0.839 | 0.899 | 0.074 |
| Bootstrap | 0.797 | 0.983 | 0.962 | 0.842 | 0.896 | 0.077 |
In addition to empirical coverage, we evaluate the same simulated
prediction intervals using the interval score (defined in Appendix A). Lower values indicate sharper and
better-calibrated intervals. The reported values are averages over test
observations and across the 1000 simulation iterations, computed with
the interval_score() function.
The tables below report the aggregate interval score, the interval score by US Census region, and the interval score within ranges of the outcome variable.
The results confirm the findings from Appendix A, indicating that the group-aware versions of the CP algorithm generally outperform the aggregated versions when comparing the mean interval score across geographic groups. Similarly in agreement with the results from Appendix A, the BCCP methods outperform the alternatives on the interval score in the tail-end bins (bin 1 and 4), but underperform in the middle of the distribution.
| Method | Mean interval score |
|---|---|
| SCP | 0.384 |
| Normal | 0.384 |
| Bootstrap | 0.385 |
| Method | Midwest | Northeast | South | West | Mean interval score |
|---|---|---|---|---|---|
| CCP | 0.357 | 0.357 | 0.377 | 0.395 | 0.371 |
| MCP | 0.357 | 0.450 | 0.382 | 0.402 | 0.397 |
| DWCP | 0.358 | 0.480 | 0.376 | 0.401 | 0.404 |
| Normal | 0.360 | 0.481 | 0.381 | 0.405 | 0.407 |
| SCP | 0.359 | 0.482 | 0.381 | 0.405 | 0.407 |
| Bootstrap | 0.359 | 0.479 | 0.385 | 0.404 | 0.407 |
| Method | Bin 1 | Bin 2 | Bin 3 | Bin 4 | Mean interval score |
|---|---|---|---|---|---|
| SCP | 0.438 | 0.318 | 0.327 | 0.447 | 0.382 |
| Normal | 0.437 | 0.320 | 0.328 | 0.445 | 0.383 |
| Bootstrap | 0.450 | 0.318 | 0.326 | 0.443 | 0.384 |
| BCCP(d) | 0.420 | 0.369 | 0.375 | 0.399 | 0.391 |
| BCCP(c) | 0.433 | 0.371 | 0.387 | 0.405 | 0.399 |
Exchangeability, requiring only that the joint distribution of the data be invariant to permutations, is a slightly weaker assumption than the data being independent and identically distributed.↩︎
The calibration set is optional for parametric prediction intervals, as parameters can be supplied directly.↩︎
In classical MCP, the taxonomy is defined either by a
known categorical feature \(X\) or, in
classification tasks, by the label \(y\). Since pintervals
targets regression settings, the taxonomy in
pinterval_mondrian() and pinterval_ccp() is
defined by a categorical feature observed in both the calibration and
test data. Conditioning on the outcome is handled by BCCP, which
partitions the continuous outcome space into bins.↩︎
For more details on guidelines for how to divide the data, see Ding et al. (2023).↩︎
The relative vote share is calculated as the party’s share of the votes in the county minus the overall vote share for the party in the electorate.↩︎
The DWCP algorithm can also use feature based distances in order to weigh observations that are similar in the feature space rather than geographical space higher.↩︎
These bins were chosen to roughly match the quartiles of relative Democratic party vote shares in the calibration set. The low values of these quartiles are explained by the fact that more populous counties tend to have a higher democratic vote share, leading the relative vote share for the Democratic party to be lower than the national aggregate in most individual counties. In practice, the bins should be defined based on domain knowledge and the specific application at hand.↩︎
In practice, users can already approximate this by
running pinterval_bccp() separately within each
feature-defined group, which yields bin-conditional coverage within each
group, at the cost of fewer calibration points per group–bin cell.↩︎