predictset

A technical working paper for this package can be found here.

predictset is an R package for model-agnostic conformal prediction and distribution-free uncertainty quantification. It constructs prediction intervals (regression) and prediction sets (classification) with finite-sample coverage guarantees - no distributional assumptions required. Pass a formula, a fitted lm, glm or ranger, or wrap anything else (xgboost, keras, lightgbm) with make_model().

Installation

# Install from CRAN
install.packages("predictset")

# Or install the development version from GitHub
# install.packages("devtools")
devtools::install_github("charlescoverdale/predictset")

library(predictset)

x <- matrix(rnorm(500 * 3), ncol = 3)
y <- x[, 1] * 2 + rnorm(500)
x_new <- matrix(rnorm(50 * 3), ncol = 3)

# 90% prediction intervals around any model
result <- conformal_split(x, y, model = y ~ ., x_new = x_new, alpha = 0.10)
result$lower  # lower bounds
result$upper  # upper bounds

Changes in 0.4.0

This release fixes six defects that changed numerical results. The three in the table below are the ones most likely to have affected your work; the other three are under smaller fixes. If you have run any of them on an earlier version, re-run it.

Two changes that may need action on your side:

• randomize = TRUE is now the default for conformal_aps() and conformal_raps(), which is the method as published. Pass seed for reproducible sets, or randomize = FALSE for the deterministic (more conservative) variant.
• conformal_weighted() takes a new weights_new argument. Supply it for the exact procedure of Tibshirani et al. (2019), in which every test point receives its own conformal quantile.

Smaller fixes: Jackknife+ and CV+ now return ±Inf where the quantile index falls outside 1..n rather than clamping to the extreme order statistic; a fitted glm now works with the classification methods; coverage_by_bin() handles tied predictions; the seed argument no longer leaves the global random stream altered; Mondrian returns an unbounded interval for a group too small to calibrate rather than silently borrowing the pooled quantile. Full detail in NEWS.md.

What is conformal prediction?

Standard machine learning models produce point predictions: a single number for regression, a single class for classification. But in practice, you almost always need to know how uncertain that prediction is. Conformal prediction is a framework for wrapping any model in a layer of calibrated uncertainty quantification. Given a target coverage level (say 90%), it produces prediction intervals or prediction sets that are guaranteed to contain the true value at least 90% of the time, regardless of the underlying model or data distribution.

The key property is that this guarantee holds in finite samples. It’s not asymptotic, and it doesn’t require distributional assumptions. The only requirement is that the calibration data and test data are exchangeable (roughly: drawn from the same distribution). This makes conformal prediction fundamentally different from parametric confidence intervals, bootstrap intervals, or Bayesian credible intervals, all of which depend on modelling assumptions that may not hold.

predictset implements the main conformal methods from the recent literature (split conformal, Jackknife+, CV+, conformalized quantile regression for regression, and APS, RAPS, and LAC for classification) in a lightweight package with a single non-base dependency (cli; the rest are base R).

How does predictset compare to other packages?

predictset is designed to complement rather than compete with probably. If you’re working in the tidymodels ecosystem and only need regression intervals, probably integrates neatly with your workflow. predictset fills the gaps: classification methods (APS, RAPS, LAC), Jackknife+/CV+ for regression, and a model-agnostic interface that works with any model, not just tidymodels workflows.

conformalInference by Ryan Tibshirani was foundational research code, but it hasn’t been updated since 2019, isn’t on CRAN, and doesn’t cover classification.

Quick start

Regression: prediction intervals with coverage verification

library(predictset)

set.seed(42)
n <- 500
x <- matrix(rnorm(n * 5), ncol = 5)
y <- x[, 1] * 2 + x[, 2] + rnorm(n)
x_new <- matrix(rnorm(100 * 5), ncol = 5)
y_new <- x_new[, 1] * 2 + x_new[, 2] + rnorm(100)  # true values (for evaluation)

