diff --git a/NEWS.md b/NEWS.md index e1ee239eb..105b44620 100644 --- a/NEWS.md +++ b/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). - 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). diff --git a/R/metrics-sample.R b/R/metrics-sample.R index 6504d7409..6680d27ba 100644 --- a/R/metrics-sample.R +++ b/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`") diff --git a/man/bias_sample.Rd b/man/bias_sample.Rd index 1364ebc9d..eadeebc42 100644 --- a/man/bias_sample.Rd +++ b/man/bias_sample.Rd @@ -40,7 +40,7 @@ the fraction equal to \eqn{x_t}. 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. diff --git a/man/ensure_sample_matrix.Rd b/man/ensure_sample_matrix.Rd new file mode 100644 index 000000000..f9c873dd4 --- /dev/null +++ b/man/ensure_sample_matrix.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/metrics-sample.R +\name{ensure_sample_matrix} +\alias{ensure_sample_matrix} +\title{Ensure that predicted samples are a matrix} +\usage{ +ensure_sample_matrix(predicted) +} +\arguments{ +\item{predicted}{A vector of size N (predictive samples for a single +observation) or an nxN matrix of predictive samples.} +} +\value{ +An nxN matrix of predictive samples. +} +\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. +} +\keyword{internal} diff --git a/tests/testthat/test-metrics-sample.R b/tests/testthat/test-metrics-sample.R index 766e190a9..d47f35bd4 100644 --- a/tests/testthat/test-metrics-sample.R +++ b/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