1197: Fix sample metrics for vector input with a single observation

NEWS.md

@@ -1,5 +1,6 @@
 # scoringutils (development version)
 
+- Fixed the sample-based metrics `bias_sample()`, `ae_median_sample()`, `se_mean_sample()` and `mad_sample()` mishandling the documented vector input for a single observation (a scalar `observed` with `predicted` given as a vector of samples). `ae_median_sample()` and `se_mean_sample()` silently treated the samples as separate one-sample forecasts and returned wrong results, while `bias_sample()` and `mad_sample()` errored. All sample metrics now treat this input as one forecast with N samples, consistent with `crps_sample()`. Also corrected the integer bias formula in the `bias_sample()` documentation, which stated P_t(x_t + 1) instead of P_t(x_t - 1) (the code was correct) (#1197).
 - `summarise_scores()` now errors when `by` contains a metric column (e.g. `by = c("model", "wis")`). Previously, such calls silently returned an unsummarised table with duplicate column names, because the score column was used both as a grouping column and as a column to summarise. This complements the fix for the empty-metrics case in #1179 (#1204).
 - Fixed `bias_quantile()` returning wrong values when quantile levels were passed unsorted: predictions were reordered by quantile level but the quantile levels themselves were not, so predictions and levels became mispaired. Also fixed a crash ("argument is of length zero") when `na.rm = TRUE` removed all quantile levels on one side of the median; `bias_quantile()` now returns `NA` in this case, consistent with `na.rm = FALSE` (#1198).
 - Forecast validation now errors when the same forecast unit has conflicting observed values. Previously, such invalid data passed validation for quantile, sample, nominal, ordinal and multivariate sample forecasts and `score()` silently returned multiple wrong score rows for a single forecast (#1201).

R/metrics-sample.R

@@ -27,6 +27,23 @@ assert_input_sample <- function(observed, predicted) {
 }
 
 
+#' @title Ensure that predicted samples are a matrix
+#' @description
+#' Converts a vector of predictive samples (allowed as input when there is
+#' only a single observation) into a 1xN matrix so that downstream
+#' row-wise operations work as expected. Matrix input is returned unchanged.
+#' @param predicted A vector of size N (predictive samples for a single
+#'   observation) or an nxN matrix of predictive samples.
+#' @returns An nxN matrix of predictive samples.
+#' @keywords internal
+ensure_sample_matrix <- function(predicted) {
+  if (is.null(dim(predicted))) {
+    dim(predicted) <- c(1, length(predicted))
+  }
+  return(predicted)
+}
+
+
 #' @title Determine bias of forecasts
 #'
 #' @description
@@ -51,7 +68,7 @@ assert_input_sample <- function(observed, predicted) {
 #' For integer valued forecasts, Bias is measured as
 #'
 #' \deqn{
-#' B_t (P_t, x_t) = 1 - (P_t (x_t) + P_t (x_t + 1))
+#' B_t (P_t, x_t) = 1 - (P_t (x_t) + P_t (x_t - 1))
 #' }
 #'
 #' to adjust for the integer nature of the forecasts.
@@ -89,6 +106,7 @@ assert_input_sample <- function(observed, predicted) {
 bias_sample <- function(observed, predicted) {
 
   assert_input_sample(observed, predicted)
+  predicted <- ensure_sample_matrix(predicted)
   prediction_type <- get_type(predicted)
 
   # empirical cdf
@@ -141,8 +159,9 @@ bias_sample <- function(observed, predicted) {
 #' @export
 ae_median_sample <- function(observed, predicted) {
   assert_input_sample(observed, predicted)
+  predicted <- ensure_sample_matrix(predicted)
   median_predictions <- apply(
-    as.matrix(predicted), MARGIN = 1, FUN = median # this is row-wise
+    predicted, MARGIN = 1, FUN = median # this is row-wise
   )
   ae_median <- abs(observed - median_predictions)
   return(ae_median)
@@ -171,7 +190,8 @@ ae_median_sample <- function(observed, predicted) {
 
 se_mean_sample <- function(observed, predicted) {
   assert_input_sample(observed, predicted)
-  mean_predictions <- rowMeans(as.matrix(predicted))
+  predicted <- ensure_sample_matrix(predicted)
+  mean_predictions <- rowMeans(predicted)
   se_mean <- (observed - mean_predictions)^2
 
   return(se_mean)
@@ -225,7 +245,7 @@ logs_sample <- function(observed, predicted, ...) {
       )
     )
   }
-  scoringRules::logs_sample(
+  scoringRules::logs_sample( # nolint: namespace_linter.
     y = observed,
     dat = predicted,
     ...
@@ -257,7 +277,7 @@ logs_sample <- function(observed, predicted, ...) {
 dss_sample <- function(observed, predicted, ...) {
   assert_input_sample(observed, predicted)
 
-  scoringRules::dss_sample(
+  scoringRules::dss_sample( # nolint: namespace_linter.
     y = observed,
     dat = predicted,
     ...
@@ -308,19 +328,16 @@ dss_sample <- function(observed, predicted, ...) {
 crps_sample <- function(observed, predicted, separate_results = FALSE, ...) {
   assert_input_sample(observed, predicted)
 
-  crps <- scoringRules::crps_sample(
+  crps <- scoringRules::crps_sample( # nolint: namespace_linter.
     y = observed,
     dat = predicted,
     ...
   )
 
   if (separate_results) {
-    if (is.null(dim(predicted))) {
-      ## if `predicted` is a vector convert to matrix
-      dim(predicted) <- c(1, length(predicted))
-    }
+    predicted <- ensure_sample_matrix(predicted)
     medians <- apply(predicted, 1, median)
-    dispersion <- scoringRules::crps_sample(
+    dispersion <- scoringRules::crps_sample( # nolint: namespace_linter.
       y = medians,
       dat = predicted,
       ...
@@ -411,7 +428,9 @@ underprediction_sample <- function(observed, predicted, ...) {
 #' @keywords metric
 mad_sample <- function(observed = NULL, predicted, ...) {
 
-  assert_input_sample(rep(NA_real_, nrow(predicted)), predicted)
+  n <- if (is.matrix(predicted)) nrow(predicted) else 1
+  assert_input_sample(rep(NA_real_, n), predicted)
+  predicted <- ensure_sample_matrix(predicted)
 
   sharpness <- apply(predicted, MARGIN = 1, mad, ...)
   return(sharpness)
@@ -544,9 +563,7 @@ pit_histogram_sample <- function(observed,
     n_replicates, null.ok = (integers != "random"),
     .var.name = paste("n_replicates with `integers` = ", integers)
   )
-  if (is.vector(predicted)) {
-    predicted <- matrix(predicted, nrow = 1)
-  }
+  predicted <- ensure_sample_matrix(predicted)
 
   if (integers != "random" && !is.null(n_replicates)) {
     cli_warn("`n_replicates` is ignored when `integers` is not `random`")

tests/testthat/test-metrics-sample.R

@@ -206,6 +206,22 @@ test_that("bias_sample() approx equals bias_quantile() for many samples", {
 })
 
 
+test_that("bias_sample() works with a vector of samples", {
+  # integer case
+  bias_vector <- expect_no_condition(bias_sample(1, c(0, 2, 4)))
+  expect_equal( # nolint: expect_identical_linter
+    bias_vector,
+    bias_sample(1, matrix(c(0, 2, 4), nrow = 1))
+  )
+  # continuous case
+  bias_vector <- expect_no_condition(bias_sample(1.5, c(0.5, 2.5, 4.5)))
+  expect_equal( # nolint: expect_identical_linter
+    bias_vector,
+    bias_sample(1.5, matrix(c(0.5, 2.5, 4.5), nrow = 1))
+  )
+})
+
+
 # `ae_median_sample` ===========================================================
 test_that("ae_median_sample works", {
   observed <- rnorm(30, mean = 1:30)
@@ -215,6 +231,19 @@ test_that("ae_median_sample works", {
   expect_equal(ae, scoringutils) # nolint: expect_identical_linter # nolint: expect_identical_linter
 })
 
+test_that("ae_median_sample() works with a vector of samples for a single observation", {
+  ae <- expect_no_condition(ae_median_sample(1, c(0, 2, 4)))
+  expect_length(ae, 1)
+  expect_equal(ae, 1) # nolint: expect_identical_linter
+})
+
+# `se_mean_sample()` ===========================================================
+test_that("se_mean_sample() works with a vector of samples", {
+  se <- expect_no_condition(se_mean_sample(1, c(0, 2, 4)))
+  expect_length(se, 1)
+  expect_equal(se, 1) # nolint: expect_identical_linter
+})
+
 # `mad_sample()` ===============================================================
 test_that("function throws an error when missing 'predicted'", {
   predicted <- replicate(50, rpois(n = 10, lambda = 1:10))
@@ -224,6 +253,35 @@ test_that("function throws an error when missing 'predicted'", {
   )
 })
 
+test_that("mad_sample() works with a vector of samples", {
+  dispersion <- expect_no_condition(mad_sample(predicted = c(0, 2, 4)))
+  expect_equal(dispersion, mad(c(0, 2, 4))) # nolint: expect_identical_linter
+})
+
+# vector input =================================================================
+test_that("sample metrics give identical results for vector and 1-row matrix", {
+  observed <- 2.3
+  predicted <- rnorm(20, mean = 2)
+  predicted_matrix <- matrix(predicted, nrow = 1)
+
+  expect_identical(
+    ae_median_sample(observed, predicted),
+    ae_median_sample(observed, predicted_matrix)
+  )
+  expect_identical(
+    se_mean_sample(observed, predicted),
+    se_mean_sample(observed, predicted_matrix)
+  )
+  expect_identical(
+    bias_sample(observed, predicted),
+    bias_sample(observed, predicted_matrix)
+  )
+  expect_identical(
+    mad_sample(predicted = predicted),
+    mad_sample(predicted = predicted_matrix)
+  )
+})
+
 
 # ============================================================================ #
 # pit_histogram_sample() # nolint: commented_code_linter