# Fit conformal intervals
result <- conformal_split(x, y, model = y ~ ., x_new = x_new, alpha = 0.10)
print(result)
#> ── Conformal Prediction Intervals (Split Conformal) ──
#> • Coverage target: "90%"
#> • Training: 250 | Calibration: 250 | Predictions: 100
#> • Conformal quantile: 1.876
#> • Median interval width: 3.7519

# Extract intervals as a data frame
head(data.frame(pred = result$pred, lower = result$lower, upper = result$upper))

# Verify coverage: should be ≥ 90%
coverage(result, y_new)

# Predict on new data later (fit once, predict many times)
future_data <- matrix(rnorm(50 * 5), ncol = 5)
predict(result, newdata = future_data)

Classification: prediction sets

set.seed(42)
n <- 400
x <- matrix(rnorm(n * 4), ncol = 4)
y <- factor(ifelse(x[, 1] + x[, 2] > 0, "A", "B"))
x_new <- matrix(rnorm(50 * 4), ncol = 4)

clf <- make_model(
  train_fun = function(x, y) glm(y ~ ., data = data.frame(y = y, x),
                                  family = "binomial"),
  predict_fun = function(object, x_new) {
    df <- as.data.frame(x_new)
    names(df) <- paste0("X", seq_len(ncol(x_new)))
    p <- predict(object, newdata = df, type = "response")
    cbind(A = 1 - p, B = p)
  },
  type = "classification"
)

# APS is randomised by default (Romano, Sesia & Candes 2020). Pass `seed`, or
# call set.seed(), for reproducible sets.
result <- conformal_aps(x, y, model = clf, x_new = x_new, alpha = 0.10, seed = 1)
result$sets[[1]]   # prediction set for the first observation, e.g. c("A")
table(set_size(result))  # distribution of set sizes

The uniform random variable in the APS score is not an optional refinement. Without it the score has an atom at exactly 1, and whenever the model ranks the true class last more often than alpha of the time the conformal quantile is exactly 1 and every set becomes the full label set. randomize = FALSE is available for deterministic output and warns when this happens.

Prediction sets are the exact inversion of the calibrated score, which can be empty for a test point the model is confident about. By default an empty set is replaced by the single most probable class. Pass allow_empty = TRUE to conformal_lac(), conformal_aps(), conformal_raps(), or conformal_mondrian_class() to get the inversion untouched: empty sets are the defining feature of LAC in Sadinle, Lei and Wasserman (2019).

Choosing a method

*Jackknife+ and CV+ have a theoretical coverage guarantee of 1-2α (Barber et al. 2021), weaker than split conformal’s 1-α. In practice, coverage is typically near 1-α. **ACI provides asymptotic (not finite-sample) coverage guarantees.

Model interface

There are three ways to specify a model. This flexibility means predictset works with anything from a simple linear model to a custom deep learning wrapper.

1. Formula shorthand (fits lm internally):

result <- conformal_split(x, y, model = y ~ ., x_new = x_new)

This is the quickest way to get started. Pass a formula and predictset handles the fitting.

2. Fitted model (auto-detected for lm, glm, and ranger):

fit <- lm(y ~ ., data = data.frame(y = y, x))
result <- conformal_split(x, y, model = fit, x_new = x_new)

If you’ve already fitted a model, pass it directly. predictset reads its formula and its call, then refits that same specification on each conformal training split. This matters: conformal prediction has to retrain on the split, so the object you pass is a template, not the model that ends up being used. Transformations in the formula (poly(), log(), interactions) and ranger hyperparameters are preserved. Anything predictset cannot refit raises an error pointing you at make_model() rather than silently substituting a default.

3. Custom model via make_model() (works with anything):

xgb_model <- make_model(
  train_fun = function(x, y) {
    dtrain <- xgboost::xgb.DMatrix(x, label = y)
    xgboost::xgb.train(list(objective = "reg:squarederror"), dtrain, nrounds = 100)
  },
  predict_fun = function(object, x_new) {
    predict(object, xgboost::xgb.DMatrix(x_new))
  },
  type = "regression"
)

result <- conformal_split(x, y, model = xgb_model, x_new = x_new)

