1202: Fix interval_coverage() floating point matching and NaN scores for quantile levels 0 and 1

NEWS.md

@@ -1,5 +1,6 @@
 # scoringutils (development version)
 
+- Fixed `interval_coverage()` erroring on quantile levels generated with `seq()` (e.g. `seq(0.05, 0.95, 0.05)`) because the required quantile levels were matched with an exact floating point comparison. Quantile levels are now rounded to 10 decimal places before matching, consistent with the rest of the package. Also fixed `wis()`, `interval_score()` and `quantile_score(weigh = FALSE)` returning `NaN` for forecasts that include the quantile levels 0 and 1 (which form a 100% prediction interval where alpha = 0). Scores are now finite when the observation falls inside the interval, restoring the identity between the WIS and the mean of the quantile scores; the unweighted scores return `Inf` when the observation falls outside a 100% prediction interval (#1202).
 - Added `filter_scores()` and `impute_missing_scores()` for handling missing forecasts before summarisation. `filter_scores()` removes target combinations with insufficient model coverage, while `impute_missing_scores()` fills in missing scores using configurable strategies (worst, mean, NA, or reference model). Both use a strategy function pattern for extensibility. See `vignette("handling-missing-forecasts")` for details (#1122).
 - Added `plot_discrimination()` to visualise the discrimination ability of binary forecasts by plotting the distribution of predicted probabilities, stratified by the observed outcome. The function requires a `forecast_binary` object (created with `as_forecast_binary()`) (#942).
 - Fixed `summarise_scores()` producing a data.table with duplicate column names when the input `scores` object had no score columns (e.g. because every metric in `score()` warned and returned nothing). `summarise_scores()` now matches metric columns by exact name rather than regex partial match, and errors with a clear message when there is nothing to summarise (#1179).

R/metrics-interval-range.R

@@ -163,17 +163,28 @@ interval_score <- function(observed,
   # calculate alpha from the interval range
   alpha <- (100 - interval_range) / 100
 
-  # calculate three components of WIS
+  # calculate three components of WIS. The penalties are computed per unit
+  # (indicator applied first, without the 2 / alpha factor) so that quantile
+  # levels 0 and 1 (interval_range = 100, i.e. alpha = 0) yield finite scores
+  # instead of Inf * 0 = NaN
   dispersion <- (upper - lower)
-  overprediction <-
-    2 / alpha * (lower - observed) * as.numeric(observed < lower)
-  underprediction <-
-    2 / alpha * (observed - upper) * as.numeric(observed > upper)
+  overprediction <- (lower - observed) * as.numeric(observed < lower)
+  underprediction <- (observed - upper) * as.numeric(observed > upper)
 
   if (weigh) {
+    # the weight alpha / 2 cancels analytically with the 2 / alpha factor of
+    # the penalties, so only the dispersion needs to be weighted
     dispersion <- dispersion * alpha / 2
-    underprediction <- underprediction * alpha / 2
-    overprediction <- overprediction * alpha / 2
+  } else {
+    # zero penalties stay zero; otherwise apply the 2 / alpha factor, which
+    # correctly yields Inf when the observation falls outside a 100%
+    # prediction interval
+    overprediction <- ifelse(
+      overprediction == 0, 0, 2 / alpha * overprediction
+    )
+    underprediction <- ifelse(
+      underprediction == 0, 0, 2 / alpha * underprediction
+    )
   }
 
   score <- dispersion + underprediction + overprediction

R/metrics-quantile.R

@@ -103,6 +103,16 @@ assert_input_quantile <- function(observed, predicted, quantile_level,
 #' for an increasing number of equally spaced quantiles, the WIS
 #' converges to the continuous ranked probability score (CRPS).
 #'
+#' Note that the equivalence between the WIS and the mean of the quantile
+#' scores (see [quantile_score()]) only holds for this default weighting
+#' (`weigh = TRUE`). With `weigh = FALSE` and the quantile levels 0 and 1
+#' present, the two differ: the 0%-100% prediction interval contributes its
+#' full dispersion to the unweighted interval score, whereas the quantile
+#' levels 0 and 1 contribute a score of 0 to the unweighted quantile scores.
+#' For example, `wis(5, matrix(c(1, 3, 5, 7, 9), nrow = 1),
+#' c(0, 0.25, 0.5, 0.75, 1), weigh = FALSE)` returns 4.8, while the
+#' corresponding call to [quantile_score()] returns 1.6.
+#'
 #' **Quantile score**
 #'
 #' In addition to the interval score, there also exists a quantile score (QS)
@@ -335,7 +345,10 @@ interval_coverage <- function(observed, predicted,
     (100 - interval_range) / 2,
     100 - (100 - interval_range) / 2
   ) / 100
-  if (!all(necessary_quantiles %in% quantile_level)) {
+  # round before matching to avoid floating point issues, e.g. with
+  # quantile levels generated by seq() (consistent with the 10 decimal
+  # places convention used elsewhere in the package)
+  if (!all(round(necessary_quantiles, 10) %in% round(quantile_level, 10))) {
     #nolint start: object_usage_linter
     cli_abort(
       c(
@@ -613,7 +626,10 @@ ae_median_quantile <- function(observed, predicted, quantile_level) {
 #'
 #' `quantile_score()` returns the average quantile score across the quantile
 #' levels provided. For a set of quantile levels that form pairwise central
-#' prediction intervals, the quantile score is equivalent to the interval score.
+#' prediction intervals, the quantile score is equivalent to the interval score
+#' for the default weighting (`weigh = TRUE`). With `weigh = FALSE`, the two
+#' can differ when the quantile levels 0 and 1 are present (see [wis()] for
+#' details).
 #' @returns Numeric vector of length n with the quantile score. The scores are
 #' averaged across quantile levels if multiple quantile levels are provided
 #' (the result of calling `rowMeans()` on the matrix of quantile scores that
@@ -689,7 +705,11 @@ quantile_score <- function(observed,
     alpha <- 1 - central_interval
 
     # unweighted score' = 2 * score / alpha
-    score <- 2 * sweep(score, 2, alpha, FUN = "/")
-    return(rowMeans(score))
+    unweighted <- 2 * sweep(score, 2, alpha, FUN = "/")
+    # a weighted score of zero corresponds to an unweighted score of zero.
+    # Setting it explicitly avoids 0 / 0 = NaN for quantile levels 0 and 1
+    # (where alpha is 0)
+    unweighted[!is.na(score) & score == 0] <- 0
+    return(rowMeans(unweighted))
   }
 }

tests/testthat/test-metrics-quantile.R

@@ -619,6 +619,81 @@ test_that("Quantlie score and interval score yield the same result, weigh = TRUE
   }
 })
 
+test_that("wis works with quantile levels 0 and 1", {
+  predicted <- matrix(c(1, 3, 5, 7, 9), nrow = 1)
+  quantile_level <- c(0, 0.25, 0.5, 0.75, 1)
+
+  # observation inside the 100% interval
+  expect_equal( # nolint: expect_identical_linter
+    wis(5, predicted, quantile_level),
+    0.4
+  )
+  expect_equal( # nolint: expect_identical_linter
+    wis(5, predicted, quantile_level),
+    quantile_score(5, predicted, quantile_level)
+  )
+
+  # observation outside the 100% interval
+  expect_equal( # nolint: expect_identical_linter
+    wis(10, predicted, quantile_level),
+    quantile_score(10, predicted, quantile_level)
+  )
+
+  # components must be finite when the observation is inside the interval
+  separate <- wis(
+    5, predicted, quantile_level,
+    separate_results = TRUE
+  )
+  expect_equal(separate$overprediction, 0) # nolint: expect_identical_linter
+  expect_equal(separate$underprediction, 0) # nolint: expect_identical_linter
+})
+
+test_that("interval_score handles interval_range = 100", {
+  # weighted (alpha / 2 = 0): all components are zero
+  weighted <- interval_score(
+    observed = 5, lower = 1, upper = 9,
+    interval_range = 100, weigh = TRUE, separate_results = TRUE
+  )
+  expect_equal( # nolint: expect_identical_linter
+    weighted,
+    list(
+      interval_score = 0, dispersion = 0,
+      underprediction = 0, overprediction = 0
+    )
+  )
+
+  # unweighted, observation inside: dispersion only
+  expect_equal( # nolint: expect_identical_linter
+    interval_score(
+      observed = 5, lower = 1, upper = 9,
+      interval_range = 100, weigh = FALSE
+    ),
+    8
+  )
+
+  # unweighted, observation outside a 100% interval: infinite penalty
+  unweighted <- interval_score(
+    observed = 10, lower = 1, upper = 9,
+    interval_range = 100, weigh = FALSE, separate_results = TRUE
+  )
+  expect_identical(unweighted$underprediction, Inf)
+  expect_equal(unweighted$overprediction, 0) # nolint: expect_identical_linter
+})
+
+test_that("quantile_score(weigh = FALSE) handles quantile levels 0 and 1", {
+  predicted <- matrix(c(1, 9), nrow = 1)
+  # observation inside: zero penalty at levels 0 and 1, not NaN
+  expect_equal( # nolint: expect_identical_linter
+    quantile_score(5, predicted, c(0, 1), weigh = FALSE),
+    0
+  )
+  # observation outside a 100% prediction interval: infinite penalty
+  expect_identical(
+    quantile_score(10, predicted, c(0, 1), weigh = FALSE),
+    Inf
+  )
+})
+
 test_that("wis works with separate results", {
   wis <- wis(
     observed = y,
@@ -658,6 +733,29 @@ test_that("interval_coverage works", {
   )
 })
 
+test_that("interval_coverage handles floating point quantile levels from seq()", {
+  # seq() produces quantile levels that are not exactly representable
+  # (e.g. 0.35 %in% seq(0.05, 0.95, 0.05) is FALSE), which must not
+  # trigger the missing-quantile error
+  quantile_level <- seq(0.05, 0.95, 0.05)
+  predicted <- t(vapply(
+    c(1, 0, 22),
+    function(x) x + qnorm(quantile_level),
+    numeric(length(quantile_level))
+  ))
+  expect_no_condition(
+    interval_coverage(observed, predicted, quantile_level, interval_range = 30)
+  )
+  expect_equal( # nolint: expect_identical_linter
+    interval_coverage(observed, predicted, quantile_level, interval_range = 30),
+    c(TRUE, FALSE, TRUE)
+  )
+  expect_equal( # nolint: expect_identical_linter
+    interval_coverage(observed, predicted, quantile_level, interval_range = 50),
+    c(TRUE, FALSE, TRUE)
+  )
+})
+
 test_that("interval_coverage rejects wrong inputs", {
   expect_error(
     interval_coverage(observed, predicted, quantile_level, interval_range = c(50, 0)),