From 30021a98316b5e7e25f255e48fb7a2bf55014d3c Mon Sep 17 00:00:00 2001 From: Matthijs Berends Date: Fri, 16 Oct 2026 06:08:26 +0200 Subject: [PATCH] (v3.0.1.9153) Speed up mdro(), as.sir(), wisca() and wisca_plot() and the binding of mo vectors, fix PDR in the CMI 2012 guideline of mdro() being assigned to the wrong rows, and limit the stratified WISCA example to the 5 most common species --- DESCRIPTION | 4 +- NEWS.md | 6 ++- R/aa_helper_functions.R | 10 +++++ R/aa_helper_pm_functions.R | 5 ++- R/aa_session_cache.R | 24 ++++++++++ R/antibiogram.R | 48 ++++++++++++-------- R/mdro.R | 92 +++++++++++++++++++------------------- R/mo.R | 15 ++++--- R/sir.R | 6 ++- R/vctrs.R | 2 +- man/antibiogram.Rd | 4 +- tests/testthat/test-mdro.R | 16 +++++++ 12 files changed, 153 insertions(+), 79 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 17833a535..4208a4c55 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: AMR -Version: 3.0.1.9152 -Date: 2026-10-15 +Version: 3.0.1.9153 +Date: 2026-10-16 Title: Antimicrobial Resistance Data Analysis Description: Functions to simplify and standardise antimicrobial resistance (AMR) data analysis and to work with microbial and antimicrobial properties by diff --git a/NEWS.md b/NEWS.md index 187310856..24074302b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# AMR 3.0.1.9152 +# AMR 3.0.1.9153 Planned as v3.1.0, mid October 2026. @@ -80,6 +80,7 @@ Planned as v3.1.0, mid October 2026. * Microorganism coercion (`as.mo()`) and properties (`mo_*()`) * Considerably faster through one consistent session cache design, which is cleared when lookup tables change, e.g., after `add_custom_microorganisms()` (#313) * `mo_*()` functions now retrieve properties of valid `mo` input directly, and only for unique values + * Subsetting, combining, and binding rows of `mo` vectors (e.g., in `rbind()`) is now much faster, as codes are checked against a hashed set instead of the whole `microorganisms` data set * The mapping of outdated taxonomic names to their current names is now remembered per code * Previously coerced input is now recognised again, e.g., for an unknown species of a known genus, instead of being matched again on every call (#313) * `mo_uncertainties()` now also lists uncertainties of previously coerced input @@ -121,6 +122,7 @@ Planned as v3.1.0, mid October 2026. * Fixed a bug where both functions failed without the `antimicrobials` argument if `tidyselect` was not available, e.g., for binary installations without `dplyr` * Isolates of unknown species are now included in grouped antibiograms and in `wisca()`, as they already were in traditional antibiograms. They were previously excluded silently, with the warning "invalid microorganism code, NA generated". This can change WISCA coverage estimates. * Grouped traditional antibiograms now show the (translated) pathogen names, as ungrouped antibiograms do, instead of microorganism codes + * `wisca()` and `wisca_plot()` are considerably faster, especially with syndromic groups, by determining intrinsic resistance and pathogen names only once per regimen and pathogen * Resistance and susceptibility calculation * New argument `guideline` (default: the `AMR_guideline` option, otherwise EUCAST) in `susceptibility()`, `resistance()`, `count_df()`, `proportion_df()`, `sir_df()`, `ggplot_sir()`, `geom_sir()`, `labels_sir_count()`, `antibiogram()`, `wisca()`, and `bug_drug_combinations()` to interpret the 'I' category correctly per guideline. Functions that return results per category also accept `guideline = NULL` to keep S, I, and R separate, while functions that return a single proportion or count give an informative error. * When EUCAST is assumed because neither the `guideline` argument nor the `AMR_guideline` option is set, a message explains how to change this, once per session and not in R Markdown or Quarto documents @@ -145,6 +147,8 @@ Planned as v3.1.0, mid October 2026. * MDRO determination (`mdro()`) * Resistance can now be inferred from a drug/inhibitor combination when the column of the base drug is absent (e.g., piperacillin inferred from piperacillin/tazobactam). This is controlled with the new argument `infer_from_combinations` (default: `TRUE`) (#209). * BRMO classification now includes bacterial complexes (#275) + * Fixed a bug in the CMI 2012 guideline where pandrug-resistant (PDR) isolates were classified as extensively drug-resistant (XDR), and another isolate could be marked as PDR instead, if the isolate was not among the first rows of the data. This can change the PDR and XDR results of `guideline = "CMI 2012"` (the default). + * Considerably faster, as the classes per isolate are now determined per column instead of per cell * `top_n_microorganisms()` * New argument `property_for_each` for sub-grouping within the top *n* groups. Rank ordering is enforced, so that only lower taxonomic ranks are allowed. * `property = NULL` is now accepted diff --git a/R/aa_helper_functions.R b/R/aa_helper_functions.R index 67d4717a0..f38813b04 100755 --- a/R/aa_helper_functions.R +++ b/R/aa_helper_functions.R @@ -724,6 +724,16 @@ return_after_integrity_check <- function(value, type, check_vector) { value } +return_after_mo_integrity_check <- function(value) { + # same as return_after_integrity_check(), but with a hashed lookup of the codes, see mo_codes_valid() + valid <- is.na(value) | mo_codes_valid(value) + if (!all(valid)) { + warning_("invalid microorganism code, NA generated") + value[!valid] <- NA + } + value +} + # transforms data set to a tibble with only ASCII values, to comply with CRAN policies dataset_UTF8_to_ASCII <- function(df) { trans <- function(vect) { diff --git a/R/aa_helper_pm_functions.R b/R/aa_helper_pm_functions.R index 87675ebba..4298e7439 100755 --- a/R/aa_helper_pm_functions.R +++ b/R/aa_helper_pm_functions.R @@ -354,8 +354,11 @@ pm_group_data_worker <- function(.data, pm_groups) { class(res) <- "data.frame" nrow_res <- nrow(res) rows <- rep(list(NA), nrow_res) + # compute the group keys only once, instead of for every group + data_keys <- as.character(interaction(.data[, pm_groups])) + res_keys <- as.character(interaction(res[, pm_groups])) for (i in seq_len(nrow_res)) { - rows[[i]] <- which(interaction(.data[, pm_groups]) %in% interaction(res[i, pm_groups])) + rows[[i]] <- which(data_keys %in% res_keys[i]) } res$`.rows` <- rows res <- res[do.call(order, lapply(pm_groups, function(x) res[, x])), , drop = FALSE] diff --git a/R/aa_session_cache.R b/R/aa_session_cache.R index 7122b1565..66d6b33bb 100644 --- a/R/aa_session_cache.R +++ b/R/aa_session_cache.R @@ -101,7 +101,31 @@ reset_ab_cache <- function() { AMR_env$ab_previously_coerced <- new_coercion_cache() } +mo_codes_valid <- function(x) { + # whether codes are in `AMR_env$MO_lookup`; a hashed set avoids scanning the whole lookup table for every check, as `[<-.mo` runs one for every piece that rbind() binds + x <- as.character(x) + if (length(x) > 1000) { + # for long vectors, base matching is faster than looking up every unique value + add_MO_lookup_to_AMR_env() + return(x %in% as.character(AMR_env$MO_lookup$mo)) + } + if (is.null(AMR_env$MO_codes_hash)) { + add_MO_lookup_to_AMR_env() + codes <- unique(as.character(AMR_env$MO_lookup$mo)) + AMR_env$MO_codes_hash <- list2env(stats::setNames(as.list(rep(TRUE, length(codes))), codes), hash = TRUE, size = length(codes)) + } + ux <- unique(x) + valid <- vapply( + FUN.VALUE = logical(1), + ux, + function(code) !is.na(code) && nzchar(code) && exists(code, envir = AMR_env$MO_codes_hash, inherits = FALSE), + USE.NAMES = FALSE + ) + valid[match(x, ux)] +} + reset_mo_cache <- function() { + AMR_env$MO_codes_hash <- NULL AMR_env$MO_index <- NULL AMR_env$MO_translated <- NULL AMR_env$mo_accepted <- NULL diff --git a/R/antibiogram.R b/R/antibiogram.R index db4d56420..94f3e3359 100755 --- a/R/antibiogram.R +++ b/R/antibiogram.R @@ -192,8 +192,8 @@ #' wisca = TRUE #' ) #' -#' # stratified by syndrome or clinical group -#' out <- wisca(example_isolates, +#' # stratified by syndrome or clinical group, here for the 5 most common species to keep the plot readable +#' out <- wisca(top_n_microorganisms(example_isolates, n = 5), #' antimicrobials = c("TZP", "TZP+TOB", "TZP+GEN"), #' syndromic_group = "ward" #' ) @@ -1230,21 +1230,31 @@ create_wisca_priors <- function(data, sep) { beta_prior_beta <- rep(0.5, n_pathogens) # strongly informative Beta(1, 9999) for intrinsically resistant bug-drug pairs (Bielicki et al., 2016) - is_intrinsic <- vapply( - FUN.VALUE = logical(1), - seq_len(nrow(data)), - function(i) { - # split by " + ", or wherever `sep` is set to - ab_components <- as.ab(trimws(strsplit(as.character(data$ab[i]), trimws(sep), fixed = TRUE)[[1]])) - ab_components <- ab_components[!is.na(ab_components)] - length(ab_components) > 0 && - all(vapply( - FUN.VALUE = logical(1), - ab_components, - function(ab) any(AMR::intrinsic_resistant$mo == data$mo[i] & AMR::intrinsic_resistant$ab == ab) - )) - } - ) + # coerce every unique regimen only once and look up the bug-drug pairs as keys, which is much faster than per row + if (is.null(AMR_env$intrinsic_resistant_keys)) { + # the intrinsic_resistant data set never changes within a session + AMR_env$intrinsic_resistant_keys <- paste(as.character(AMR::intrinsic_resistant$mo), as.character(AMR::intrinsic_resistant$ab)) + } + regimens <- unique(as.character(data$ab)) + regimen_components <- lapply(regimens, function(regimen) { + # split by " + ", or wherever `sep` is set to + ab_components <- as.ab(trimws(strsplit(regimen, trimws(sep), fixed = TRUE)[[1]])) + as.character(ab_components[!is.na(ab_components)]) + }) + components <- regimen_components[match(as.character(data$ab), regimens)] + # look up all bug-drug pairs at once, and require all components of a regimen to be intrinsically resistant + pair_row <- rep(seq_len(nrow(data)), lengths(components)) + pair_intrinsic <- paste(as.character(data$mo)[pair_row], unlist(components)) %in% AMR_env$intrinsic_resistant_keys + is_intrinsic <- lengths(components) > 0 & !seq_len(nrow(data)) %in% pair_row[!pair_intrinsic] + for (i in which(is.na(data$mo))) { + # keep the original behaviour for a missing microorganism + is_intrinsic[i] <- length(components[[i]]) > 0 && + all(vapply( + FUN.VALUE = logical(1), + components[[i]], + function(ab) any(AMR::intrinsic_resistant$mo == data$mo[i] & AMR::intrinsic_resistant$ab == ab) + )) + } beta_prior_alpha[is_intrinsic] <- 1 beta_prior_beta[is_intrinsic] <- 9999 @@ -1635,7 +1645,9 @@ plot_wisca_susceptibility_incidence <- function(wisca_model, sep) { ) }) df <- do.call(rbind, df_list) - df$pathogen <- mo_shortname(df$pathogen, keep_synonyms = TRUE, info = FALSE) + # only translate the unique pathogens, as there is a row for every simulation + pathogens <- unique(df$pathogen) + df$pathogen <- mo_shortname(pathogens, keep_synonyms = TRUE, info = FALSE)[match(df$pathogen, pathogens)] # order pathogens by median incidence across all regimens med_inc <- tapply(df$incidence, df$pathogen, stats::median) diff --git a/R/mdro.R b/R/mdro.R index d2afebeca..4c115497e 100755 --- a/R/mdro.R +++ b/R/mdro.R @@ -743,6 +743,19 @@ mdro <- function(x = NULL, fluoroquinolones <- c(CIP, ENX, FLE, GAT, GEM, GRX, LVX, LOM, MFX, NOR, OFX, PAZ, PEF, PRU, RFL, SPX, TMX, TVA) # nolint end + # helper function returning a logical matrix (rows x cols) of which cells contain any of `values` + column_matches <- function(df, rows, cols, values) { + matrix( + vapply( + FUN.VALUE = logical(length(rows)), + cols, + function(col) as.character(df[rows, col, drop = TRUE]) %in% values + ), + nrow = length(rows), + dimnames = list(NULL, cols) + ) + } + # helper function for editing the table trans_tbl <- function(to, rows, cols, any_all, reason = NULL) { cols.bak <- cols @@ -759,19 +772,18 @@ mdro <- function(x = NULL, ), stringsAsFactors = FALSE ) + # determine the non-susceptible columns per column instead of per cell, which is much faster + group_vct <- unique(unname(cols_ab)) + nonsus <- column_matches(x, rows, group_vct, search_result) + current <- x[rows, "all_nonsusceptible_columns", drop = TRUE] x[rows, "all_nonsusceptible_columns"] <<- vapply( FUN.VALUE = character(1), - rows, - function(row, group_vct = cols_ab) { - cols_nonsus <- vapply( - FUN.VALUE = logical(1), - x[row, group_vct, drop = FALSE], - function(y) y %in% search_result - ) + seq_along(rows), + function(i) { paste( unique(sort(c( - unlist(strsplit(x[row, "all_nonsusceptible_columns", drop = TRUE], ", ", fixed = TRUE)), - names(cols_nonsus)[cols_nonsus] + unlist(strsplit(current[i], ", ", fixed = TRUE)), + group_vct[nonsus[i, ]] ))), collapse = ", " ) @@ -827,51 +839,39 @@ mdro <- function(x = NULL, ), stringsAsFactors = FALSE ) + # determine the matches per column instead of per cell, which is much faster + lst_cols <- unique(unname(lst_vector)) + available <- column_matches(x, rows, lst_cols, VALID_SIR_LEVELS[VALID_SIR_LEVELS != "NI"]) + nonsus <- column_matches(x, rows, lst_cols, search_result) + count_groups <- function(matches) { + rowSums(matrix( + vapply( + FUN.VALUE = logical(length(rows)), + lst, + function(group) { + group <- group[!is.na(group)] + rowSums(matches[, lst_cols %in% group, drop = FALSE]) > 0 + } + ), + nrow = length(rows) + )) + } + x[rows, "classes_in_guideline"] <<- length(lst) - x[rows, "classes_available"] <<- vapply( - FUN.VALUE = double(1), - rows, - function(row, group_tbl = lst) { - sum(vapply( - FUN.VALUE = logical(1), - group_tbl, - function(group) any(unlist(x[row, group[!is.na(group)], drop = TRUE]) %in% VALID_SIR_LEVELS[VALID_SIR_LEVELS != "NI"]) - )) - } - ) + x[rows, "classes_available"] <<- count_groups(available) if (isTRUE(verbose)) { x[rows, "all_nonsusceptible_columns"] <<- vapply( FUN.VALUE = character(1), - rows, - function(row, group_vct = lst_vector) { - cols_nonsus <- vapply(FUN.VALUE = logical(1), x[row, group_vct, drop = FALSE], function(y) y %in% search_result) - paste(unique(sort(names(cols_nonsus)[cols_nonsus])), collapse = ", ") - } + seq_along(rows), + function(i) paste(sort(lst_cols[nonsus[i, ]]), collapse = ", ") ) } - x[rows, "classes_affected"] <<- vapply( - FUN.VALUE = double(1), - rows, - function(row, group_tbl = lst) { - sum( - vapply( - FUN.VALUE = logical(1), - group_tbl, - function(group) { - any(unlist(x[row, group[!is.na(group)], drop = TRUE]) %in% search_result, na.rm = TRUE) - } - ), - na.rm = TRUE - ) - } - ) + x[rows, "classes_affected"] <<- count_groups(nonsus) # for PDR; all drugs are R (or I if combine_SI = FALSE) - x_transposed <- as.list(as.data.frame(t(x[rows, lst_vector, drop = FALSE]), - stringsAsFactors = FALSE - )) - row_filter <- vapply(FUN.VALUE = logical(1), x_transposed, function(y) all(y %in% search_result, na.rm = TRUE)) - x[which(row_filter), "classes_affected"] <<- 999 + row_filter <- rowSums(!nonsus) == 0 + # `row_filter` is relative to `rows`, so it must be mapped back to the rows of `x` + x[rows[row_filter], "classes_affected"] <<- 999 } if (isTRUE(info)) { diff --git a/R/mo.R b/R/mo.R index 5f709f939..576f1d40d 100755 --- a/R/mo.R +++ b/R/mo.R @@ -767,7 +767,8 @@ pillar_shaft.mo <- function(x, ...) { mo_cols <- NULL } - all_mos <- c(AMR_env$MO_lookup$mo, NA) + # as.character() prevents that c.mo() runs its integrity check on the full lookup table + all_mos <- c(as.character(AMR_env$MO_lookup$mo), NA) if (!all(x %in% all_mos) || (!is.null(df) && !all(unlist(df[, which(mo_cols), drop = FALSE]) %in% all_mos))) { # markup old mo codes @@ -887,7 +888,8 @@ print.mo <- function(x, print.shortnames = FALSE, ...) { } x <- as.character(x) names(x) <- x_names - if (!all(x %in% c(AMR_env$MO_lookup$mo, NA))) { + # as.character() prevents that c.mo() runs its integrity check on the full lookup table + if (!all(as.character(x) %in% c(as.character(AMR_env$MO_lookup$mo), NA))) { warning_( "Some MO codes are from another AMR package version. ", "Please update the MO codes with {.help [{.fun as.mo}](AMR::as.mo)}.", @@ -921,7 +923,8 @@ summary.mo <- function(object, ...) { #' @noRd as.data.frame.mo <- function(x, ...) { add_MO_lookup_to_AMR_env() - if (!all(x %in% c(AMR_env$MO_lookup$mo, NA))) { + # as.character() prevents that c.mo() runs its integrity check on the full lookup table + if (!all(as.character(x) %in% c(as.character(AMR_env$MO_lookup$mo), NA))) { warning_( "The data contains old MO codes (from another AMR package version). ", "Please update your MO codes with {.help [{.fun as.mo}](AMR::as.mo)}." @@ -959,7 +962,7 @@ as.data.frame.mo <- function(x, ...) { attributes(y) <- attributes(i) # must only contain valid MOs add_MO_lookup_to_AMR_env() - return_after_integrity_check(y, "microorganism code", as.character(AMR_env$MO_lookup$mo)) + return_after_mo_integrity_check(y) } #' @method [[<- mo #' @export @@ -969,7 +972,7 @@ as.data.frame.mo <- function(x, ...) { attributes(y) <- attributes(i) # must only contain valid MOs add_MO_lookup_to_AMR_env() - return_after_integrity_check(y, "microorganism code", as.character(AMR_env$MO_lookup$mo)) + return_after_mo_integrity_check(y) } #' @method c mo #' @export @@ -979,7 +982,7 @@ c.mo <- function(...) { y <- NextMethod() attributes(y) <- attributes(x) add_MO_lookup_to_AMR_env() - return_after_integrity_check(y, "microorganism code", as.character(AMR_env$MO_lookup$mo)) + return_after_mo_integrity_check(y) } #' @method unique mo diff --git a/R/sir.R b/R/sir.R index 3d4577b8b..3dbe7916d 100755 --- a/R/sir.R +++ b/R/sir.R @@ -1178,8 +1178,10 @@ get_guideline <- function(guideline, reference_data) { return(guideline) } guideline_param <- trimws2(toupper(guideline)) - latest_clsi <- rev(sort(subset(reference_data, guideline %like% "CLSI")$guideline))[1L] - latest_eucast <- rev(sort(subset(reference_data, guideline %like% "EUCAST")$guideline))[1L] + # only the unique guidelines are needed, which is much faster than subsetting the whole table + guidelines_available <- unique(reference_data$guideline) + latest_clsi <- rev(sort(guidelines_available[guidelines_available %like% "CLSI"]))[1L] + latest_eucast <- rev(sort(guidelines_available[guidelines_available %like% "EUCAST"]))[1L] guideline_param[guideline_param == "CLSI"] <- latest_clsi guideline_param[guideline_param == "EUCAST"] <- latest_eucast # like 'EUCAST2020', should be 'EUCAST 2020' diff --git a/R/vctrs.R b/R/vctrs.R index 9dc425111..e9cbf816a 100755 --- a/R/vctrs.R +++ b/R/vctrs.R @@ -160,7 +160,7 @@ vec_cast.character.mo <- function(x, to, ...) { } vec_cast.mo.character <- function(x, to, ...) { add_MO_lookup_to_AMR_env() - return_after_integrity_check(x, "microorganism code", as.character(AMR_env$MO_lookup$mo)) + return_after_mo_integrity_check(x) } # S3: disk ---- diff --git a/man/antibiogram.Rd b/man/antibiogram.Rd index 69da2a47c..44c6b947e 100644 --- a/man/antibiogram.Rd +++ b/man/antibiogram.Rd @@ -363,8 +363,8 @@ antibiogram(example_isolates, wisca = TRUE ) -# stratified by syndrome or clinical group -out <- wisca(example_isolates, +# stratified by syndrome or clinical group, here for the 5 most common species to keep the plot readable +out <- wisca(top_n_microorganisms(example_isolates, n = 5), antimicrobials = c("TZP", "TZP+TOB", "TZP+GEN"), syndromic_group = "ward" ) diff --git a/tests/testthat/test-mdro.R b/tests/testthat/test-mdro.R index 721812921..7e003c695 100755 --- a/tests/testthat/test-mdro.R +++ b/tests/testthat/test-mdro.R @@ -348,3 +348,19 @@ test_that("test-mdro.R", { "AMC, AMX (inferred from AMC/AXS), AXS, DAP, LNZ, TEC, VAN" ) }) + +test_that("test-mdro.R - PDR does not depend on the row position", { + # a pandrug-resistant S. aureus must be PDR regardless of the rows before it, which other microorganisms take + abx <- c("GEN", "RIF", "CPT", "OXA", "FOX", "CIP", "MFX", "SXT", "FUS", "VAN", "TEC", "TLV", "TGC", "CLI", "DAP", "ERY", "LNZ", "CHL", "FOS", "QDA", "TCY", "DOX", "MNO") + pdr_first <- data.frame(mo = as.mo(c("S. aureus", "E. coli", "E. coli", "E. coli"))) + pdr_last <- data.frame(mo = as.mo(c("E. coli", "E. coli", "E. coli", "S. aureus"))) + for (ab in abx) { + pdr_first[[ab]] <- as.sir(c("R", "S", "S", "S")) + pdr_last[[ab]] <- as.sir(c("S", "S", "S", "R")) + } + out_first <- suppressMessages(suppressWarnings(mdro(pdr_first, guideline = "CMI 2012", info = FALSE))) + out_last <- suppressMessages(suppressWarnings(mdro(pdr_last, guideline = "CMI 2012", info = FALSE))) + expect_identical(as.character(out_first[1]), "Pandrug-resistant (PDR)") + expect_identical(as.character(out_last[4]), "Pandrug-resistant (PDR)") + expect_identical(as.character(out_first[2:4]), as.character(out_last[1:3])) +})