make_model() takes a training function, a prediction function, and a type ("regression" or "classification"). This is how you use conformal prediction with xgboost, keras, lightgbm, or any other model.

Methods

More examples

Jackknife+ with ranger

Jackknife+ uses leave-one-out refitting to produce prediction intervals without splitting the data. This gives tighter intervals than split conformal (because it uses all the data for both training and calibration) at the cost of refitting the model n times.

set.seed(42)
n <- 200
x <- matrix(rnorm(n * 3), ncol = 3)
y <- x[, 1]^2 + x[, 2] + rnorm(n, sd = 0.5)
x_new <- matrix(rnorm(50 * 3), ncol = 3)

rf <- make_model(
  train_fun = function(x, y) {
    ranger::ranger(y ~ ., data = data.frame(y = y, x))
  },
  predict_fun = function(object, x_new) {
    predict(object, data = as.data.frame(x_new))$predictions
  },
  type = "regression"
)

result <- conformal_jackknife(x, y, model = rf, x_new = x_new, alpha = 0.10)
print(result)
plot(result)

Normalised conformal for heteroscedastic data

When the noise varies across the input space (e.g. predictions are more uncertain at extreme values), standard conformal intervals are too wide in low-noise regions and too narrow in high-noise ones. Normalised conformal scoring fixes this by scaling residuals by a local estimate of variability.

set.seed(42)
n <- 500
x <- matrix(runif(n, 0, 10), ncol = 1)
y <- sin(x[, 1]) + rnorm(n, sd = 0.1 + 0.3 * x[, 1])  # noise grows with x
x_new <- matrix(seq(0, 10, length.out = 100), ncol = 1)

result <- conformal_split(
  x, y, model = y ~ ., x_new = x_new, alpha = 0.10,
  score_type = "normalized"
)

# Intervals are narrower near x = 0, wider near x = 10
plot(result)

Mondrian conformal: group-conditional coverage

Standard conformal prediction guarantees marginal coverage (across all test points), but coverage can vary wildly across subgroups. Mondrian conformal computes a separate quantile for each group, guaranteeing coverage within each subgroup. This is critical for fairness and regulatory compliance.

To our knowledge no actively maintained CRAN package offers Mondrian conformal prediction for both regression and classification. (conformalClassification implemented label-wise Mondrian ICP but was archived from CRAN in May 2026.)

set.seed(42)
n <- 600
x <- matrix(rnorm(n * 3), ncol = 3)
groups <- factor(ifelse(x[, 1] > 0, "high", "low"))
y <- x[, 1] * 2 + ifelse(groups == "high", 3, 0.5) * rnorm(n)
x_new <- matrix(rnorm(200 * 3), ncol = 3)
groups_new <- factor(ifelse(x_new[, 1] > 0, "high", "low"))

result <- conformal_mondrian(x, y, model = y ~ ., x_new = x_new,
                              groups = groups, groups_new = groups_new)
print(result)

# Check per-group coverage
y_new <- x_new[, 1] * 2 + ifelse(groups_new == "high", 3, 0.5) * rnorm(200)
coverage_by_group(result, y_new, groups_new)

Weighted conformal: handling covariate shift

When the test distribution differs from the training distribution (covariate shift), standard conformal coverage guarantees break down. Weighted conformal prediction uses importance weights to correct for this shift.

set.seed(42)
n <- 500
x <- matrix(rnorm(n * 3), ncol = 3)
y <- x[, 1] * 2 + rnorm(n)
x_new <- matrix(rnorm(100 * 3, mean = 1), ncol = 3)  # shifted test data

# Importance weights (likelihood ratio of test vs training distributions),
# evaluated on both the calibration covariates and the test covariates
weights <- dnorm(x[, 1], mean = 1) / dnorm(x[, 1], mean = 0)
weights_new <- dnorm(x_new[, 1], mean = 1) / dnorm(x_new[, 1], mean = 0)

result <- conformal_weighted(x, y, model = y ~ ., x_new = x_new,
                              weights = weights, weights_new = weights_new)
