Keep only EUCAST blocking rows that have an effect on as.sir() (#326)

The new internal blocking_row_is_needed() defines this: as.sir() must give another result without the row, so the definition follows the fallback path of as.sir() itself. The build removes all rows without an effect, and the new strict test checks all rows of the latest guideline, and for all years that a blocking row has a breakpoint to block.

This replaces the separate taxonomic comparison, handles the points of the review on this pruning, and keeps the sites of blocking rows, which `site` in as.sir() now uses.

R/sir.R

@@ -1348,6 +1348,44 @@ sir_site_rank <- function(x) {
   min(match(tokens, SIR_SITE_ORDER))
 }
 
+# Whether a blocking row has any effect, which is the definition used for clinical_breakpoints (#326): as.sir() must give another result without it. This takes one row number and the rows of one guideline, agent and method.
+blocking_row_is_needed <- function(i, breakpoints) {
+  row <- breakpoints[i, , drop = FALSE]
+  others <- breakpoints[-i, , drop = FALSE]
+  values <- c(others$breakpoint_S, others$breakpoint_R)
+  values <- values[!is.na(values)]
+  if (length(values) == 0) {
+    return(FALSE)
+  }
+  if (row$method == "MIC") {
+    input <- as.mic(sort(unique(pmax(0.001, c(0.001, 512, values / 2, values, values * 2)))))
+  } else {
+    input <- as.disk(sort(unique(pmin(50, c(6, 50, values - 1, values, values + 1)))))
+  }
+  sir <- function(reference_data, uti, site) {
+    as.character(suppressMessages(suppressWarnings(as.sir(input,
+      mo = as.character(row$mo), ab = as.character(row$ab), guideline = row$guideline, uti = uti, site = site,
+      reference_data = reference_data, info = FALSE
+    ))))
+  }
+  row_site <- if (is.na(row$site) || length(sir_site_tokens(row$site)) == 0) NA_character_ else c(NA_character_, paste(setdiff(sir_site_tokens(row$site), "screen"), collapse = ","))
+  row_site <- unique(row_site[is.na(row_site) | nzchar(row_site)])
+  for (uti in list(NA, FALSE, TRUE)) {
+    for (site in row_site) {
+      if (isTRUE(uti) && !is.na(site) && !grepl("uti", site, fixed = TRUE) && !isTRUE(row$uti)) {
+        next
+      }
+      if (isFALSE(uti) && !is.na(site) && grepl("uti", site, fixed = TRUE)) {
+        next
+      }
+      if (!identical(sir(breakpoints, uti, site), sir(others, uti, site))) {
+        return(TRUE)
+      }
+    }
+  }
+  FALSE
+}
+
 as_sir_method <- function(method_short,
                           method_long,
                           x,

data-raw/_reproduction_scripts/reproduction_of_clinical_breakpoints.R

@@ -753,6 +753,37 @@ if (any(!clinical_breakpoints$ab %in% AMR::antimicrobials$ab)) {
   stop("agent(s) missing from the antimicrobials data set: ", toString(unique(clinical_breakpoints$ab[!clinical_breakpoints$ab %in% AMR::antimicrobials$ab])))
 }
 clinical_breakpoints$ab <- as.ab(clinical_breakpoints$ab)
+
+# Blocking rows ----
+# A blocking row of EUCAST only stays if it has an effect: as.sir() must give another result without it (#326). This uses as.sir() itself, so that it follows the same fallback path. The most specific rows are compared first, so that of two blocking rows on one path, the row that is needed stays. The build stops if this leaves any row without an effect.
+is_pruned_group <- clinical_breakpoints$type == "human" & grepl("EUCAST 20(19|2[0-9])", clinical_breakpoints$guideline)
+blocking_key <- paste(clinical_breakpoints$guideline, clinical_breakpoints$ab, clinical_breakpoints$method)
+keys_with_blocking <- unique(blocking_key[is_pruned_group & is.na(clinical_breakpoints$breakpoint_S) & is.na(clinical_breakpoints$breakpoint_R)])
+keep_per_key <- parallel::mclapply(keys_with_blocking, function(key) {
+  rows <- which(is_pruned_group & blocking_key == key)
+  group <- clinical_breakpoints[rows, , drop = FALSE]
+  kept <- rep(TRUE, nrow(group))
+  for (i in order(group$rank_index)) {
+    if (is.na(group$breakpoint_S[i]) && is.na(group$breakpoint_R[i])) {
+      position <- sum(kept[seq_len(i)])
+      if (!AMR:::blocking_row_is_needed(position, group[kept, , drop = FALSE])) {
+        kept[i] <- FALSE
+      }
+    }
+  }
+  rows[!kept]
+}, mc.cores = 8L, mc.preschedule = FALSE)
+if (any(vapply(keep_per_key, inherits, logical(1), what = "try-error"))) {
+  print(keep_per_key[vapply(keep_per_key, inherits, logical(1), what = "try-error")][[1]])
+  stop("pruning blocking rows failed, see above")
+}
+rows_without_effect <- unlist(keep_per_key)
+message("NOTE: ", length(rows_without_effect), " blocking row(s) of EUCAST removed, since as.sir() gives the same results without them")
+if (length(rows_without_effect) > 0) {
+  clinical_breakpoints <- clinical_breakpoints[-rows_without_effect, , drop = FALSE]
+}
+rm(is_pruned_group, blocking_key, keys_with_blocking, keep_per_key, rows_without_effect)
+
 usethis::use_data(clinical_breakpoints, overwrite = TRUE, compress = "xz", version = 2)
 rm(clinical_breakpoints)
 devtools::load_all(".")

data-raw/_reproduction_scripts/reproduction_of_clinical_breakpoints_eucast.R

@@ -499,10 +499,7 @@ parse_bp <- function(txt) {
   t <- trimws(txt)
   # placeholders such as "-" (susceptibility testing not recommended), "IE"
   # or "Note" are kept as they are: the cleanup turns them into blocking rows
-  # Parenthesised values like "(2)" are EUCAST's "breakpoints in brackets"
-  # (see is_bracketed_bp()). The raw data retain these values and their
-  # EUCAST definition. They are excluded before clinical breakpoint cleanup,
-  # so they do not generate blocking rows or prevent PK-PD fallback.
+  # Parenthesised values like "(2)" are EUCAST's "breakpoints in brackets" (see is_bracketed_bp()). The raw data retain these values and their EUCAST definition. They are excluded before clinical breakpoint cleanup, so they do not generate blocking rows or prevent PK-PD fallback.
   t
 }
 
@@ -2389,10 +2386,10 @@ taxon_contains <- function(broad, narrow) {
     broad == as.character(as.mo(mo_order(narrow, keep_synonyms = TRUE), info = FALSE)) |
     broad %in% AMR::microorganisms.groups$mo_group[AMR::microorganisms.groups$mo == as.character(narrow)]
 }
-to_blocker <- function(df, reason) {
+to_blocker <- function(df, reason, site = NA_character_) {
   df |> mutate(
     is_blocker = TRUE, breakpoint_S = NA_character_, breakpoint_R = NA_character_,
-    site = NA_character_, note = blocking_note(reason, note)
+    site = site, note = blocking_note(reason, note)
   )
 }
 for (i in seq_len(nrow(eucast_note_rules))) {
@@ -2640,16 +2637,8 @@ scope_blockers <- breakpoints_eucast |>
   distinct(guideline, sheet, table_organism, ab, method, site, .keep_all = TRUE) |>
   inner_join(table_scope, by = c("guideline", "sheet", "table_organism"), relationship = "many-to-many") |>
   mutate(mo = table_mo, note = NA_character_) |>
-  select(-table_mo) |>
-  mutate(
-    is_blocker = TRUE,
-    breakpoint_S = NA_character_,
-    breakpoint_R = NA_character_,
-    note = blocking_note(
-      "EUCAST lists this agent in this table for other organisms only",
-      note
-    )
-  )
+  select(-table_mo)
+scope_blockers <- to_blocker(scope_blockers, "EUCAST lists this agent in this table for other organisms only", site = scope_blockers$site)
 
 breakpoints_eucast <- bind_rows(breakpoints_eucast, scope_blockers)
 
@@ -2775,12 +2764,7 @@ if (nrow(unresolved_mo) > 0 || nrow(unresolved_ab) > 0) {
 }
 
 # Pruning blocking rows ----
-# A blocking row is only needed where its organism has no breakpoint of
-# its own for the same agent, method and UTI status. Blocking rows for the
-# PK-PD table ("UNKNOWN") have no function. UTI status is compared separately,
-# so a UTI-only breakpoint does not remove a blocker for other infections.
-# For duplicate blockers with these keys, the first is kept. UTI blockers
-# retain their site to distinguish them from blockers for other infections.
+# A blocking row is not needed where its organism has a breakpoint of its own for the same agent, method and UTI status, and blocking rows for the PK-PD table ("UNKNOWN") have no function. UTI status is part of the comparison, so that a breakpoint for UTI only does not remove the blocking row for other infections. Of several blocking rows for the same organism, agent, method, UTI status and site, the first is kept. All other blocking rows without any effect are removed in reproduction_of_clinical_breakpoints.R, by comparing the results of as.sir() with and without them (#326).
 has_own_bp <- breakpoints_eucast |>
   filter(!is_blocker) |>
   distinct(guideline, type, host, mo, ab, method, uti)
@@ -2789,69 +2773,10 @@ breakpoints_eucast <- bind_rows(
   breakpoints_eucast |>
     filter(is_blocker, mo != "UNKNOWN") |>
     anti_join(has_own_bp, by = c("guideline", "type", "host", "mo", "ab", "method", "uti")) |>
-    distinct(guideline, type, host, mo, ab, method, uti, .keep_all = TRUE) |>
+    distinct(guideline, type, host, mo, ab, method, uti, site, .keep_all = TRUE) |>
     mutate(disk_dose = NA_character_)
 )
 
-# Keep only blockers that prevent inheritance of a broader breakpoint.
-real_bp <- breakpoints_eucast |>
-  filter(
-    !is_blocker,
-    !is.na(breakpoint_S) | !is.na(breakpoint_R)
-  ) |>
-  distinct(guideline, type, host, ab, method, uti, mo)
-
-blockers <- breakpoints_eucast |>
-  filter(is_blocker) |>
-  mutate(block_id = row_number())
-
-# Find potentially applicable broader breakpoints.
-candidates <- blockers |>
-  select(block_id, guideline, type, host, ab, method, uti,
-         block_mo = mo) |>
-  inner_join(
-    real_bp |>
-      rename(parent_mo = mo),
-    by = c(
-      "guideline", "type", "host",
-      "ab", "method", "uti"
-    ),
-    relationship = "many-to-many",
-    na_matches = "na"
-  ) |>
-  filter(block_mo != parent_mo)
-
-# Calculate each taxonomic relationship only once.
-taxon_pairs <- candidates |>
-  distinct(parent_mo, block_mo)
-
-taxon_pairs$is_broader <- mapply(
-  function(parent, child) {
-    if (is.na(parent) || is.na(child)) return(FALSE)
-    
-    parent == "UNKNOWN" ||
-      isTRUE(taxon_contains(parent, child))
-  },
-  as.character(taxon_pairs$parent_mo),
-  as.character(taxon_pairs$block_mo)
-)
-
-# Identify blockers with at least one applicable broader breakpoint.
-needed_blockers <- candidates |>
-  inner_join(
-    taxon_pairs |> filter(is_broader),
-    by = c("parent_mo", "block_mo"),
-    relationship = "many-to-many"
-  ) |>
-  distinct(block_id)
-
-breakpoints_eucast <- bind_rows(
-  breakpoints_eucast |> filter(!is_blocker),
-  blockers |>
-    semi_join(needed_blockers, by = "block_id") |>
-    select(-block_id)
-)
-
 breakpoints_eucast <- breakpoints_eucast |>
   distinct(guideline, type, host, ab, mo, method, site, breakpoint_S, .keep_all = TRUE) |>
   mutate(

tests/testthat/test-data.R

@@ -617,6 +617,36 @@ test_that("clinical_breakpoints: screening tests are under screening agents (#32
   expect_identical(unname(as.list(ery("ERY-S"))), unname(as.list(ery("ERY"))))
 })
 
+test_that("clinical_breakpoints: every blocking row of EUCAST has an effect (#326)", {
+  skip_on_cran()
+  cb <- as.data.frame(clinical_breakpoints, stringsAsFactors = FALSE)
+  cb <- cb[cb$type == "human" & grepl("EUCAST 20(19|2[0-9])", cb$guideline), , drop = FALSE]
+  is_blocking <- is.na(cb$breakpoint_S) & is.na(cb$breakpoint_R)
+  key <- paste(cb$guideline, cb$ab, cb$method)
+  # a blocking row needs a breakpoint to block, for all years
+  has_breakpoint <- key %in% key[!is_blocking]
+  wrong <- unique(key[is_blocking & !has_breakpoint])
+  expect_true(length(wrong) == 0, info = paste("no breakpoints to block:", toString(utils::head(wrong, 20))))
+  # and as.sir() must give another result without it, for the latest guideline; all years are checked during the build
+  latest <- paste0("EUCAST ", max(as.integer(gsub("[^0-9]", "", cb$guideline))))
+  for (k in unique(key[is_blocking & cb$guideline == latest])) {
+    group <- cb[key == k, , drop = FALSE]
+    for (i in which(is.na(group$breakpoint_S) & is.na(group$breakpoint_R))) {
+      expect_true(blocking_row_is_needed(i, group), info = paste(k, group$mo[i], group$site[i]))
+    }
+  }
+  # the definition itself, on a constructed set of rows
+  rows <- cb[cb$guideline == latest & as.character(cb$ab) == "CTX" & cb$method == "MIC" & as.character(cb$mo) == "B_VIBRI", , drop = FALSE][1, ]
+  blocker <- rows
+  blocker$mo <- as.mo("Vibrio fluvialis", info = FALSE)
+  blocker$rank_index <- 2.1
+  blocker$breakpoint_S <- NA_real_
+  blocker$breakpoint_R <- NA_real_
+  group <- rbind(rows, blocker)
+  expect_true(blocking_row_is_needed(2, group))
+  expect_false(blocking_row_is_needed(1, blocker))
+})
+
 test_that("clinical_breakpoints: EUCAST breakpoints stated in the Notes (#329)", {
   skip_on_cran()
   cb <- as.data.frame(clinical_breakpoints, stringsAsFactors = FALSE)