1209: Skip subgroups with too few comparators in get_pairwise_comparisons()

NEWS.md

@@ -1,5 +1,6 @@
 # scoringutils (development version)
 
+- Fixed `get_pairwise_comparisons()` and `add_relative_skill()` aborting entirely when any subgroup defined by `by` had fewer than two comparators, without saying which subgroup was the problem. Such subgroups are now skipped with a warning that names them, results are returned for the remaining subgroups, and `add_relative_skill()` fills `NA` for the skipped subgroups. If no subgroup has at least two comparators, the previous error is kept (#1209).
 - 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/pairwise-comparisons.R

@@ -67,7 +67,9 @@
 #'   will be one relative skill score per distinct entry of the column selected
 #'   in `compare`. If further columns are given here, for example, `by =
 #'   "location"` with `compare = "model"`, then one separate relative skill
-#'   score is calculated for every model in every location.
+#'   score is calculated for every model in every location. Subgroups with
+#'   fewer than two comparators are skipped with a warning; if no subgroup
+#'   has at least two comparators, an error is thrown.
 #' @param metric A string with the name of the metric for which
 #'   a relative skill shall be computed. By default this is either "crps",
 #'   "wis" or "brier_score" if any of these are available.
@@ -217,8 +219,45 @@ get_pairwise_comparisons <- function(
   }
 
   # do the pairwise comparison -------------------------------------------------
-  # split data set into groups determined by 'by'
-  split_scores <- split(scores, by = by)
+  # split data set into groups determined by 'by'. `drop = TRUE` silently
+  # drops empty groups (e.g. from unused factor levels in a `by` column) so
+  # that only subgroups actually present in the data are compared.
+  split_scores <- split(scores, by = by, drop = TRUE)
+
+  # exclude groups with fewer than two comparators, as no pairwise comparison
+  # is possible there. Error if no group has enough comparators.
+  n_comparators <- vapply(
+    split_scores,
+    function(x) length(unique(x[[compare]])),
+    integer(1)
+  )
+  if (all(n_comparators < 2)) {
+    cli_abort(
+      c(`!` = "There are not enough comparators to do any comparison")
+    )
+  }
+  if (any(n_comparators < 2)) {
+    #nolint start: object_usage_linter
+    too_few <- vapply(
+      split_scores[n_comparators < 2],
+      function(x) {
+        toString(
+          paste0(by, "=", vapply(by, function(col) {
+            format(x[[col]][1])
+          }, character(1)))
+        )
+      },
+      character(1)
+    )
+    cli_warn(
+      c(
+        `!` = "Some groups have fewer than two comparators and were excluded
+        from the pairwise comparisons: {.val {too_few}}"
+      )
+    )
+    #nolint end
+    split_scores <- split_scores[n_comparators >= 2]
+  }
 
   results <- lapply(split_scores,
     FUN = function(scores) {
@@ -249,6 +288,11 @@ get_pairwise_comparisons <- function(
 #' that subgroup is managed from [pairwise_comparison_one_group()]. In order to
 #' actually do the comparison between two models over a subset of common
 #' forecasts it calls [compare_forecasts()].
+#' @param by Character vector with column names that define further grouping
+#'   levels for the pairwise comparisons. Unlike
+#'   [get_pairwise_comparisons()], this function does not skip subgroups with
+#'   fewer than two comparators: it operates on a single subgroup and throws
+#'   an error if `scores` contains fewer than two comparators.
 #' @inherit get_pairwise_comparisons params return
 #' @importFrom cli cli_abort
 #' @importFrom data.table setnames

tests/testthat/test-pairwise_comparison.R

@@ -448,6 +448,104 @@ test_that("get_pairwise_comparisons() throws errors with wrong inputs", {
   )
 })
 
+test_that("get_pairwise_comparisons() skips subgroups with fewer than two comparators", {
+  # only one model left for Deaths, but Cases still has three models
+  scores_sub <- scores_quantile[
+    !(target_type == "Deaths" & model != "EuroCOVIDhub-ensemble")
+  ]
+
+  expect_warning(
+    res <- get_pairwise_comparisons(scores_sub, by = "target_type"),
+    "fewer than two comparators"
+  )
+
+  # results are returned for the valid subgroup only
+  expect_identical(unique(res$target_type), "Cases")
+  expect_identical(nrow(res), 9L)
+
+  # results are identical to computing on the valid subgroup alone
+  res_cases <- get_pairwise_comparisons(
+    scores_quantile[target_type == "Cases"],
+    by = "target_type"
+  )
+  expect_equal(res, res_cases) # nolint: expect_identical_linter
+
+  # warning names the offending subgroup
+  expect_warning(
+    get_pairwise_comparisons(scores_sub, by = "target_type"),
+    "target_type=Deaths"
+  )
+
+  # if all subgroups have fewer than two comparators, an error is thrown
+  expect_error(
+    get_pairwise_comparisons(
+      scores_quantile[model == "EuroCOVIDhub-ensemble"],
+      by = "target_type"
+    ),
+    "not enough comparators"
+  )
+})
+
+test_that("get_pairwise_comparisons() silently drops empty subgroups from unused factor levels", {
+  # only one model left for Deaths, but Cases still has three models
+  scores_sub <- scores_quantile[
+    !(target_type == "Deaths" & model != "EuroCOVIDhub-ensemble")
+  ]
+
+  # same data, but with `target_type` as a factor with an unused level
+  scores_factor <- data.table::copy(scores_sub)
+  scores_factor[, target_type := factor(
+    target_type,
+    levels = c("Cases", "Deaths", "Hospitalisations")
+  )]
+  scores_nolevel <- data.table::copy(scores_factor)
+  scores_nolevel[, target_type := droplevels(target_type)]
+
+  # the empty subgroup from the unused level is dropped silently: the only
+  # warning names the real skipped subgroup, with no "NA" label
+  warnings <- testthat::capture_warnings(
+    res_factor <- get_pairwise_comparisons(scores_factor, by = "target_type")
+  )
+  expect_length(warnings, 1)
+  expect_match(warnings, "target_type=Deaths")
+  expect_false(any(grepl("NA", warnings, fixed = TRUE)))
+
+  # results are identical to the same data without the unused level
+  suppressWarnings(
+    res_nolevel <- get_pairwise_comparisons(scores_nolevel, by = "target_type")
+  )
+  expect_identical(
+    droplevels(res_factor$target_type),
+    droplevels(res_nolevel$target_type)
+  )
+  res_factor[, target_type := NULL]
+  res_nolevel[, target_type := NULL]
+  expect_identical(res_factor, res_nolevel)
+
+  # an unused factor level alone triggers no warning at all
+  scores_valid <- data.table::copy(scores_quantile)
+  scores_valid[, target_type := factor(
+    target_type,
+    levels = c("Cases", "Deaths", "Hospitalisations")
+  )]
+  expect_no_warning(
+    get_pairwise_comparisons(scores_valid, by = "target_type")
+  )
+})
+
+test_that("add_relative_skill() fills NA for subgroups with too few comparators", {
+  scores_sub <- scores_quantile[
+    !(target_type == "Deaths" & model != "EuroCOVIDhub-ensemble")
+  ]
+
+  expect_warning(
+    rs <- add_relative_skill(scores_sub, by = "target_type"),
+    "fewer than two comparators"
+  )
+  expect_true(all(is.na(rs[target_type == "Deaths"]$wis_relative_skill)))
+  expect_false(anyNA(rs[target_type == "Cases"]$wis_relative_skill))
+})
+
 test_that("pairwise_comparison_one_group() throws error with wrong inputs", {
   test <- data.table::copy(scores_sample_continuous)
   test <- test[, "model" := NULL]