print(result)

# Each test point gets its own conformal quantile
summary(result$quantile_by_point)

Conformalized Quantile Regression

CQR combines conformal prediction with quantile regression to produce intervals that naturally adapt to heteroscedasticity. Instead of fitting a model for the mean and adding symmetric bands, CQR fits models for the lower and upper quantiles and then adjusts them to guarantee coverage.

set.seed(42)
n <- 500
x <- matrix(rnorm(n * 3), ncol = 3)
y <- x[, 1] + x[, 2]^2 + rnorm(n, sd = 0.5 + abs(x[, 1]))
x_new <- matrix(rnorm(100 * 3), ncol = 3)

# In practice, use quantile regression (e.g. quantreg::rq).
# Here we approximate with shifted linear models for illustration.
model_lo <- make_model(
  train_fun = function(x, y) lm(y ~ ., data = data.frame(y = y, x)),
  predict_fun = function(obj, x_new) predict(obj, newdata = as.data.frame(x_new)) - 1.5,
  type = "regression"
)
model_hi <- make_model(
  train_fun = function(x, y) lm(y ~ ., data = data.frame(y = y, x)),
  predict_fun = function(obj, x_new) predict(obj, newdata = as.data.frame(x_new)) + 1.5,
  type = "regression"
)

result <- conformal_cqr(x, y, model_lo, model_hi, x_new = x_new, alpha = 0.10)
print(result)
plot(result)

Adaptive Conformal Inference (sequential prediction)

ACI adapts the miscoverage level online based on observed coverage, maintaining long-run coverage even under distribution shift. See also conformalForecast for adaptive conformal prediction and conformal PID control aimed specifically at multistep time series forecasting, and AdaptiveConformal (GitHub) for a wider family of online conformal algorithms.

set.seed(42)
n <- 500
y_true <- cumsum(rnorm(n, sd = 0.1)) + rnorm(n)  # drifting process
y_pred <- c(0, y_true[-n])                         # naive lag-1 predictor

result <- conformal_aci(y_pred, y_true, alpha = 0.10, gamma = 0.01)
print(result)
plot(result)  # intervals + adaptive alpha trace

Comparing methods

Benchmark multiple conformal methods side-by-side:

set.seed(42)
n <- 500
x <- matrix(rnorm(n * 3), ncol = 3)
y <- x[, 1] * 2 + rnorm(n)
x_new <- matrix(rnorm(200 * 3), ncol = 3)
y_new <- x_new[, 1] * 2 + rnorm(200)

comp <- conformal_compare(x, y, model = y ~ ., x_new = x_new, y_new = y_new,
                           methods = c("split", "cv", "jackknife"))
print(comp)

Diagnostics

After producing predictions, use the diagnostic functions to evaluate calibration and efficiency.

# Suppose y_test contains the true values for x_new
coverage(result, y_test)
#> [1] 0.91

# Average interval width (regression)
mean(interval_width(result))
#> [1] 3.42

# Prediction set sizes (classification)
table(set_size(result))
#>  1  2  3
#> 74 21  5

# Coverage within subgroups
coverage_by_group(result, y_test, groups = groups_test)
#>   group coverage   n target
#> 1  high    0.920  98    0.9
#> 2   low    0.891 102    0.9

# Coverage by prediction quantile bin
coverage_by_bin(result, y_test, bins = 5)

# Conformal p-values for outlier detection
pvals <- conformal_pvalue(result$scores, new_scores)

coverage() should be close to 1 - alpha from either side. Coverage substantially below target usually means exchangeability is violated; coverage far above it means the method is over-covering, which is a real cost, since it buys nothing and widens every interval. interval_width() and set_size() measure efficiency: narrower intervals and smaller sets are better, conditional on achieving the target coverage.

coverage_by_group() and coverage_by_bin() diagnose conditional coverage. Marginal coverage can mask severe under-coverage in subgroups (e.g. by demographic group or by prediction magnitude). These diagnostics help identify where intervals fail, and are essential for fairness evaluation.

Theory and references

