Quantile coverage probability

Quantile coverage probability (QCP) evaluates the calibration of individual predictive quantiles. For a predicted quantile at nominal probability p, QCP is the empirical proportion of observations that are less than or equal to that predicted quantile.

Usage

qcp(
  obs,
  quantiles = NULL,
  levels = NULL,
  na.rm = TRUE,
  pred = NULL,
  predictive_sd = NULL,
  distribution = NULL
)

Arguments

obs

Numeric observation vector.

quantiles

Numeric matrix or data frame with observations in rows and predicted quantiles in columns. Supply this with levels, without another input representation.

levels

Strictly increasing quantile probabilities, one per column of quantiles. Values must lie strictly between zero and one. For generated quantiles, defaults to seq(0.05, 0.95, by = 0.05) and is sorted with duplicates removed.

na.rm

Logical; remove incomplete observation/quantile rows?

pred, predictive_sd

Optional numeric vectors of predictive means and predictive standard deviations, supplied together and of the same length as obs. Assumes normal predictive distributions. Non-missing predictive standard deviations must be finite and strictly positive.

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

Named numeric vector containing one QCP value for each supplied quantile level, on the probability scale from zero to one.

Details

$$ \mathrm{QCP}(p) = \frac{1}{n} \sum_{i=1}^{n} I(obs_i \leq q_{i,p}) $$

For a calibrated predictive distribution, QCP should be close to the nominal quantile probability p. For example, approximately 5% of observations should fall below predicted 0.05 quantiles, approximately 50% below predicted medians, and approximately 95% below predicted 0.95 quantiles.

Values above p indicate that observations fall below the predicted quantile more frequently than expected, whereas values below p indicate that they do so less frequently than expected.

QCP differs from prediction interval coverage probability (picp()). PICP evaluates the joint coverage of a lower and upper prediction-interval bound, whereas QCP evaluates each predictive quantile separately. QCP can therefore reveal asymmetric or one-sided miscalibration that may be hidden by apparently satisfactory central prediction-interval coverage.

Quantile calibration is generally most informative when evaluated across several probability levels rather than at a single quantile. A dense sequence of quantiles provides a more complete view of distributional calibration, although selected levels remain useful when only specific predictive quantiles are available.

Instead of explicit quantiles, supply predictive mean and standard deviation (pred and predictive_sd) for normal quantiles, or distribution for empirical quantiles (stats::quantile(), type = 7). These representations are mutually exclusive. Complete rows are selected across obs and all supplied predictive values, so every level uses the same cases. If na.rm = FALSE and any case is incomplete, or no complete cases remain, every QCP value is NA.

References

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

See also

gg_qcp(), picp(), gg_coverage(), pit(), gg_pit()

Examples

set.seed(123)

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

levels <- seq(0.05, 0.95, by = 0.05)

quantiles <- vapply(
  levels,
  function(p) pred + stats::qnorm(p) * predictive_sd,
  numeric(n)
)

qcp(
  obs,
  quantiles = quantiles,
  levels = levels
)
#>  0.05   0.1  0.15   0.2  0.25   0.3  0.35   0.4  0.45   0.5  0.55   0.6  0.65 
#> 0.042 0.086 0.138 0.178 0.220 0.272 0.334 0.392 0.440 0.490 0.556 0.608 0.658 
#>   0.7  0.75   0.8  0.85   0.9  0.95 
#> 0.706 0.746 0.796 0.846 0.896 0.942 
qcp(obs, pred = pred, predictive_sd = predictive_sd, levels = levels)
#>  0.05   0.1  0.15   0.2  0.25   0.3  0.35   0.4  0.45   0.5  0.55   0.6  0.65 
#> 0.042 0.086 0.138 0.178 0.220 0.272 0.334 0.392 0.440 0.490 0.556 0.608 0.658 
#>   0.7  0.75   0.8  0.85   0.9  0.95 
#> 0.706 0.746 0.796 0.846 0.896 0.942 
qcp(1:3, distribution = cbind(0:2, 1:3, 2:4), levels = c(0.25, 0.75))
#> 0.25 0.75 
#>    0    1