Plot prediction-interval calibration

Produces a prediction interval coverage probability (PICP) reliability plot, also known in geostatistics as an accuracy plot. The plot compares empirical prediction-interval coverage with the corresponding nominal coverage over one or more central prediction-interval levels.

Usage

gg_coverage(
  obs,
  lower = NULL,
  upper = NULL,
  level = NULL,
  pred = NULL,
  predictive_sd = NULL,
  levels = NULL,
  na.rm = TRUE,
  point_size = 3,
  line_width = 0.6,
  distribution = NULL
)

Arguments

obs: Numeric observation vector.
lower, upper: Named lists of lower and upper prediction-interval bounds. Names must represent nominal coverage levels such as "0.50" or "0.95". A single pair of numeric vectors is also accepted when level is supplied.
level: Nominal central prediction-interval coverage for a single numeric lower/upper pair. It must be NULL when named lists are supplied.
pred, predictive_sd: Optional numeric vectors of predictive means and predictive standard deviations. When supplied together, central prediction intervals are generated assuming normal predictive distributions. Do not also supply lower or upper.
levels: Nominal central prediction-interval coverage probabilities used when pred and predictive_sd, or distribution, are supplied. Values must lie strictly between zero and one. Defaults to every percentage from 1% to 99%.
na.rm: Logical; remove incomplete observation/interval combinations? With interval lists, only cases complete in obs and both bounds at every supplied level are used, so all levels share the same validation sample. If FALSE, any incomplete case makes coverage missing at every level. With predictive samples, a row missing any draw or obs is incomplete.
point_size: Positive numeric point size.
line_width: Positive numeric width of the 1:1 reference line.
distribution: Optional numeric matrix or data frame of equally weighted predictive samples: one row per observation and one column per predictive draw. Supply this instead of bounds or predictive means and standard deviations. At least one draw is required; infinite values are not allowed.

Value

A ggplot2 object. Its plotting data contain:

nominal: Nominal prediction-interval coverage probability.
picp: Empirical prediction interval coverage probability.

Details

For a nominal central prediction interval with coverage probability $p$, a calibrated predictive uncertainty model should contain approximately a proportion $p$ of independent validation observations. Consequently, empirical PICP should satisfy

$$\mathrm{PICP}(p) \approx p,$$

and a well-calibrated model should follow the dashed 1:1 reference line.

Points below the 1:1 line indicate under-coverage: fewer observations are contained in the prediction intervals than expected. This generally indicates prediction intervals that are too narrow and predictive uncertainty that is underestimated.

Points above the 1:1 line indicate over-coverage: more observations are contained in the prediction intervals than expected. This generally indicates prediction intervals that are wider than required and predictive uncertainty that is overestimated.

Evaluating several interval levels provides more information than evaluating a single PICP value because it shows how calibration changes across the predictive distribution. When pred and predictive_sd are supplied and levels is left NULL, central normal prediction intervals are evaluated from 1% to 99% nominal coverage. Alternatively, distribution generates central empirical intervals from equally weighted predictive samples, using stats::quantile() with type = 7. The same default levels apply. Supply only one input representation.

The accuracy-plot approach was developed for direct assessment of local uncertainty in geostatistics by Deutsch (1997) and subsequently applied to uncertainty evaluation in soil science by Goovaerts (2001). Similar reliability diagnostics have also been used in digital soil mapping, including Wadoux, Brus, and Heuvelink (2018).

The graphical departures from the 1:1 line can be summarized numerically with accuracy_plot_metrics(), which calculates the total absolute area between the empirical coverage curve and the reference line and separates this departure into over-coverage and under-coverage components.

Central PICP evaluates the joint coverage of lower and upper prediction-interval bounds. It therefore does not identify how non-coverage is distributed between the two tails. A model can have approximately correct central interval coverage while having too many observations below one bound and too few above the other. Use gg_qcp() or gg_pit() when tail-specific or distributional calibration is also of interest.

Calibration should also be distinguished from sharpness. Good coverage can be obtained using unnecessarily wide prediction intervals, so PICP calibration should generally be interpreted together with interval width or a proper scoring rule such as interval_score().

References

Deutsch, C. V. (1997). Direct assessment of local accuracy and precision. In E. Y. Baafi and N. A. Schofield (Eds.), Geostatistics Wollongong '96, pp. 115-125.

Goovaerts, P. (2001). Geostatistical modelling of uncertainty in soil science. Geoderma, 103, 3-26. doi:10.1016/S0016-7061(01)00067-2

Wadoux, A. M. J.-C., Brus, D. J. and Heuvelink, G. B. M. (2018). Accounting for non-stationary variance in geostatistical mapping of soil properties. Geoderma, 324, 138-147.

Schmidinger, J. and Heuvelink, G. B. M. (2023). Validation of uncertainty predictions in digital soil mapping. Geoderma, 437, 116585. doi:10.1016/j.geoderma.2023.116585

Examples

set.seed(123)

n <- 200
pred <- seq(0, 10, length.out = n)
predictive_sd <- rep(1, n)
obs <- stats::rnorm(n, mean = pred, sd = predictive_sd)

# Reliability curve generated directly from a normal predictive distribution
p_normal <- gg_coverage(
  obs,
  pred = pred,
  predictive_sd = predictive_sd
)

# Selected prediction intervals can also be supplied directly
lower <- list(
  `0.50` = pred + stats::qnorm(0.25) * predictive_sd,
  `0.90` = pred + stats::qnorm(0.05) * predictive_sd
)

upper <- list(
  `0.50` = pred + stats::qnorm(0.75) * predictive_sd,
  `0.90` = pred + stats::qnorm(0.95) * predictive_sd
)

p_intervals <- gg_coverage(obs, lower = lower, upper = upper)
p_samples <- gg_coverage(1:3, distribution = cbind(0:2, 1:3, 2:4),
                         levels = c(0.5, 0.9))
# Print p_normal, p_intervals or p_samples to display a plot.