Conformal prediction was introduced by Vovk, Gammerman, and Shafer in the early 2000s. The key insight is that if calibration and test data are exchangeable (i.e. their joint distribution is invariant to permutation), then the conformal p-value is uniformly distributed, which gives an exact finite-sample coverage guarantee. Unlike bootstrap or Bayesian intervals, this guarantee holds regardless of model misspecification.

The recent explosion of interest in conformal prediction has been driven by several methodological advances that make it practical for modern machine learning:

• Vovk, Gammerman, Shafer (2005). Algorithmic Learning in a Random World. Springer. The foundational book introducing conformal prediction.
• Lei, G’Sell, Rinaldo, Tibshirani, Wasserman (2018). Distribution-free predictive inference for regression. Journal of the American Statistical Association, 113(523), 1094–1111. Formalises split conformal prediction.
• Barber, Candes, Ramdas, Tibshirani (2021). Predictive inference with the Jackknife+. Annals of Statistics, 49(1), 486–507. Introduces Jackknife+ and CV+, which avoid the efficiency loss from data splitting.
• Romano, Patterson, Candes (2019). Conformalized quantile regression. NeurIPS 2019. Combines quantile regression with conformal calibration for adaptive intervals.
• Romano, Sesia, Candes (2020). Classification with valid and adaptive coverage. NeurIPS 2020. Introduces Adaptive Prediction Sets (APS) for classification.
• Angelopoulos, Bates, Malik, Jordan (2021). Uncertainty sets for image classifiers using conformal prediction. ICLR 2021. Introduces Regularized APS (RAPS) to reduce set sizes.
• Sadinle, Lei, Wasserman (2019). Least ambiguous set-valued classifiers with bounded error levels. Journal of the American Statistical Association, 114(525), 223–234. The LAC method for classification.
• Tibshirani, Barber, Candes, Ramdas (2019). Conformal prediction under covariate shift. NeurIPS 2019. Weighted conformal prediction for distribution shift.
• Gibbs, Candes (2021). Adaptive conformal inference under distribution shift. NeurIPS 2021. Online alpha adjustment for sequential prediction.

For an accessible introduction to the field, see Angelopoulos and Bates (2023), A Gentle Introduction to Conformal Prediction and Distribution-Free Uncertainty Quantification.

Limitations

• Split methods halve the training data. Split conformal, APS, RAPS, and LAC all divide the data into a training set and a calibration set. With small datasets, this can noticeably reduce model quality. Jackknife+ and CV+ avoid this at the cost of refitting the model multiple times.
• Jackknife+ and CV+ are computationally expensive. Jackknife+ refits the model n times; CV+ refits it K times. For large datasets or expensive models, this may be impractical.
• The coverage guarantee requires exchangeability. If the calibration data and test data come from different distributions (for example, if there is temporal drift) the coverage guarantee does not hold. This means conformal prediction is not directly applicable to time series forecasting without modification (e.g. conformal methods for time series exist but are not implemented here).
• Classification methods depend on probability estimates. APS, RAPS, and LAC require the model to output well-calibrated class probabilities. If the probabilities are poorly calibrated, the prediction sets will still have valid coverage but may be unnecessarily large.
• Small samples give unbounded intervals, by design. The conformal quantile is the ceiling((n+1)(1-alpha))-th calibration score. When that index exceeds the number of calibration points, the correct bound is infinite, and predictset returns Inf rather than silently substituting the largest observed score. At alpha = 0.10 this means fewer than 9 calibration points. The same applies per group in Mondrian conformal.
• Predictors must be numeric. Encode factors yourself, for example with stats::model.matrix(~ . - 1, data = x).
• Weighted conformal needs test-point weights. Supply weights_new to get the exact procedure of Tibshirani et al. (2019), in which each test point receives its own quantile.

Related packages

Issues

Found a bug or have a feature request? Please open an issue on GitHub.

Keywords

conformal prediction, prediction intervals, prediction sets, uncertainty quantification, coverage guarantee, distribution-free, model-agnostic, machine learning, statistics, R package