1211: Fix sample/quantile validation type checks and PIT histogram/rounding messages

NEWS.md

@@ -1,5 +1,6 @@
 # scoringutils (development version)
 
+- Fixed several validation and messaging issues (#1211): `as_forecast_sample()`, `as_forecast_quantile()` and `as_forecast_multivariate_sample()` now error at validation time when `observed` or `predicted` are not numeric, instead of failing later inside `score()`; the rounding warning in `as_forecast_quantile()` now correctly states that quantile levels are rounded to 9 digits; and `get_pit_histogram()` for quantile-based forecasts now displays its full warning message and actually falls back to the quantiles present in the forecast when requested quantiles are missing, instead of returning an empty or incorrect result.
 - 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/class-forecast-multivariate-sample.R

@@ -80,13 +80,15 @@ as_forecast_multivariate_sample.default <- function(data,
 #' @export
 #' @rdname assert_forecast
 #' @importFrom cli cli_abort qty
-#' @importFrom checkmate assert_subset
+#' @importFrom checkmate assert_subset assert_numeric
 #' @keywords validate-forecast-object
 assert_forecast.forecast_multivariate_sample <- function(
   forecast, forecast_type = NULL, verbose = TRUE, ...
 ) {
   assert_subset(c("sample_id", ".mv_group_id"), colnames(forecast))
   forecast <- assert_forecast_generic(forecast, verbose)
+  assert_numeric(forecast$observed, .var.name = "observed")
+  assert_numeric(forecast$predicted, .var.name = "predicted")
 
   # make sure that for every .mv_group_id, the number of samples per
   # forecast unit is the same

R/class-forecast-quantile.R

@@ -61,9 +61,9 @@ as_forecast_quantile.default <- function(data,
     cli_warn(
       "The {.code quantile_level} column in your data
       seems to have a rounding issue
-      (run {.code diff(sort(unique(data$quantile_level)))} to see this.
+      (run {.code diff(sort(unique(data$quantile_level)))} to see this).
       As {.code scoringutils} does not support arbitrarily fine quantile level
-      increments, we're going to run {.code round(x, digits = 10)} on
+      increments, we're going to run {.code round(x, digits = 9)} on
       the {.code quantile_level} column."
     )
     data$quantile_level <- round(data$quantile_level, digits = 9)
@@ -77,7 +77,7 @@ as_forecast_quantile.default <- function(data,
 
 #' @export
 #' @rdname assert_forecast
-#' @importFrom checkmate assert_subset
+#' @importFrom checkmate assert_subset assert_numeric
 #' @keywords validate-forecast-object
 assert_forecast.forecast_quantile <- function(
   forecast, forecast_type = NULL, verbose = TRUE, ...
@@ -86,6 +86,8 @@ assert_forecast.forecast_quantile <- function(
   forecast <- assert_forecast_generic(forecast, verbose)
   assert_forecast_type(forecast, actual = "quantile", desired = forecast_type)
   assert_numeric(forecast$quantile_level, lower = 0, upper = 1)
+  assert_numeric(forecast$observed, .var.name = "observed")
+  assert_numeric(forecast$predicted, .var.name = "predicted")
   return(invisible(NULL))
 }
 
@@ -232,7 +234,7 @@ get_pit_histogram.forecast_quantile <- function(forecast, num_bins = NULL,
   assert_numeric(breaks, lower = 0, upper = 1, null.ok = TRUE)
   forecast <- clean_forecast(forecast, copy = TRUE, na.omit = TRUE)
   forecast <- as.data.table(forecast)
-  present_quantiles <- unique(c(0, forecast$quantile_level, 1))
+  present_quantiles <- sort(unique(c(0, forecast$quantile_level, 1)))
   present_quantiles <- round(present_quantiles, 10)
 
   if (!is.null(breaks)) {
@@ -244,14 +246,22 @@ get_pit_histogram.forecast_quantile <- function(forecast, num_bins = NULL,
   }
   ## avoid rounding errors
   quantiles <- round(quantiles, 10)
-  diffs <- round(diff(quantiles), 10)
 
   if (length(setdiff(quantiles, present_quantiles)) > 0) {
-    cli_warn(
-      "Some requested quantiles are missing in the forecast. ",
-      "The PIT histogram will be based on the quantiles present in the forecast."
-    )
+    cli_warn(c(
+      "Some requested quantiles are missing in the forecast.",
+      i = "The PIT histogram will be based on the quantiles present in the
+      forecast."
+    ))
+    # fall back to the quantiles present in the forecast: keep the requested
+    # quantiles that are present; if none are (except 0 and 1), use all
+    # present quantiles
+    quantiles <- quantiles[quantiles %in% present_quantiles]
+    if (length(quantiles) <= 2) {
+      quantiles <- present_quantiles
+    }
   }
+  diffs <- round(diff(quantiles), 10)
 
   forecast <- forecast[quantile_level %in% quantiles]
   forecast[, quantile_coverage := (observed <= predicted)]

R/class-forecast-sample.R

@@ -54,14 +54,16 @@ as_forecast_sample.default <- function(data,
 
 #' @export
 #' @rdname assert_forecast
-#' @importFrom checkmate assert_subset
+#' @importFrom checkmate assert_subset assert_numeric
 #' @keywords validate-forecast-object
 assert_forecast.forecast_sample <- function(
   forecast, forecast_type = NULL, verbose = TRUE, ...
 ) {
   forecast <- assert_forecast_generic(forecast, verbose)
   assert_subset("sample_id", colnames(forecast))
   assert_forecast_type(forecast, actual = "sample", desired = forecast_type)
+  assert_numeric(forecast$observed, .var.name = "observed")
+  assert_numeric(forecast$predicted, .var.name = "predicted")
   return(invisible(NULL))
 }
 

tests/testthat/test-class-forecast-multivariate-sample.R

@@ -20,6 +20,26 @@ test_that("as_forecast_multivariate_sample() works as expected", {
   )
 })
 
+test_that("as_forecast_multivariate_sample() errors on non-numeric observed and predicted", {
+  df <- data.frame(
+    observed = as.character(c(5, 5, 5, 6, 6, 6)),
+    predicted = as.character(c(4, 5, 6, 5, 6, 7)),
+    sample_id = rep(1:3, 2),
+    location = rep(c("A", "B"), each = 3),
+    model = "m1"
+  )
+  expect_error(
+    as_forecast_multivariate_sample(df, joint_across = "location"),
+    "Must be of type 'numeric', not 'character'"
+  )
+
+  df$observed <- c(5, 5, 5, 6, 6, 6)
+  expect_error(
+    as_forecast_multivariate_sample(df, joint_across = "location"),
+    "Must be of type 'numeric', not 'character'"
+  )
+})
+
 test_that("as_forecast_multivariate_sample() creates expected structure", {
   test <- na.omit(data.table::copy(example_sample_continuous))
   data.table::setnames(test,

tests/testthat/test-class-forecast-quantile.R

@@ -190,11 +190,35 @@ test_that("as_forecast_quantile handles rounding issues correctly", {
     as_forecast_quantile(quantile_data),
     "rounding issue"
   )
+  # the warning should advertise the rounding that is actually applied
+  expect_warning(
+    as_forecast_quantile(quantile_data),
+    "digits = 9"
+  )
   expect_no_condition(
     score(quantile_forecast, metrics = c(wis = wis))
   )
 })
 
+test_that("as_forecast_quantile() errors on non-numeric observed and predicted", {
+  df <- data.frame(
+    observed = as.character(c(1, 1)),
+    predicted = as.character(c(0.5, 1.5)),
+    quantile_level = c(0.25, 0.75),
+    model = "m1"
+  )
+  expect_error(
+    as_forecast_quantile(df),
+    "Must be of type 'numeric', not 'character'"
+  )
+
+  df$observed <- c(1, 1)
+  expect_error(
+    as_forecast_quantile(df),
+    "Must be of type 'numeric', not 'character'"
+  )
+})
+
 
 
 
@@ -426,3 +450,45 @@ test_that("get_pit_histogram.forecast_quantile() works as expected", {
   # check printing works
   expect_output(print(pit_quantile))
 })
+
+test_that("get_pit_histogram.forecast_quantile() falls back to present quantiles", {
+  w <- capture_warnings(
+    res <- get_pit_histogram(example_quantile, num_bins = 7, by = "model")
+  )
+  # a single warning whose message includes both sentences
+  expect_length(w, 1)
+  expect_match(w, "Some requested quantiles are missing in the forecast")
+  expect_match(
+    w, "The PIT histogram will be based on the quantiles present in the forecast"
+  )
+
+  # none of the requested quantiles (except 0 and 1) are present, so the
+  # histogram should fall back to all quantiles present in the forecast
+  present <- sort(unique(na.omit(example_quantile)$quantile_level))
+  n_models <- length(unique(na.omit(example_quantile)$model))
+  expect_identical(nrow(res), n_models * (length(present) + 1L))
+
+  # densities integrate to one for every model
+  widths <- diff(c(0, present, 1))
+  integrals <- res[, sum(density * widths), by = "model"]$V1
+  expect_equal(integrals, rep(1, n_models))
+})
+
+test_that("get_pit_histogram.forecast_quantile() keeps requested quantiles that are present", {
+  w <- capture_warnings(
+    res <- get_pit_histogram(
+      example_quantile, breaks = c(0.25, 0.33, 0.5), by = "model"
+    )
+  )
+  # only the missing-quantiles warning, no recycling warnings
+  expect_length(w, 1)
+  expect_match(w, "Some requested quantiles are missing in the forecast")
+
+  # 0.33 is missing in the forecast and should be dropped from the bins
+  expect_identical(unique(res$bin), c("[0,0.25)", "[0.25,0.5)", "[0.5,1)"))
+
+  # densities integrate to one for every model
+  n_models <- length(unique(na.omit(example_quantile)$model))
+  integrals <- res[, sum(density * c(0.25, 0.25, 0.5)), by = "model"]$V1
+  expect_equal(integrals, rep(1, n_models))
+})

tests/testthat/test-class-forecast-sample.R

@@ -28,6 +28,25 @@ test_that("Running `as_forecast_sample()` twice returns the same object", {
   )
 })
 
+test_that("as_forecast_sample() errors on non-numeric observed and predicted", {
+  df <- data.frame(
+    observed = as.character(c(5, 5, 5)),
+    predicted = as.character(c(4, 5, 6)),
+    sample_id = 1:3,
+    model = "m1"
+  )
+  expect_error(
+    as_forecast_sample(df),
+    "Must be of type 'numeric', not 'character'"
+  )
+
+  df$observed <- c(5, 5, 5)
+  expect_error(
+    as_forecast_sample(df),
+    "Must be of type 'numeric', not 'character'"
+  )
+})
+
 
 # ==============================================================================
 # is_forecast_sample() # nolint: commented_code_linter