From 8a00987f251c763e57b7b187fb7f96521b3bb61c Mon Sep 17 00:00:00 2001 From: Emily Boyer Date: Fri, 16 Oct 2026 10:49:42 -0600 Subject: [PATCH 01/14] docs: add MDR design spec --- design/mdr.md | 178 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 178 insertions(+) create mode 100644 design/mdr.md diff --git a/design/mdr.md b/design/mdr.md new file mode 100644 index 0000000..13ce9b2 --- /dev/null +++ b/design/mdr.md @@ -0,0 +1,178 @@ +# MDR: modeling combinations of resistant classes + +Status: draft spec. The MDR yes/no target is implemented; the combination target is not yet. + +## Goal + +Two MDR models, both going through the same stages as every other target +(profile, eligibility, matrices, fitting) and running in the rewrite's modes, +not only on the whole dataset: + +1. **Combination model** (`target = "resistant_classes"`): port of old amRml's + MDR model (`.parquet2MDRMatrix()` in `R/generate_matrices_ml.R` on `main`). + Each combination of resistant classes is one label of a multi-class model. +2. **MDR yes/no model** (`target = "mdr"`): Resistant if the genome is resistant + to **2 or more** drug classes, Susceptible otherwise. The lab chose 2, not + the clinical standard of 3 (Magiorakos et al., 2012). + +Multi-label models (one yes/no per drug, nested in classes), random forests and +neural networks are later work and out of scope here. + +## What old amRml did + +- **Label:** amRdata's `resistant_classes`, the genome's resistant classes + joined with `_` (e.g. `MAC_PEN_TMD`), or `Susceptible`. Each combination is + one class of a multi-class model. +- **Filter:** keep combinations with at least `min_n` genomes (default 25) over + the whole dataset, and drop `Intermediate`. +- **Features:** variance > 0 across all kept genomes, test genomes included. + Counts and binary per feature type, binary only for struct (13 matrices). +- **Checks:** none. On S. argenteus it builds a model with a single label (`PEN`). +- **Fitting:** `multinom_reg()` with glmnet, the same tuning grid, chosen by MCC. + Only MCC and nMCC are reported. Importance is one coefficient table per label. +- **Modes:** whole dataset only (`runMLmodels()` stops if MDR is combined with + stratification, LOO or cross-testing). + +## Data check (amRdata main, October 2026) + +`resistant_classes` equals the combination recomputed from the drug-level +phenotypes (Resistant to any drug in a class means that class is resistant) for +every genome: S. flexneri 379/379, S. argenteus 73/73, E. faecium 2,907/2,907. +No genome has more than one value, and none is `Intermediate`. + +With `min_n = 25`: + +| Dataset | Labels kept | Genomes kept | Years / countries among kept genomes | +|---|---|---|---| +| S. flexneri | 4 | 288 of 379 | one year bin (2015–2019), one country (AUS) | +| S. argenteus | 1 (`PEN`) | 65 of 73 | one year bin, one country (THA) | +| E. faecium | 19 | 2,537 of 2,907 | several; many genomes have no year | + +MDR yes/no (2+ classes). "Can't tell" means the genome was tested in fewer +than 2 classes, so it can't be called not MDR: + +| Dataset | MDR | Not MDR | Can't tell | +|---|---|---|---| +| S. flexneri | 300 | 75 | 4 | +| S. argenteus | 5 | 68 | 0 | +| E. faecium | 1,413 | 563 | 931 | + +So only E. faecium exercises the stratified, cross and leave-one-out modes. Its +label mix differs strongly by country (e.g. `GLY_MAC_PEN`: 90 in CAN, 0 in GBR), +which matters for the cross modes below. + +## Design, by stage + +### 1. Profile (`profileORB()`) + +Add two kinds of `unit = "mdr"` rows to `profile$phenotypes`: + +- `target = "mdr"` (**implemented**): each genome gets a status in the new + `profile$mdr` table, with `n_classes_tested` and `n_classes_resistant`: + - `mdr`: Resistant in at least `mdr_min_classes` classes (new + `profileORB()` argument, default 2). Phenotype `Resistant`. + - `not_mdr`: tested in at least `mdr_min_classes` classes but Resistant in + fewer. Phenotype `Susceptible`. + - `undetermined`: tested in fewer than `mdr_min_classes` classes, so MDR + can't be ruled out (931 genomes in E. faecium). Kept as its own label in + `profile$mdr`, but given no `mdr` phenotype, so it never enters the yes/no + model. (Agreed.) + + This is a plain two-phenotype target, so every later stage handles it + unchanged. +- `target = "resistant_classes"`: `phenotype` set to the genome's combination + (or `Susceptible`). + +For the combination target, Build it from +the class phenotypes the profile already computes, sorting class names with +radix order and joining them with `_`. That matches amRdata's string exactly (see the data check) +and keeps a single source of labels. The genome exclusions and contradiction +checks then apply to MDR too. + +`targets` counts stay per target. For the combination target, `n_resistant`/`n_susceptible` don't +apply. Use `n_resistant` = genomes with any resistant class, and +`n_susceptible` = `Susceptible` genomes, so the counts still add up. + +### 2. Eligibility (`eligibleScopes()`) + +The `mdr` yes/no target needs no change: it is judged like a drug or class. +The rest of this section is about the combination target. + +MDR uses the existing scope builders, so it gets scopes in every mode with a +grouping of `none`, `year` or `country`. The drug-grouped modes (`cross_drug`, +`lodo`) only take `unit == "drug"` rows, so they leave MDR out with no change. + +New argument `min_label_genomes = 25` (old `min_n`). Per MDR scope: + +1. Keep the labels with at least `min_label_genomes` **training** genomes, and + drop the scope's genomes (train and test) with any other label. On the whole + dataset this gives exactly old amRml's genome set. Doing it per scope keeps a + label with 3 genomes in one year from breaking that year's folds. +2. Then apply the usual rules, with "phenotype" meaning label: + - `single_phenotype`: fewer than 2 labels kept. + - `too_few_genomes`: fewer than `min_genomes` genomes kept. + - `too_few_for_cv`: the rarest kept label can't fill `n_fold` folds after the + holdout. With 25 genomes per label this always passes. + - `test_too_small`: the test set has fewer than `min_test_genomes` genomes + with kept labels. +3. Record the dropped labels and genomes per scope, so the eligibility figure + can show them. + +The rest of eligibility assumes two phenotypes (`pmin(n_resistant, n_susceptible)`). +MDR needs its own count of the rarest label. The cleanest fix is a per-scope +`n_labels` and `n_rarest` computed for every unit; for drugs, `n_rarest` is the +existing `pmin()`. + +In cross and leave-one-out modes, test genomes whose label never appears in +training can't be predicted correctly, so they are dropped from the test set +(rule 1 already does this) and the number dropped is recorded. (Agreed.) + +The combination label has the same blind spot as the yes/no one: a genome +tested in only one class is labelled from that class alone. Old amRml did the +same, so the port keeps it, but the per-scope record should count these genomes. + +### 3. Matrices (`planMatrices()`, `buildMatrices()`) + +No change expected. The indexed store already keeps each scope's rows with +their `phenotype`, and features are chosen from training genomes only. That is a +deliberate difference from old MDR (it used all genomes), as in #51. + +### 4. Fitting (`fitModels()`) + +The `mdr` yes/no target fits exactly like a drug. For the combination target, +use `multinom_reg()` instead of `logistic_reg()`. +Same engine, grid and MCC selection. Changes needed in `R/fit-models.R`: + +- `.modelData()`: factor levels from the scope's labels, not `.PHENOTYPES`. +- `.stratifiedFolds()`: loop over the labels present, not `.PHENOTYPES`. + Otherwise MDR rows are never assigned a fold. +- `.splitTask()`: already splits per phenotype; no change. +- `.taskMetrics()`: for MDR, MCC and nMCC only, as old amRml did. Columns that + only make sense for Resistant vs Susceptible are `NA`. +- `.taskImportance()`: one coefficient table per label, with a new `label` + column (`NA` for binary models). +- Predictions: `prob_resistant` is `NA` for MDR; add `prob_predicted`. +- Performance: `res_prop` and `res_prop_test` are `NA` for MDR; add `n_labels`. +- Tuning metric set: drop `pr_auc` for MDR, since it needs one probability + column per label. + +The new columns (`label`, `prob_predicted`, `n_labels`) go in the shared +results tables, not in separate MDR tables as old amRml had. (Agreed.) They +should be included in the global rename. + +### 5. Checking against old amRml + +1. S. flexneri, `modes = "vanilla"`: same 288 genomes and 4 labels as the old + MDR matrices, and the same features where the old code's all-genome variance + filter agrees. Compare fits on the old and new matrices with seed 123. +2. S. argenteus: old builds a 1-label MDR model; new records `single_phenotype`. +3. E. faecium: run every mode and check the dropped-label counts by hand for + one cross scope. + +## Order of work + +1. Profile: `mdr` yes/no rows, tests. **Done.** Nothing after the profile needs to change for it, so it runs in every mode. +2. Profile: combination rows, tests. +3. Eligibility: label filter, rules, tests. +4. Fitting: multinomial path, metrics, importance, tests. +5. Comparison with old amRml on S. flexneri. From 376a8b80b7978384f9ed54aa3b6922eeadb0dd8c Mon Sep 17 00:00:00 2001 From: Emily Boyer Date: Fri, 16 Oct 2026 10:50:02 -0600 Subject: [PATCH 02/14] feat: add MDR yes/no target and undetermined status to profileORB() --- R/orb-profile.R | 68 +++++++++++++++++++++++++++---- man/profileORB.Rd | 24 ++++++++--- tests/testthat/test-orb-profile.R | 68 ++++++++++++++++++++++++++++++- 3 files changed, 145 insertions(+), 15 deletions(-) diff --git a/R/orb-profile.R b/R/orb-profile.R index c02a9f8..f94a28e 100644 --- a/R/orb-profile.R +++ b/R/orb-profile.R @@ -99,6 +99,29 @@ ) } +# One row per genome and its MDR status: "mdr" if Resistant in at least `min_classes` +# classes; "not_mdr" if tested in at least that many but Resistant in fewer; otherwise +# "undetermined", since too few classes were tested to rule MDR out. +.mdrStatus <- function(genome_ids, class_phenotypes, min_classes) { + counts <- dplyr::summarise( + dplyr::group_by(class_phenotypes, .data$genome_id), + n_classes_tested = dplyr::n(), + n_classes_resistant = sum(.data$phenotype == "Resistant"), + .groups = "drop" + ) + status <- dplyr::left_join( + tibble::tibble(genome_id = sort(unique(genome_ids), method = "radix")), counts, + by = "genome_id" + ) + status$n_classes_tested[is.na(status$n_classes_tested)] <- 0L + status$n_classes_resistant[is.na(status$n_classes_resistant)] <- 0L + + status$mdr <- "undetermined" + status$mdr[status$n_classes_tested >= min_classes] <- "not_mdr" + status$mdr[status$n_classes_resistant >= min_classes] <- "mdr" + status +} + # One row per genome with its strata, which must have one value per genome. .genomeTable <- function(metadata, call) { columns <- setdiff(names(metadata), c("drug_abbr", "class_abbr", "phenotype")) @@ -158,9 +181,9 @@ #' Profile an amRdata ORB #' -#' Labels each genome Resistant or Susceptible per drug and per drug class, and -#' counts the genomes behind each target, over the whole dataset and within each -#' value of each stratum (year and country). +#' Labels each genome Resistant or Susceptible per drug, per drug class and for +#' multidrug resistance (MDR), and counts the genomes behind each target, over +#' the whole dataset and within each value of each stratum (year and country). #' #' @details #' A genome is excluded, and listed with the reason, if it has no Resistant or @@ -168,6 +191,12 @@ #' any member drug is, and otherwise Susceptible. Repeated rows for a genome and #' drug count once if they agree; if they disagree, profiling stops. #' +#' A genome is MDR (target `mdr`, Resistant) if it is Resistant in at least +#' `mdr_min_classes` drug classes, and not MDR (Susceptible) if it was tested in +#' at least that many classes but is Resistant in fewer. A genome tested in +#' fewer classes can't be ruled out, so it is `"undetermined"` in `$mdr` and has +#' no `mdr` phenotype. +#' #' @section Errors: #' All errors have class `amrml_error`, plus one of: #' @@ -180,13 +209,18 @@ #' * `amrml_stratum_inconsistent`: a genome has two values for a stratum. #' #' @param orb An `amr_orb` from [readORB()]. +#' @param mdr_min_classes Fewest drug classes a genome must be Resistant in to +#' be MDR; at least 2. #' #' @return An `amr_orb_profile` list: #' * `orb`: the ORB profiled. #' * `genomes`: each profiled genome and its strata values. -#' * `phenotypes`: each genome's phenotype per `unit` (drug or drug class) and -#' `target`. +#' * `phenotypes`: each genome's phenotype per `unit` (drug, drug class or +#' mdr) and `target`. #' * `drug_classes`: each profiled drug and its `class`. +#' * `mdr`: each profiled genome's `n_classes_tested`, `n_classes_resistant` +#' and `mdr` status (`"mdr"`, `"not_mdr"` or `"undetermined"`). +#' * `settings`: `mdr_min_classes`. #' * `targets`: genome counts per target and stratum value; `stratum` is NA #' for the whole dataset. #' * `excluded_genomes`: each excluded genome and its `reason` @@ -199,11 +233,12 @@ #' profile$targets #' } #' @export -profileORB <- function(orb) { +profileORB <- function(orb, mdr_min_classes = 2) { call <- rlang::current_env() - # Check the argument is an ORB from readORB(). + # Check the arguments. .checkArgClass(orb, "orb", "amr_orb", "readORB()", call = call) + .checkArgCount(mdr_min_classes, "mdr_min_classes", min = 2, call = call) # Read the metadata and keep its Resistant and Susceptible rows. metadata <- .readMetadata(orb$metadata_parquet, call = call) @@ -246,8 +281,10 @@ profileORB <- function(orb) { )) } - # Each genome's phenotype per drug class. + # Each genome's phenotype per drug class, and its MDR status from those. class_phenotypes <- .classPhenotypes(drug_phenotypes) + mdr <- .mdrStatus(drug_phenotypes$genome_id, class_phenotypes, mdr_min_classes) + modeled_mdr <- mdr[mdr$mdr != "undetermined", , drop = FALSE] phenotypes <- dplyr::bind_rows( tibble::tibble( genome_id = drug_phenotypes$genome_id, @@ -260,6 +297,13 @@ profileORB <- function(orb) { unit = "drug_class", target = class_phenotypes$class_abbr, phenotype = class_phenotypes$phenotype + ), + # Undetermined genomes are left out: the yes/no model needs both phenotypes to mean something. + tibble::tibble( + genome_id = modeled_mdr$genome_id, + unit = "mdr", + target = "mdr", + phenotype = unname(c(mdr = "Resistant", not_mdr = "Susceptible")[modeled_mdr$mdr]) ) ) @@ -284,6 +328,8 @@ profileORB <- function(orb) { genomes = genomes, phenotypes = phenotypes, drug_classes = drug_classes, + mdr = mdr, + settings = list(mdr_min_classes = mdr_min_classes), targets = targets, excluded_genomes = excluded ), @@ -317,5 +363,11 @@ print.amr_orb_profile <- function(x, ...) { cat(counts("drug"), sep = "\n") cat(" per drug class:\n") cat(counts("drug_class"), sep = "\n") + mdr <- table(factor(x$mdr$mdr, c("mdr", "not_mdr", "undetermined"))) + cat( + " MDR (", x$settings$mdr_min_classes, "+ classes): ", mdr[["mdr"]], " MDR, ", + mdr[["not_mdr"]], " not MDR, ", mdr[["undetermined"]], " undetermined\n", + sep = "" + ) invisible(x) } diff --git a/man/profileORB.Rd b/man/profileORB.Rd index e1e53b7..aec3ae4 100644 --- a/man/profileORB.Rd +++ b/man/profileORB.Rd @@ -4,19 +4,25 @@ \alias{profileORB} \title{Profile an amRdata ORB} \usage{ -profileORB(orb) +profileORB(orb, mdr_min_classes = 2) } \arguments{ \item{orb}{An \code{amr_orb} from \code{\link[=readORB]{readORB()}}.} + +\item{mdr_min_classes}{Fewest drug classes a genome must be Resistant in to +be MDR; at least 2.} } \value{ An \code{amr_orb_profile} list: \itemize{ \item \code{orb}: the ORB profiled. \item \code{genomes}: each profiled genome and its strata values. -\item \code{phenotypes}: each genome's phenotype per \code{unit} (drug or drug class) and -\code{target}. +\item \code{phenotypes}: each genome's phenotype per \code{unit} (drug, drug class or +mdr) and \code{target}. \item \code{drug_classes}: each profiled drug and its \code{class}. +\item \code{mdr}: each profiled genome's \code{n_classes_tested}, \code{n_classes_resistant} +and \code{mdr} status (\code{"mdr"}, \code{"not_mdr"} or \code{"undetermined"}). +\item \code{settings}: \code{mdr_min_classes}. \item \code{targets}: genome counts per target and stratum value; \code{stratum} is NA for the whole dataset. \item \code{excluded_genomes}: each excluded genome and its \code{reason} @@ -25,15 +31,21 @@ it is \code{missing_from}. } } \description{ -Labels each genome Resistant or Susceptible per drug and per drug class, and -counts the genomes behind each target, over the whole dataset and within each -value of each stratum (year and country). +Labels each genome Resistant or Susceptible per drug, per drug class and for +multidrug resistance (MDR), and counts the genomes behind each target, over +the whole dataset and within each value of each stratum (year and country). } \details{ A genome is excluded, and listed with the reason, if it has no Resistant or Susceptible row or is missing from a feature table. A class is Resistant if any member drug is, and otherwise Susceptible. Repeated rows for a genome and drug count once if they agree; if they disagree, profiling stops. + +A genome is MDR (target \code{mdr}, Resistant) if it is Resistant in at least +\code{mdr_min_classes} drug classes, and not MDR (Susceptible) if it was tested in +at least that many classes but is Resistant in fewer. A genome tested in +fewer classes can't be ruled out, so it is \code{"undetermined"} in \verb{$mdr} and has +no \code{mdr} phenotype. } \section{Errors}{ diff --git a/tests/testthat/test-orb-profile.R b/tests/testthat/test-orb-profile.R index c7edc32..1de2c70 100644 --- a/tests/testthat/test-orb-profile.R +++ b/tests/testthat/test-orb-profile.R @@ -217,6 +217,68 @@ test_that("strata count observed values; genomes without one count in the whole expect_equal(fxTarget(profile, "drug", "DRA")$n_genomes, 4) }) +# MDR + +# g1: Resistant in 3 classes. g2: Resistant in 2. g3: tested in 3, Resistant in 1. +# g4: tested in 1 class only. g5: tested in 2 drugs of one class only. +fxMdrMetadata <- function() { + fxMetadataRows( + genome = c("g1", "g1", "g1", "g2", "g2", "g2", "g3", "g3", "g3", "g4", "g5", "g5"), + drug = c("DRA", "DRB", "DRC", "DRA", "DRB", "DRC", "DRA", "DRB", "DRC", "DRA", "DRA", "DRD"), + phenotype = c( + "Resistant", "Resistant", "Resistant", + "Resistant", "Resistant", "Susceptible", + "Resistant", "Susceptible", "Susceptible", + "Resistant", + "Resistant", "Resistant" + ), + class = c( + "CLX", "CLY", "CLZ", "CLX", "CLY", "CLZ", "CLX", "CLY", "CLZ", "CLX", "CLX", "CLX" + ) + ) +} + +test_that("MDR status counts classes tested and resistant per genome", { + profile <- fxProfile(fxMdrMetadata(), parquets = fxDefaultParquets()) + + expect_equal(profile$mdr, tibble::tibble( + genome_id = c("g1", "g2", "g3", "g4", "g5"), + n_classes_tested = c(3L, 3L, 3L, 1L, 1L), + n_classes_resistant = c(3L, 2L, 1L, 1L, 1L), + mdr = c("mdr", "mdr", "not_mdr", "undetermined", "undetermined") + )) +}) + +test_that("undetermined genomes get no mdr phenotype; the rest are Resistant or Susceptible", { + profile <- fxProfile(fxMdrMetadata(), parquets = fxDefaultParquets()) + mdr <- profile$phenotypes[profile$phenotypes$unit == "mdr", ] + + expect_equal(mdr, tibble::tibble( + genome_id = c("g1", "g2", "g3"), + unit = "mdr", + target = "mdr", + phenotype = c("Resistant", "Resistant", "Susceptible") + )) + expect_equal(fxTarget(profile, "mdr", "mdr")$n_resistant, 2) + expect_equal(fxTarget(profile, "mdr", "mdr")$n_susceptible, 1) + expect_output(print(profile), "2 MDR, 1 not MDR, 2 undetermined") +}) + +test_that("a higher mdr_min_classes needs more classes tested and resistant", { + orb <- readORB(fxOrb(metadata = fxMdrMetadata())) + profile <- profileORB(orb, mdr_min_classes = 3) + + expect_equal(profile$mdr$mdr, c("mdr", "not_mdr", "not_mdr", "undetermined", "undetermined")) + expect_equal(profile$settings$mdr_min_classes, 3) +}) + +test_that("mdr_min_classes must be a whole number of at least 2", { + orb <- readORB(fxOrb()) + + expect_error(profileORB(orb, mdr_min_classes = 1), class = "amrml_invalid_argument") + expect_error(profileORB(orb, mdr_min_classes = 2.5), class = "amrml_invalid_argument") +}) + # The profile object test_that("the profile has its parts and prints a summary", { @@ -225,12 +287,16 @@ test_that("the profile has its parts and prints a summary", { expect_s3_class(profile, "amr_orb_profile") expect_named( profile, - c("orb", "genomes", "phenotypes", "drug_classes", "targets", "excluded_genomes") + c( + "orb", "genomes", "phenotypes", "drug_classes", "mdr", "settings", "targets", + "excluded_genomes" + ) ) expect_output(print(profile), "5 profiled, 0 excluded") expect_output(print(profile), "country \\(values: 1, missing: 0\\)") expect_output(print(profile), "DRA:5/0, DRB:0/5") expect_output(print(profile), "CLX:5/0") + expect_output(print(profile), "MDR \\(2\\+ classes\\): 0 MDR, 0 not MDR, 5 undetermined") }) test_that("profiling does not write to the ORB", { From 2b0f6782d74e7a4cb7a3db5e297561f5fa325a12 Mon Sep 17 00:00:00 2001 From: Emily Boyer Date: Fri, 16 Oct 2026 11:31:42 -0600 Subject: [PATCH 03/14] feat: add MDR combination target to profile and eligibility Each genome tested in a class gets its combination of resistant classes as a label (target resistant_classes). Eligibility keeps labels with at least min_label_genomes (default 25) training genomes per scope and drops genomes, training or test, with other labels. --- R/eligibility.R | 80 +++++++++++++++++++++----- R/orb-profile.R | 41 +++++++++++++- R/workflow.R | 6 +- man/eligibleScopes.Rd | 14 ++++- man/profileORB.Rd | 6 ++ man/runModelingWorkflow.Rd | 9 ++- tests/testthat/test-eligibility.R | 94 +++++++++++++++++++++++++++++-- tests/testthat/test-orb-profile.R | 23 +++++++- 8 files changed, 245 insertions(+), 28 deletions(-) diff --git a/R/eligibility.R b/R/eligibility.R index 07ddafc..6073e4a 100644 --- a/R/eligibility.R +++ b/R/eligibility.R @@ -122,6 +122,25 @@ .asMembers(train, test) } +# Combination scopes keep the labels with at least `min_label_genomes` training genomes, and +# drop every genome, training or test, with another label: a test genome whose label the +# model never trained on can't be predicted. Other scopes keep every genome. +.keepLabels <- function(members, min_label_genomes) { + members$kept <- TRUE + combination <- members$unit == "mdr" & members$target == .MDR_COMBINATION + if (!any(combination)) { + return(members) + } + + keys <- c(.SCOPE_KEYS, "phenotype") + train <- members[combination & members$role == "train", keys] + n <- dplyr::count(train, dplyr::across(dplyr::all_of(keys))) + labels <- n[n$n >= min_label_genomes, keys] + key <- function(x) do.call(paste, c(unname(as.list(x[keys])), sep = "\r")) + members$kept[combination] <- key(members[combination, ]) %in% key(labels) + members +} + # Each target tested in fewer of the mode's groups than it needs, with its genome counts in # them. A target whose genomes all lack a group has none, and is listed too. .targetsInTooFewGroups <- function(phenotypes, grouped, loo) { @@ -129,7 +148,7 @@ dplyr::group_by(grouped, .data$unit, .data$target), n_groups = dplyr::n_distinct(.data$group), n_genomes = dplyr::n(), - n_resistant = sum(.data$phenotype == "Resistant"), + n_resistant = sum(.data$phenotype != "Susceptible"), n_susceptible = sum(.data$phenotype == "Susceptible"), .groups = "drop" ) @@ -174,17 +193,21 @@ list(members = members, too_few = too_few) } -# Per scope, the training and test genome counts. +# Per scope, the training and test genome counts of kept genomes, with how many labels they +# have, the size of the rarest, and how many genomes the label filter dropped. .scopeCounts <- function(members, has_test) { counts <- dplyr::summarise( dplyr::group_by(members, dplyr::across(dplyr::all_of(c(.SCOPE_KEYS, "role")))), - n_genomes = dplyr::n(), - n_resistant = sum(.data$phenotype == "Resistant"), - n_susceptible = sum(.data$phenotype == "Susceptible"), + n_genomes = sum(.data$kept), + n_resistant = sum(.data$kept & .data$phenotype != "Susceptible"), + n_susceptible = sum(.data$kept & .data$phenotype == "Susceptible"), + n_labels = dplyr::n_distinct(.data$phenotype[.data$kept]), + n_rarest = if (any(.data$kept)) min(table(.data$phenotype[.data$kept])) else 0L, + n_dropped = sum(!.data$kept), .groups = "drop" ) - sides <- c("n_genomes", "n_resistant", "n_susceptible") + sides <- c("n_genomes", "n_resistant", "n_susceptible", "n_labels", "n_rarest", "n_dropped") train <- counts[counts$role == "train", c(.SCOPE_KEYS, sides)] test <- counts[counts$role == "test", c(.SCOPE_KEYS, sides)] names(test) <- c(.SCOPE_KEYS, paste0("test_", sides)) @@ -224,7 +247,13 @@ # A class with data for only one of its drugs is that drug's model again. duplicate <- scopes$unit == "drug_class" & scopes$n_drugs %in% 1 - rarer <- pmin(scopes$n_resistant, scopes$n_susceptible) + # A combination scope has a label per combination, so its "rarer phenotype" is its + # rarest label, and it needs two labels rather than both phenotypes. + combination <- scopes$unit == "mdr" & scopes$target == .MDR_COMBINATION + rarer <- ifelse( + combination, scopes$n_rarest, pmin(scopes$n_resistant, scopes$n_susceptible) + ) + single <- ifelse(combination, scopes$n_labels < 2, rarer == 0) enough_for_cv <- if (has_test) { rarer >= settings$n_fold } else { @@ -232,13 +261,14 @@ } test_ok <- !has_test | ( scopes$test_n_genomes >= settings$min_test_genomes & - pmin(scopes$test_n_resistant, scopes$test_n_susceptible) >= settings$min_test_minority + (combination | + pmin(scopes$test_n_resistant, scopes$test_n_susceptible) >= settings$min_test_minority) ) # Check the rules in order and keep each scope's first failure. failed <- list( duplicate_class = duplicate, - single_phenotype = rarer == 0, + single_phenotype = single, too_few_genomes = scopes$n_genomes < settings$min_genomes, too_few_for_cv = !enough_for_cv, test_too_small = !test_ok @@ -263,8 +293,9 @@ # in too few groups get one decision each. .judgeScopes <- function(built, mode, profile, settings) { counts <- c( - "n_genomes", "n_resistant", "n_susceptible", - "test_n_genomes", "test_n_resistant", "test_n_susceptible" + "n_genomes", "n_resistant", "n_susceptible", "n_labels", "n_dropped", + "test_n_genomes", "test_n_resistant", "test_n_susceptible", "test_n_labels", + "test_n_dropped" ) has_test <- mode$LOO || mode$cross_test scopes <- .scopeCounts(built$members, has_test) @@ -275,6 +306,8 @@ too_few <- built$too_few if (NROW(too_few)) { too_few[c("train_group", "test_group")] <- NA_character_ + too_few$n_labels <- NA_integer_ + too_few$n_dropped <- NA_integer_ too_few$rule_id <- "too_few_groups_tested" decisions <- dplyr::bind_rows(decisions, too_few[names(too_few) != "n_groups"]) } @@ -296,7 +329,12 @@ #' built or fitted. #' #' @details -#' Counts are of genomes, each with one phenotype. A drug or class needs 2 of a +#' Counts are of genomes, each with one phenotype. In the MDR combination target +#' (`resistant_classes`), each combination is a label: a scope keeps the labels +#' with at least `min_label_genomes` training genomes and drops the genomes, +#' training or test, with any other (`n_dropped`, `test_n_dropped`). It then +#' needs 2 labels, and its rarest label must fill the folds; its test set needs +#' only `min_test_genomes`. A drug or class needs 2 of a #' mode's groups (3 for leave-one-out) to get scopes, and a stratified one needs #' 2 eligible groups. A scope fails, in order, if it's a class with data for only #' one of its drugs, has one phenotype, has fewer than `min_genomes` genomes, @@ -314,12 +352,15 @@ #' @param n_fold Number of cross-validation folds. #' @param min_test_genomes Fewest genomes a separate test set needs. #' @param min_test_minority Fewest of the rarer phenotype a test set needs. +#' @param min_label_genomes Fewest training genomes an MDR combination needs to +#' stay a label in a scope (old amRml's `min_n`). #' #' @return An `amr_eligibility` list: #' * `profile` and `settings` (the thresholds, and the `holdout` share). #' * `modes`: per mode, `feasible`, `rule_id`, `n_groups`, and how many #' scopes it has and how many are eligible (`n_scopes`, `n_eligible`). #' * `scopes`: per scope, training and test counts (`NA` without a test set), +#' the number of labels kept and genomes dropped by the label filter, #' `eligible` and `rule_id`. #' * `members`: each eligible scope's genomes, by `role` (train or test). #' @@ -334,7 +375,8 @@ eligibleScopes <- function(profile, min_genomes = 40, n_fold = 5, min_test_genomes = 5, - min_test_minority = 2) { + min_test_minority = 2, + min_label_genomes = 25) { call <- rlang::current_env() # Check the arguments. @@ -343,6 +385,7 @@ eligibleScopes <- function(profile, .checkArgCount(n_fold, "n_fold", min = 2, call = call) .checkArgCount(min_test_genomes, "min_test_genomes", call = call) .checkArgCount(min_test_minority, "min_test_minority", call = call) + .checkArgCount(min_label_genomes, "min_label_genomes", call = call) modes <- modes %||% .MODES$mode_id .checkArgModes(modes, call = call) @@ -351,6 +394,7 @@ eligibleScopes <- function(profile, n_fold = n_fold, min_test_genomes = min_test_genomes, min_test_minority = min_test_minority, + min_label_genomes = min_label_genomes, holdout = .CV_HOLDOUT ) @@ -364,13 +408,17 @@ eligibleScopes <- function(profile, judged <- lapply(seq_len(nrow(runnable)), function(i) { mode <- runnable[i, , drop = FALSE] built <- .buildScopes(mode, profile) + built$members <- .keepLabels(built$members, settings$min_label_genomes) decisions <- .judgeScopes(built, mode, profile, settings) eligible <- decisions[decisions$eligible, .SCOPE_KEYS] list( decisions = decisions, members = tibble::tibble( mode_id = mode$mode_id, - dplyr::semi_join(built$members, eligible, by = .SCOPE_KEYS) + dplyr::semi_join( + built$members[built$members$kept, names(built$members) != "kept"], + eligible, by = .SCOPE_KEYS + ) ) ) }) @@ -383,8 +431,10 @@ eligibleScopes <- function(profile, mode_id = character(), unit = character(), target = character(), train_group = character(), test_group = character(), n_genomes = integer(), n_resistant = integer(), n_susceptible = integer(), + n_labels = integer(), n_dropped = integer(), test_n_genomes = integer(), test_n_resistant = integer(), - test_n_susceptible = integer(), eligible = logical(), rule_id = character() + test_n_susceptible = integer(), test_n_labels = integer(), + test_n_dropped = integer(), eligible = logical(), rule_id = character() ) members <- tibble::tibble( mode_id = character(), unit = character(), target = character(), diff --git a/R/orb-profile.R b/R/orb-profile.R index f94a28e..268ad37 100644 --- a/R/orb-profile.R +++ b/R/orb-profile.R @@ -1,5 +1,9 @@ # Profiling an amRdata ORB: each genome's phenotypes, and the counts per target. +# The MDR target labelling each genome with its combination of resistant classes, as +# amRdata's `resistant_classes` does. +.MDR_COMBINATION <- "resistant_classes" + # amRdata's metadata.parquet columns that profiling needs, named as they are used here. .METADATA_COLUMNS <- c( genome_id = "genome.genome_id", @@ -122,6 +126,21 @@ status } +# One row per genome tested in any class: its resistant classes in radix order, joined +# with "_" (e.g. "MAC_PEN_TMD"), or "Susceptible" if it has none. This matches amRdata's +# `resistant_classes`. +.classCombinations <- function(class_phenotypes) { + dplyr::summarise( + dplyr::group_by(class_phenotypes, .data$genome_id), + phenotype = if (any(.data$phenotype == "Resistant")) { + paste(sort(.data$class_abbr[.data$phenotype == "Resistant"], method = "radix"), collapse = "_") + } else { + "Susceptible" + }, + .groups = "drop" + ) +} + # One row per genome with its strata, which must have one value per genome. .genomeTable <- function(metadata, call) { columns <- setdiff(names(metadata), c("drug_abbr", "class_abbr", "phenotype")) @@ -168,7 +187,8 @@ .data$unit, .data$target, .data$stratum, .data$stratum_value ), n_genomes = dplyr::n(), - n_resistant = sum(.data$phenotype == "Resistant"), + # Anything not Susceptible is resistant, so a combination counts as resistant too. + n_resistant = sum(.data$phenotype != "Susceptible"), n_susceptible = sum(.data$phenotype == "Susceptible"), .groups = "drop" ) @@ -197,6 +217,12 @@ #' fewer classes can't be ruled out, so it is `"undetermined"` in `$mdr` and has #' no `mdr` phenotype. #' +#' Target `resistant_classes` labels each genome with its combination of +#' resistant classes, joined with `_` in radix order (e.g. `"MAC_PEN_TMD"`), or +#' `"Susceptible"`, as amRdata's column of that name does. Every genome tested in +#' a class gets one, as in old amRml. In `targets`, `n_resistant` counts its +#' genomes with any resistant class. +#' #' @section Errors: #' All errors have class `amrml_error`, plus one of: #' @@ -285,6 +311,7 @@ profileORB <- function(orb, mdr_min_classes = 2) { class_phenotypes <- .classPhenotypes(drug_phenotypes) mdr <- .mdrStatus(drug_phenotypes$genome_id, class_phenotypes, mdr_min_classes) modeled_mdr <- mdr[mdr$mdr != "undetermined", , drop = FALSE] + combinations <- .classCombinations(class_phenotypes) phenotypes <- dplyr::bind_rows( tibble::tibble( genome_id = drug_phenotypes$genome_id, @@ -304,6 +331,12 @@ profileORB <- function(orb, mdr_min_classes = 2) { unit = "mdr", target = "mdr", phenotype = unname(c(mdr = "Resistant", not_mdr = "Susceptible")[modeled_mdr$mdr]) + ), + tibble::tibble( + genome_id = combinations$genome_id, + unit = "mdr", + target = .MDR_COMBINATION, + phenotype = combinations$phenotype ) ) @@ -369,5 +402,11 @@ print.amr_orb_profile <- function(x, ...) { mdr[["not_mdr"]], " not MDR, ", mdr[["undetermined"]], " undetermined\n", sep = "" ) + combinations <- x$phenotypes$phenotype[x$phenotypes$target == .MDR_COMBINATION] + cat( + " MDR class combinations: ", dplyr::n_distinct(combinations), " across ", + length(combinations), " genomes\n", + sep = "" + ) invisible(x) } diff --git a/R/workflow.R b/R/workflow.R index 3fdc063..47397c9 100644 --- a/R/workflow.R +++ b/R/workflow.R @@ -16,7 +16,7 @@ #' `/amRml//`. #' @param modes Mode IDs from [modelingModes()], or `NULL` for all of them. #' @param seeds,labels Passed to [expandTasks()]. -#' @param min_genomes,n_fold,min_test_genomes,min_test_minority Passed to +#' @param min_genomes,n_fold,min_test_genomes,min_test_minority,min_label_genomes Passed to #' [eligibleScopes()]. #' @param overwrite Replace an earlier run's matrices and results in `out_dir`; #' if `FALSE`, an `out_dir` holding earlier output is refused. The default run @@ -45,6 +45,7 @@ runModelingWorkflow <- function(orb_dir, n_fold = 5, min_test_genomes = 5, min_test_minority = 2, + min_label_genomes = 25, overwrite = FALSE, BPPARAM = BiocParallel::SnowParam(workers = 8L, tasks = 128L)) { call <- rlang::current_env() @@ -59,7 +60,8 @@ runModelingWorkflow <- function(orb_dir, eligibility <- eligibleScopes( profileORB(orb), modes = modes, min_genomes = min_genomes, n_fold = n_fold, - min_test_genomes = min_test_genomes, min_test_minority = min_test_minority + min_test_genomes = min_test_genomes, min_test_minority = min_test_minority, + min_label_genomes = min_label_genomes ) built <- buildMatrices(planMatrices(eligibility), out_dir, overwrite = overwrite) arrow::write_parquet(eligibility$scopes, file.path(out_dir, "eligibility.parquet")) diff --git a/man/eligibleScopes.Rd b/man/eligibleScopes.Rd index ecc53c8..8989873 100644 --- a/man/eligibleScopes.Rd +++ b/man/eligibleScopes.Rd @@ -10,7 +10,8 @@ eligibleScopes( min_genomes = 40, n_fold = 5, min_test_genomes = 5, - min_test_minority = 2 + min_test_minority = 2, + min_label_genomes = 25 ) } \arguments{ @@ -25,6 +26,9 @@ eligibleScopes( \item{min_test_genomes}{Fewest genomes a separate test set needs.} \item{min_test_minority}{Fewest of the rarer phenotype a test set needs.} + +\item{min_label_genomes}{Fewest training genomes an MDR combination needs to +stay a label in a scope (old amRml's \code{min_n}).} } \value{ An \code{amr_eligibility} list: @@ -33,6 +37,7 @@ An \code{amr_eligibility} list: \item \code{modes}: per mode, \code{feasible}, \code{rule_id}, \code{n_groups}, and how many scopes it has and how many are eligible (\code{n_scopes}, \code{n_eligible}). \item \code{scopes}: per scope, training and test counts (\code{NA} without a test set), +the number of labels kept and genomes dropped by the label filter, \code{eligible} and \code{rule_id}. \item \code{members}: each eligible scope's genomes, by \code{role} (train or test). } @@ -44,7 +49,12 @@ can, records whether it has enough data to model, or why not. Nothing is built or fitted. } \details{ -Counts are of genomes, each with one phenotype. A drug or class needs 2 of a +Counts are of genomes, each with one phenotype. In the MDR combination target +(\code{resistant_classes}), each combination is a label: a scope keeps the labels +with at least \code{min_label_genomes} training genomes and drops the genomes, +training or test, with any other (\code{n_dropped}, \code{test_n_dropped}). It then +needs 2 labels, and its rarest label must fill the folds; its test set needs +only \code{min_test_genomes}. A drug or class needs 2 of a mode's groups (3 for leave-one-out) to get scopes, and a stratified one needs 2 eligible groups. A scope fails, in order, if it's a class with data for only one of its drugs, has one phenotype, has fewer than \code{min_genomes} genomes, diff --git a/man/profileORB.Rd b/man/profileORB.Rd index aec3ae4..4421b65 100644 --- a/man/profileORB.Rd +++ b/man/profileORB.Rd @@ -46,6 +46,12 @@ A genome is MDR (target \code{mdr}, Resistant) if it is Resistant in at least at least that many classes but is Resistant in fewer. A genome tested in fewer classes can't be ruled out, so it is \code{"undetermined"} in \verb{$mdr} and has no \code{mdr} phenotype. + +Target \code{resistant_classes} labels each genome with its combination of +resistant classes, joined with \verb{_} in radix order (e.g. \code{"MAC_PEN_TMD"}), or +\code{"Susceptible"}, as amRdata's column of that name does. Every genome tested in +a class gets one, as in old amRml. In \code{targets}, \code{n_resistant} counts its +genomes with any resistant class. } \section{Errors}{ diff --git a/man/runModelingWorkflow.Rd b/man/runModelingWorkflow.Rd index a416e05..eafb92e 100644 --- a/man/runModelingWorkflow.Rd +++ b/man/runModelingWorkflow.Rd @@ -14,7 +14,9 @@ runModelingWorkflow( n_fold = 5, min_test_genomes = 5, min_test_minority = 2, - overwrite = FALSE + min_label_genomes = 25, + overwrite = FALSE, + BPPARAM = BiocParallel::SnowParam(workers = 8L, tasks = 128L) ) } \arguments{ @@ -27,12 +29,15 @@ runModelingWorkflow( \item{seeds, labels}{Passed to \code{\link[=expandTasks]{expandTasks()}}.} -\item{min_genomes, n_fold, min_test_genomes, min_test_minority}{Passed to +\item{min_genomes, n_fold, min_test_genomes, min_test_minority, min_label_genomes}{Passed to \code{\link[=eligibleScopes]{eligibleScopes()}}.} \item{overwrite}{Replace an earlier run's matrices and results in \code{out_dir}; if \code{FALSE}, an \code{out_dir} holding earlier output is refused. The default run folder is always new.} + +\item{BPPARAM}{BiocParallel backend passed to \code{\link[=fitModels]{fitModels()}}. Defaults to a +eight-worker cluster with 128 chunks of work to chew through.} } \value{ The \code{amr_fit} from \code{\link[=fitModels]{fitModels()}}; every earlier stage's result can be diff --git a/tests/testthat/test-eligibility.R b/tests/testthat/test-eligibility.R index 5702672..7420c82 100644 --- a/tests/testthat/test-eligibility.R +++ b/tests/testthat/test-eligibility.R @@ -51,7 +51,8 @@ test_that("grouped modes need 2 groups, and leave-one-out needs 3", { expect_equal(modes$rule_id, c(NA, NA, "too_few_groups")) expect_equal(modes$n_groups, c(2, 2, 2)) # Feasible modes can still have no eligible scope: Y1 is all Resistant, Y2 all Susceptible. - expect_equal(modes$n_scopes, c(4, 4, NA)) + # DRA, its class and the MDR combination target each get a scope per year. + expect_equal(modes$n_scopes, c(6, 6, NA)) expect_equal(modes$n_eligible, c(0, 0, NA)) one_drug <- fxEligibility(metadata, modes = "cross_drug")$modes @@ -133,11 +134,16 @@ test_that("vanilla scopes are the profile's whole-dataset counts", { profile <- profileORB(readORB(fxOrb())) scopes <- fxScopes(eligibleScopes(profile, modes = "vanilla"), "vanilla") whole <- profile$targets[is.na(profile$targets$stratum), ] - expect_equal(scopes$target, whole$target) - expect_equal(scopes$n_genomes, whole$n_genomes) - expect_equal(scopes$n_resistant, whole$n_resistant) expect_true(all(is.na(scopes$test_n_genomes))) + + # The MDR combination target counts only the genomes its label filter keeps: here none, + # since its one label has 5 genomes, under min_label_genomes. + combination <- scopes$target == "resistant_classes" + expect_equal(scopes$n_genomes[!combination], whole$n_genomes[!combination]) + expect_equal(scopes$n_resistant[!combination], whole$n_resistant[!combination]) + expect_equal(scopes$n_genomes[combination], 0) + expect_equal(scopes$n_dropped[combination], 5) }) test_that("stratified scopes are per value, leaving out genomes without one", { @@ -353,6 +359,84 @@ test_that("a stratified drug needs at least two groups with enough data", { expect_equal(drug$rule_id, c("too_few_eligible_groups", "too_few_genomes")) }) +# The MDR combination target + +# Genomes tested on DRA (class CLX) and DRB (class CLY), numbered from `from`: `n_x` resistant +# to CLX only, `n_xy` to both, and `n_s` to neither, all in `year_bin`. +fxCombinations <- function(n_x, n_xy, n_s, from = 1, year_bin = "Y1") { + n <- n_x + n_xy + n_s + genomes <- sprintf("g%03d", from - 1 + seq_len(n)) + label <- rep(c("x", "xy", "s"), c(n_x, n_xy, n_s)) + fxMetadataRows( + genome = rep(genomes, each = 2), + drug = rep(c("DRA", "DRB"), n), + class = rep(c("CLX", "CLY"), n), + phenotype = ifelse( + rep(label, each = 2) == "xy" | (rep(label, each = 2) == "x" & rep(c(TRUE, FALSE), n)), + "Resistant", "Susceptible" + ), + year_bin = year_bin + ) +} + +fxCombinationScopes <- function(eligibility, mode_id) { + scopes <- fxScopes(eligibility, mode_id) + scopes[scopes$target == "resistant_classes", ] +} + +test_that("a combination scope drops labels with too few training genomes, and their genomes", { + eligibility <- fxEligibility( + fxCombinations(4, 3, 2), + modes = "vanilla", min_genomes = 4, n_fold = 2, min_label_genomes = 3 + ) + scope <- fxCombinationScopes(eligibility, "vanilla") + + expect_equal(scope$n_genomes, 7) + expect_equal(scope$n_labels, 2) + expect_equal(scope$n_dropped, 2) + expect_true(scope$eligible) + kept <- eligibility$members[eligibility$members$target == "resistant_classes", ] + expect_equal(sort(unique(kept$phenotype)), c("CLX", "CLX_CLY")) + expect_false(any(kept$genome_id %in% c("g008", "g009"))) +}) + +test_that("a combination scope with one label left is single_phenotype", { + eligibility <- fxEligibility( + fxCombinations(4, 1, 1), + modes = "vanilla", min_genomes = 2, n_fold = 2, min_label_genomes = 3 + ) + scope <- fxCombinationScopes(eligibility, "vanilla") + + expect_equal(scope$n_labels, 1) + expect_equal(scope$rule_id, "single_phenotype") + expect_false(any(eligibility$members$target == "resistant_classes")) +}) + +test_that("a cross test set drops genomes whose label the training set doesn't keep", { + metadata <- rbind( + fxCombinations(3, 3, 0, from = 1, year_bin = "Y1"), + fxCombinations(2, 0, 2, from = 7, year_bin = "Y2") + ) + eligibility <- fxEligibility( + metadata, + modes = "cross_year", min_genomes = 4, n_fold = 2, min_test_genomes = 2, + min_label_genomes = 3 + ) + scopes <- fxCombinationScopes(eligibility, "cross_year") + y1_to_y2 <- scopes[scopes$train_group == "Y1", ] + + # Y2's two Susceptible genomes have a label Y1 never trained on. + expect_equal(y1_to_y2$test_n_genomes, 2) + expect_equal(y1_to_y2$test_n_dropped, 2) + expect_true(y1_to_y2$eligible) + expect_equal( + fxMemberIds(eligibility, "cross_year", "resistant_classes", "Y1", "Y2", "test"), + c("g007", "g008") + ) + # Y2 has no label with 3 training genomes, so nothing trains on it. + expect_equal(scopes$rule_id[scopes$train_group == "Y2"], "single_phenotype") +}) + # The result test_that("the result keeps the profile and settings, and every rule is an error code", { @@ -414,7 +498,7 @@ test_that("the print shows each mode's result", { eligibility <- fxEligibility(modes = c("vanilla", "stratified_year", "lodo")) expect_output(print(eligibility), "40 training genomes, 5 folds") - expect_output(print(eligibility), "vanilla +: 0 of 3 eligible") + expect_output(print(eligibility), "vanilla +: 0 of 4 eligible") expect_output(print(eligibility), "stratified_year +: not run \\(too_few_groups\\)") expect_output(print(eligibility), "lodo +: not run \\(mode_not_supported\\)") }) diff --git a/tests/testthat/test-orb-profile.R b/tests/testthat/test-orb-profile.R index 1de2c70..4dda3f5 100644 --- a/tests/testthat/test-orb-profile.R +++ b/tests/testthat/test-orb-profile.R @@ -251,7 +251,7 @@ test_that("MDR status counts classes tested and resistant per genome", { test_that("undetermined genomes get no mdr phenotype; the rest are Resistant or Susceptible", { profile <- fxProfile(fxMdrMetadata(), parquets = fxDefaultParquets()) - mdr <- profile$phenotypes[profile$phenotypes$unit == "mdr", ] + mdr <- profile$phenotypes[profile$phenotypes$target == "mdr", ] expect_equal(mdr, tibble::tibble( genome_id = c("g1", "g2", "g3"), @@ -272,6 +272,27 @@ test_that("a higher mdr_min_classes needs more classes tested and resistant", { expect_equal(profile$settings$mdr_min_classes, 3) }) +test_that("each genome tested in a class gets its combination of resistant classes", { + metadata <- rbind( + fxMdrMetadata(), + fxMetadataRows(genome = "g6", drug = c("DRA", "DRB"), phenotype = "Susceptible", + class = c("CLX", "CLY")) + ) + profile <- fxProfile(metadata, parquets = fxDefaultParquets(c(fxGenomes(), "g6"))) + combinations <- profile$phenotypes[profile$phenotypes$target == "resistant_classes", ] + + expect_equal(combinations, tibble::tibble( + genome_id = c("g1", "g2", "g3", "g4", "g5", "g6"), + unit = "mdr", + target = "resistant_classes", + phenotype = c("CLX_CLY_CLZ", "CLX_CLY", "CLX", "CLX", "CLX", "Susceptible") + )) + # Every combination counts as resistant. + expect_equal(fxTarget(profile, "mdr", "resistant_classes")$n_resistant, 5) + expect_equal(fxTarget(profile, "mdr", "resistant_classes")$n_susceptible, 1) + expect_output(print(profile), "MDR class combinations: 4 across 6 genomes") +}) + test_that("mdr_min_classes must be a whole number of at least 2", { orb <- readORB(fxOrb()) From 4be5691970d6e662003231a02851ca7a0a84af95 Mon Sep 17 00:00:00 2001 From: Emily Boyer Date: Fri, 16 Oct 2026 11:31:42 -0600 Subject: [PATCH 04/14] feat: fit MDR combination target with multinomial regression Tasks with labels other than Resistant/Susceptible use multinom_reg (glmnet), report MCC and nMCC only, and keep coefficients per label. Adds label to importance, n_labels to performance and prob_predicted to predictions. --- R/fit-models.R | 123 +++++++++++++++++++++---------- man/fitModels.Rd | 23 ++++-- tests/testthat/test-fit-models.R | 81 ++++++++++++++++++++ 3 files changed, 183 insertions(+), 44 deletions(-) create mode 100644 tests/testthat/test-fit-models.R diff --git a/R/fit-models.R b/R/fit-models.R index 41ebebb..eb98cc0 100644 --- a/R/fit-models.R +++ b/R/fit-models.R @@ -1,5 +1,6 @@ # Fitting a model for each task, with old amRml's method: glmnet logistic regression, -# tuned by stratified cross-validation and chosen by MCC. +# tuned by stratified cross-validation and chosen by MCC. A task whose labels aren't +# Resistant and Susceptible (the MDR combination target) gets multinomial regression. # The phenotypes, with Resistant first so it is the event yardstick scores. .PHENOTYPES <- c("Resistant", "Susceptible") @@ -56,7 +57,7 @@ # without pooling rare phenotypes (`pool = 0.1`), which can leave a fold without them. .stratifiedFolds <- function(phenotype, n_fold) { fold <- integer(length(phenotype)) - for (level in intersect(.PHENOTYPES, phenotype)) { + for (level in intersect(levels(phenotype), phenotype)) { i <- which(phenotype == level) fold[i] <- sample(rep_len(seq_len(n_fold), length(i))) } @@ -64,8 +65,8 @@ } # The training data as a data frame: the label, then features renamed x1, x2, ..., so any -# feature ID works in a formula. -.modelData <- function(values, phenotype, sparse = FALSE) { +# feature ID works in a formula. `labels` are the label's levels, in order. +.modelData <- function(values, phenotype, sparse = FALSE, labels = .PHENOTYPES) { if (sparse) { colnames(values) <- paste0("x", seq_len(ncol(values))) data <- sparsevctrs::coerce_to_sparse_tibble(values) @@ -73,7 +74,7 @@ data <- as.data.frame(as.matrix(values)) } names(data) <- paste0("x", seq_len(ncol(data))) - data$phenotype <- factor(phenotype, levels = .PHENOTYPES) + data$phenotype <- factor(phenotype, levels = labels) data } @@ -99,10 +100,15 @@ } } -# Ported from old amRml's buildLRModel(). -.modelSpec <- function() { - parsnip::logistic_reg(penalty = tune::tune(), mixture = tune::tune()) |> - parsnip::set_engine("glmnet") +# Ported from old amRml's buildLRModel(), which also switched to multinomial regression +# for more than two labels. +.modelSpec <- function(multi = FALSE) { + spec <- if (multi) { + parsnip::multinom_reg(penalty = tune::tune(), mixture = tune::tune()) + } else { + parsnip::logistic_reg(penalty = tune::tune(), mixture = tune::tune()) + } + parsnip::set_engine(spec, "glmnet") } # Ported from old amRml's buildTuningGrid(): 10 penalties x 6 mixtures. @@ -117,10 +123,10 @@ # Ported from old amRml's tuneGrid(), selectBestModel() and fitBestModel(): tune over # stratified folds, choose by MCC, and fit on all training rows. -.fitModel <- function(data, n_fold, sparse = FALSE) { +.fitModel <- function(data, n_fold, sparse = FALSE, multi = FALSE) { workflow <- workflows::workflow() |> workflows::add_recipe(.buildRecipe(data, sparse = sparse)) |> - workflows::add_model(.modelSpec()) + workflows::add_model(.modelSpec(multi)) # After the recipe, whose step IDs use random numbers, so the folds match rsample's. fold <- .stratifiedFolds(data$phenotype, n_fold) @@ -137,10 +143,15 @@ workflow, resamples = resamples, grid = .tuningGrid(), - metrics = yardstick::metric_set( - yardstick::f_meas, yardstick::pr_auc, yardstick::spec, yardstick::sens, - yardstick::bal_accuracy, yardstick::mcc - ), + # pr_auc, spec and sens need a Resistant event; only MCC chooses the model. + metrics = if (multi) { + yardstick::metric_set(yardstick::f_meas, yardstick::bal_accuracy, yardstick::mcc) + } else { + yardstick::metric_set( + yardstick::f_meas, yardstick::pr_auc, yardstick::spec, yardstick::sens, + yardstick::bal_accuracy, yardstick::mcc + ) + }, control = tune::control_grid(allow_par = FALSE) ) best <- tune::select_best(tuned, metric = "mcc") @@ -154,11 +165,18 @@ # Ported from old amRml's .calculate*() helpers, keeping AUPRC. Resistant is the event. # Each metric is rounded to 2, and nMCC and log2(AUPRC / prior) use the rounded values. -.taskMetrics <- function(scored) { +# With more than two labels, old amRml reported only MCC and nMCC; the rest are NA. +.taskMetrics <- function(scored, multi = FALSE) { estimate <- function(metric, ...) { round(as.numeric(metric(scored, truth = "phenotype", ...)$.estimate), 2) } mcc <- estimate(yardstick::mcc, estimate = ".pred_class") + if (multi) { + return(tibble::tibble( + mcc = mcc, nmcc = round((mcc + 1) / 2, 2), spec = NA_real_, sens = NA_real_, + log2_apop = NA_real_, f1 = NA_real_, bal_acc = NA_real_, auprc = NA_real_ + )) + } auprc <- estimate(yardstick::pr_auc, ".pred_Resistant") tibble::tibble( @@ -175,16 +193,28 @@ # Ported from old amRml's .viGlmnet(): the non-zero coefficients at the chosen penalty, on # normalized features. glmnet models the second level, so "POS" points to Susceptible. -.taskImportance <- function(fit, penalty, features) { - coefs <- stats::coef(parsnip::extract_fit_engine(fit), s = penalty)[, 1] - coefs <- coefs[names(coefs) != "(Intercept)" & coefs != 0] - coefs <- coefs[order(-abs(coefs))] +# A multinomial model has coefficients per label, as old amRml's multi-class branch kept; +# there "POS" points to that row's `label`. +.taskImportance <- function(fit, penalty, features, multi = FALSE) { + table <- function(coefs, label) { + coefs <- coefs[names(coefs) != "(Intercept)" & coefs != 0] + coefs <- coefs[order(-abs(coefs))] + tibble::tibble( + label = rep(label, length(coefs)), + Variable = features[as.integer(sub("^x", "", names(coefs)))], + Importance = abs(unname(coefs)), + Sign = ifelse(unname(coefs) > 0, "POS", "NEG") + ) + } - tibble::tibble( - Variable = features[as.integer(sub("^x", "", names(coefs)))], - Importance = abs(unname(coefs)), - Sign = ifelse(unname(coefs) > 0, "POS", "NEG") - ) + coefs <- stats::coef(parsnip::extract_fit_engine(fit), s = penalty) + if (!multi) { + return(table(coefs[, 1], NA_character_)) + } + dplyr::bind_rows(lapply(names(coefs), function(label) { + per_label <- coefs[[label]][, 1] + table(per_label, label) + })) } # One task: shuffle if asked, split, tune, fit, and score on the test genomes. @@ -192,6 +222,10 @@ started <- Sys.time() matrix <- .readMatrixWide(path, rows, columns, encoding, cache) rows <- matrix$rows + # Labels other than Resistant and Susceptible are MDR combinations. Eligibility keeps + # only labels with training genomes, so every test label is among them. + multi <- !all(rows$phenotype %in% .PHENOTYPES) + labels <- if (multi) sort(unique(rows$phenotype), method = "radix") else .PHENOTYPES if (task$labels == "shuffled") { rows <- .shuffleLabels(rows) } @@ -204,23 +238,27 @@ sparsity <- 1 - Matrix::nnzero(train_values) / (as.double(nrow(train_values)) * ncol(train_values)) # Is that fraction greater than or equal to our cutoff? If so, that's sparse use_sparse <- sparsity >= .SPARSE_THRESHOLD - train <- .modelData(train_values, rows$phenotype[split$train], sparse = use_sparse) + train <- .modelData(train_values, rows$phenotype[split$train], + sparse = use_sparse, labels = labels) test <- .modelData(matrix$values[split$test, , drop = FALSE], - rows$phenotype[split$test], sparse = use_sparse) - fitted <- .fitModel(train, task$n_fold, sparse = use_sparse) + rows$phenotype[split$test], sparse = use_sparse, labels = labels) + fitted <- .fitModel(train, task$n_fold, sparse = use_sparse, multi = multi) scored <- parsnip::augment(fitted$fit, test) - metrics <- .taskMetrics(scored) - importance <- .taskImportance(fitted$fit, fitted$penalty, matrix$features) + metrics <- .taskMetrics(scored, multi = multi) + importance <- .taskImportance(fitted$fit, fitted$penalty, matrix$features, multi = multi) input <- rows$role == "train" + # The predicted label is the most probable, so its probability is the row's largest. + probabilities <- as.matrix(scored[setdiff(grep("^[.]pred_", names(scored), value = TRUE), ".pred_class")]) performance <- tibble::tibble( task_id = task$task_id, num_obs = sum(input), - res_prop = round(mean(rows$phenotype[input] == "Resistant"), 2), + res_prop = if (multi) NA_real_ else round(mean(rows$phenotype[input] == "Resistant"), 2), n_test = length(split$test), - res_prop_test = round(mean(test$phenotype == "Resistant"), 2), + res_prop_test = if (multi) NA_real_ else round(mean(test$phenotype == "Resistant"), 2), n_feat = length(matrix$features), + n_labels = length(labels), model = "LR", n_feats_returned = nrow(importance), n_fold = task$n_fold, @@ -236,11 +274,13 @@ genome_id = rows$genome_id[split$test], truth = as.character(scored$phenotype), predicted = as.character(scored$.pred_class), - prob_resistant = scored$.pred_Resistant + prob_resistant = if (multi) NA_real_ else scored$.pred_Resistant, + prob_predicted = apply(probabilities, 1, max) ) # Metrics that can't be computed, such as MCC when every test genome gets one prediction. - undefined <- names(metrics)[vapply(metrics, is.na, logical(1))] + computed <- if (multi) metrics[c("mcc", "nmcc")] else metrics + undefined <- names(computed)[vapply(computed, is.na, logical(1))] if (length(undefined) && dplyr::n_distinct(scored$.pred_class) == 1) { undefined <- paste0( toString(undefined), "; every test genome was predicted ", scored$.pred_class[[1]] @@ -280,7 +320,12 @@ #' keeps old amRml's columns and rounding, adding `auprc`; #' `n_feat` counts the matrix's features (the recipe may drop more). #' `importance` has old amRml's `Variable`, `Importance` and `Sign` ("POS" -#' points to Susceptible). +#' points to Susceptible), plus `label`: for the MDR combination target, +#' which multinomial regression fits with coefficients per label, the label +#' a row's coefficient belongs to, and "POS" points to it (`NA` otherwise). +#' For that target, `performance` has only `mcc` and `nmcc`, as in old amRml, +#' and `res_prop` and `prob_resistant` are `NA`. `n_labels` counts each +#' task's labels and `prob_predicted` is the predicted label's probability. #' #' @section Errors: #' All errors have class `amrml_error`, plus `amrml_invalid_argument`: `tasks` @@ -367,19 +412,19 @@ fitModels <- function(tasks, overwrite = FALSE, fit$performance <- tibble::tibble( task_id = character(), num_obs = integer(), res_prop = numeric(), n_test = integer(), res_prop_test = numeric(), n_feat = integer(), - model = character(), n_feats_returned = integer(), n_fold = integer(), + n_labels = integer(), model = character(), n_feats_returned = integer(), n_fold = integer(), fit_penalty = numeric(), fit_mixture = numeric(), mcc = numeric(), nmcc = numeric(), spec = numeric(), sens = numeric(), log2_apop = numeric(), f1 = numeric(), bal_acc = numeric(), auprc = numeric(), run_time_sec = numeric(), seed = integer(), date = character() ) fit$importance <- tibble::tibble( - task_id = character(), Variable = character(), Importance = numeric(), - Sign = character() + task_id = character(), label = character(), Variable = character(), + Importance = numeric(), Sign = character() ) fit$predictions <- tibble::tibble( task_id = character(), genome_id = character(), truth = character(), - predicted = character(), prob_resistant = numeric() + predicted = character(), prob_resistant = numeric(), prob_predicted = numeric() ) } diff --git a/man/fitModels.Rd b/man/fitModels.Rd index 4e59377..aa60323 100644 --- a/man/fitModels.Rd +++ b/man/fitModels.Rd @@ -4,12 +4,20 @@ \alias{fitModels} \title{Fit the listed models} \usage{ -fitModels(tasks, overwrite = FALSE) +fitModels( + tasks, + overwrite = FALSE, + BPPARAM = BiocParallel::SnowParam(workers = 8L, tasks = 128L) +) } \arguments{ \item{tasks}{An \code{amr_tasks} from \code{\link[=expandTasks]{expandTasks()}}.} \item{overwrite}{Replace results already in the output folder.} + +\item{BPPARAM}{A \link[BiocParallel:BiocParallelParam]{BiocParallel::BiocParallelParam} backend; by default, +a two-worker socket cluster with 32 scheduling chunks. Use +\code{\link[BiocParallel:SerialParam]{BiocParallel::SerialParam()}} to run sequentially.} } \value{ An \code{amr_fit} list, also written to the output folder: \code{tasks} (with @@ -18,10 +26,15 @@ An \code{amr_fit} list, also written to the output folder: \code{tasks} (with keeps old amRml's columns and rounding, adding \code{auprc}; \code{n_feat} counts the matrix's features (the recipe may drop more). \code{importance} has old amRml's \code{Variable}, \code{Importance} and \code{Sign} ("POS" -points to Susceptible). +points to Susceptible), plus \code{label}: for the MDR combination target, +which multinomial regression fits with coefficients per label, the label +a row's coefficient belongs to, and "POS" points to it (\code{NA} otherwise). +For that target, \code{performance} has only \code{mcc} and \code{nmcc}, as in old amRml, +and \code{res_prop} and \code{prob_resistant} are \code{NA}. \code{n_labels} counts each +task's labels and \code{prob_predicted} is the predicted label's probability. } \description{ -Fits every task from \code{\link[=expandTasks]{expandTasks()}}, one at a time, with old amRml's method: +Fits tasks from \code{\link[=expandTasks]{expandTasks()}} in parallel with old amRml's method: glmnet logistic regression over 10 penalties and 6 mixtures, tuned by stratified \code{n_fold} cross-validation and chosen by MCC, then scored on the task's test genomes. @@ -35,8 +48,8 @@ split. \section{Errors}{ All errors have class \code{amrml_error}, plus \code{amrml_invalid_argument}: \code{tasks} -is not from \code{\link[=expandTasks]{expandTasks()}}, or results already exist and \code{overwrite} is -\code{FALSE}. +is not from \code{\link[=expandTasks]{expandTasks()}}, \code{BPPARAM} is not a BiocParallel backend, +or results already exist and \code{overwrite} is \code{FALSE}. } \examples{ diff --git a/tests/testthat/test-fit-models.R b/tests/testthat/test-fit-models.R new file mode 100644 index 0000000..4492666 --- /dev/null +++ b/tests/testthat/test-fit-models.R @@ -0,0 +1,81 @@ +# Tests for fitting the MDR combination target, whose labels aren't Resistant and Susceptible. + +# `n` genomes per label, tested on DRA (class CLX) and DRB (class CLY): resistant to CLX only, +# to both, or to neither. geneX marks CLX resistance, geneY CLY resistance, and geneN alternates +# with no signal. +fxLabelled <- function(n) { + genomes <- sprintf("g%03d", seq_len(3 * n)) + label <- rep(c("x", "xy", "s"), each = n) + metadata <- fxMetadataRows( + genome = rep(genomes, each = 2), + drug = rep(c("DRA", "DRB"), 3 * n), + class = rep(c("CLX", "CLY"), 3 * n), + phenotype = ifelse( + c(rbind(label %in% c("x", "xy"), label == "xy")), "Resistant", "Susceptible" + ) + ) + genes <- tibble::tibble( + genome_id = rep(genomes, 3), + gene = rep(c("geneX", "geneY", "geneN"), each = 3 * n), + value = c( + as.integer(label %in% c("x", "xy")), as.integer(label == "xy"), rep_len(0:1, 3 * n) + ) + ) + list(metadata = metadata, parquets = list(gene_count = genes)) +} + +# The combination target's matrices from `fxLabelled(n)`, built and fitted. +fxCombinationFit <- function(n = 20, env = parent.frame()) { + data <- fxLabelled(n) + eligibility <- fxEligibility( + data$metadata, modes = "vanilla", parquets = data$parquets, + min_genomes = n, n_fold = 2, min_label_genomes = n, env = env + ) + eligibility$scopes <- eligibility$scopes[eligibility$scopes$target == "resistant_classes", ] + eligibility$members <- eligibility$members[eligibility$members$target == "resistant_classes", ] + plan <- planMatrices(eligibility) + built <- buildMatrices(plan, withr::local_tempdir(.local_envir = env)) + fitModels(expandTasks(built), BPPARAM = BiocParallel::SerialParam()) +} + +test_that("folds are assigned for every label, not only Resistant and Susceptible", { + phenotype <- factor(rep(c("CLX", "CLX_CLY", "Susceptible"), each = 4)) + fold <- withr::with_seed(1, .stratifiedFolds(phenotype, 2)) + + expect_true(all(fold %in% 1:2)) + expect_equal(as.vector(table(phenotype, fold)), rep(2, 6)) +}) + +test_that("Resistant and Susceptible fold assignment is unchanged", { + phenotype <- factor(rep(c("Resistant", "Susceptible"), c(5, 7)), levels = .PHENOTYPES) + fold <- withr::with_seed(1, .stratifiedFolds(phenotype, 3)) + + expected <- withr::with_seed(1, { + expected <- integer(12) + expected[1:5] <- sample(rep_len(1:3, 5)) + expected[6:12] <- sample(rep_len(1:3, 7)) + expected + }) + expect_equal(fold, expected) +}) + +test_that("the combination target is fitted by multinomial regression, one task per matrix", { + fit <- fxCombinationFit() + + expect_equal(fit$tasks$status, rep("fitted", nrow(fit$tasks))) + expect_true(all(is.na(fit$tasks$rule_id))) + expect_equal(unique(fit$performance$n_labels), 3) + expect_false(anyNA(fit$performance$mcc)) + expect_true(all(is.na(fit$performance$res_prop))) + expect_true(all(is.na(fit$performance$auprc))) +}) + +test_that("importance is per label, and predictions carry the predicted label's probability", { + fit <- fxCombinationFit() + + expect_true(all(fit$importance$label %in% c("CLX", "CLX_CLY", "Susceptible"))) + expect_true("geneX" %in% fit$importance$Variable) + expect_true(all(is.na(fit$predictions$prob_resistant))) + expect_true(all(fit$predictions$prob_predicted >= 1 / 3 & fit$predictions$prob_predicted <= 1)) + expect_true(all(fit$predictions$truth %in% c("CLX", "CLX_CLY", "Susceptible"))) +}) From 803ba0207e9dec04b302d6be9c415f36d086c17b Mon Sep 17 00:00:00 2001 From: Emily Boyer Date: Fri, 16 Oct 2026 11:31:42 -0600 Subject: [PATCH 05/14] docs: record MDR combination results against old amRml --- design/mdr.md | 38 ++++++++++++++++++++++++++------------ 1 file changed, 26 insertions(+), 12 deletions(-) diff --git a/design/mdr.md b/design/mdr.md index 13ce9b2..0af31ea 100644 --- a/design/mdr.md +++ b/design/mdr.md @@ -1,6 +1,6 @@ # MDR: modeling combinations of resistant classes -Status: draft spec. The MDR yes/no target is implemented; the combination target is not yet. +Status: both MDR targets are implemented. See "Checking against old amRml" for results. ## Goal @@ -162,17 +162,31 @@ should be included in the global rename. ### 5. Checking against old amRml -1. S. flexneri, `modes = "vanilla"`: same 288 genomes and 4 labels as the old - MDR matrices, and the same features where the old code's all-genome variance - filter agrees. Compare fits on the old and new matrices with seed 123. -2. S. argenteus: old builds a 1-label MDR model; new records `single_phenotype`. -3. E. faecium: run every mode and check the dropped-label counts by hand for - one cross scope. +Done on S. flexneri (amRdata main, October 2026), with old amRml from #42: + +1. **Matrices:** all 13 old MDR matrices (`.parquet2MDRMatrix()`, `min_n = 25`) + are identical to the new whole-dataset combination matrices: the same 288 + genomes and labels, the same features, and the same non-zero values. In the + whole-dataset mode every scope genome is a training genome at matrix time + (the holdout happens at fitting), so selecting features from training genomes + gives the old all-genome filter. +2. **Fit, AMRFinder binary, seed 123:** old `runMLPipeline()` and new + `fitModels()` choose the same penalty (1e-4) and mixture (0.2), test on the + same 59 genomes, predict the same label for every one, and score MCC 0.85. + All 521 non-zero coefficients match exactly. + +Still to check: + +3. **S. argenteus:** old builds a 1-label MDR model (`PEN`, 65 genomes); new + records `single_phenotype` and drops the other 8 genomes. Its yes/no target + (5 MDR, 68 not) fails the fold check. + +Still to check: E. faecium in every mode, once its ORB is ready. ## Order of work -1. Profile: `mdr` yes/no rows, tests. **Done.** Nothing after the profile needs to change for it, so it runs in every mode. -2. Profile: combination rows, tests. -3. Eligibility: label filter, rules, tests. -4. Fitting: multinomial path, metrics, importance, tests. -5. Comparison with old amRml on S. flexneri. +1. Profile: `mdr` yes/no rows, tests. **Done.** +2. Profile: combination rows, tests. **Done.** +3. Eligibility: label filter, rules, tests. **Done.** +4. Fitting: multinomial path, metrics, importance, tests. **Done.** +5. Comparison with old amRml on S. flexneri. **Done.** From 3dd51bb70e9a004d985d61e57b56cbc6f0c6094c Mon Sep 17 00:00:00 2001 From: Emily Boyer Date: Mon, 19 Oct 2026 09:15:05 -0400 Subject: [PATCH 06/14] chore: exclude design/ from the package build --- .Rbuildignore | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.Rbuildignore b/.Rbuildignore index 1ec22b7..91cd97b 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -7,4 +7,4 @@ ^CODE_OF_CONDUCT\.md$ ^doc$ ^Meta$ -^results$ \ No newline at end of file +^results$^design$ From 191cf56efc8ea7b713b368da04c472efa2e34c07 Mon Sep 17 00:00:00 2001 From: Emily Boyer Date: Mon, 19 Oct 2026 09:16:57 -0400 Subject: [PATCH 07/14] fix: put design/ on its own .Rbuildignore line --- .Rbuildignore | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.Rbuildignore b/.Rbuildignore index 91cd97b..ebfabc7 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -7,4 +7,5 @@ ^CODE_OF_CONDUCT\.md$ ^doc$ ^Meta$ -^results$^design$ +^results$ +^design$ From 5bed52480781a535115a2329177aae8a72d925b2 Mon Sep 17 00:00:00 2001 From: Emily Boyer Date: Mon, 19 Oct 2026 09:55:33 -0400 Subject: [PATCH 08/14] refactor: declare multi-class targets in the registry Adds .MULTICLASS_TARGETS and .isMulticlass() next to .MODES. Eligibility and fitting read it instead of matching the MDR combination target by name or inferring it from labels; fitting now errors if a Resistant/Susceptible target has other phenotypes. Results on S. flexneri are unchanged. --- R/eligibility.R | 22 +++++++++++----------- R/fit-models.R | 22 ++++++++++++++-------- R/orb-profile.R | 4 ---- R/registry.R | 18 ++++++++++++++++++ design/mdr.md | 9 ++++++++- tests/testthat/test-fit-models.R | 7 +++++++ 6 files changed, 58 insertions(+), 24 deletions(-) diff --git a/R/eligibility.R b/R/eligibility.R index 6073e4a..fc8e114 100644 --- a/R/eligibility.R +++ b/R/eligibility.R @@ -122,22 +122,22 @@ .asMembers(train, test) } -# Combination scopes keep the labels with at least `min_label_genomes` training genomes, and +# Multi-class scopes keep the labels with at least `min_label_genomes` training genomes, and # drop every genome, training or test, with another label: a test genome whose label the # model never trained on can't be predicted. Other scopes keep every genome. .keepLabels <- function(members, min_label_genomes) { members$kept <- TRUE - combination <- members$unit == "mdr" & members$target == .MDR_COMBINATION - if (!any(combination)) { + multiclass <- .isMulticlass(members$unit, members$target) + if (!any(multiclass)) { return(members) } keys <- c(.SCOPE_KEYS, "phenotype") - train <- members[combination & members$role == "train", keys] + train <- members[multiclass & members$role == "train", keys] n <- dplyr::count(train, dplyr::across(dplyr::all_of(keys))) labels <- n[n$n >= min_label_genomes, keys] key <- function(x) do.call(paste, c(unname(as.list(x[keys])), sep = "\r")) - members$kept[combination] <- key(members[combination, ]) %in% key(labels) + members$kept[multiclass] <- key(members[multiclass, ]) %in% key(labels) members } @@ -247,13 +247,13 @@ # A class with data for only one of its drugs is that drug's model again. duplicate <- scopes$unit == "drug_class" & scopes$n_drugs %in% 1 - # A combination scope has a label per combination, so its "rarer phenotype" is its - # rarest label, and it needs two labels rather than both phenotypes. - combination <- scopes$unit == "mdr" & scopes$target == .MDR_COMBINATION + # A multi-class scope's "rarer phenotype" is its rarest label, and it needs two labels + # rather than both phenotypes. + multiclass <- .isMulticlass(scopes$unit, scopes$target) rarer <- ifelse( - combination, scopes$n_rarest, pmin(scopes$n_resistant, scopes$n_susceptible) + multiclass, scopes$n_rarest, pmin(scopes$n_resistant, scopes$n_susceptible) ) - single <- ifelse(combination, scopes$n_labels < 2, rarer == 0) + single <- ifelse(multiclass, scopes$n_labels < 2, rarer == 0) enough_for_cv <- if (has_test) { rarer >= settings$n_fold } else { @@ -261,7 +261,7 @@ } test_ok <- !has_test | ( scopes$test_n_genomes >= settings$min_test_genomes & - (combination | + (multiclass | pmin(scopes$test_n_resistant, scopes$test_n_susceptible) >= settings$min_test_minority) ) diff --git a/R/fit-models.R b/R/fit-models.R index 081e95b..ac93195 100644 --- a/R/fit-models.R +++ b/R/fit-models.R @@ -1,6 +1,6 @@ # Fitting a model for each task, with old amRml's method: glmnet logistic regression, -# tuned by stratified cross-validation and chosen by MCC. A task whose labels aren't -# Resistant and Susceptible (the MDR combination target) gets multinomial regression. +# tuned by stratified cross-validation and chosen by MCC. A multi-class target (the MDR +# combination target; see `.MULTICLASS_TARGETS`) gets multinomial regression. # The phenotypes, with Resistant first so it is the event yardstick scores. .PHENOTYPES <- c("Resistant", "Susceptible") @@ -269,15 +269,17 @@ split(seq_len(nrow(listed)), factor(leader, levels = unique(leader))) } -# One task: shuffle if asked, split, tune, fit, and score on the test genomes. +# One task: shuffle if asked, split, tune, fit, and score on the test genomes. `multi` says +# whether the task's target is multi-class. .fitTask <- function(task, path, rows, columns, encoding, cache, - fitted = NULL, return_fitted = FALSE) { + fitted = NULL, return_fitted = FALSE, multi = FALSE) { started <- Sys.time() matrix <- .readMatrixWide(path, rows, columns, encoding, cache) rows <- matrix$rows - # Labels other than Resistant and Susceptible are MDR combinations. Eligibility keeps - # only labels with training genomes, so every test label is among them. - multi <- !all(rows$phenotype %in% .PHENOTYPES) + # Eligibility keeps only labels with training genomes, so every test label is among them. + if (!multi && !all(rows$phenotype %in% .PHENOTYPES)) { + stop("A Resistant/Susceptible target has other phenotypes.") + } labels <- if (multi) sort(unique(rows$phenotype), method = "radix") else .PHENOTYPES if (task$labels == "shuffled") { rows <- .shuffleLabels(rows) @@ -426,6 +428,9 @@ fitModels <- function(tasks, overwrite = FALSE, index <- readRDS(file.path(out_dir, "matrix_index.rds")) cache <- new.env(parent = emptyenv()) matrix_positions <- match(listed$matrix_id, tasks$built$matrices$matrix_id) + multiclass <- .isMulticlass( + tasks$built$matrices$unit[matrix_positions], tasks$built$matrices$target[matrix_positions] + ) groups <- if (reuse_fits) { .fitTaskGroups(listed, tasks$built$matrices, index) @@ -453,7 +458,8 @@ fitModels <- function(tasks, overwrite = FALSE, encoding = tasks$built$matrices$encoding[[matrix_positions[[i]]]], cache = cache, fitted = shared_fit, - return_fitted = keep_fit + return_fitted = keep_fit, + multi = multiclass[[i]] ), .rng_kind = "Mersenne-Twister", .rng_normal_kind = "Inversion", diff --git a/R/orb-profile.R b/R/orb-profile.R index 268ad37..868a507 100644 --- a/R/orb-profile.R +++ b/R/orb-profile.R @@ -1,9 +1,5 @@ # Profiling an amRdata ORB: each genome's phenotypes, and the counts per target. -# The MDR target labelling each genome with its combination of resistant classes, as -# amRdata's `resistant_classes` does. -.MDR_COMBINATION <- "resistant_classes" - # amRdata's metadata.parquet columns that profiling needs, named as they are used here. .METADATA_COLUMNS <- c( genome_id = "genome.genome_id", diff --git a/R/registry.R b/R/registry.R index 6b9140e..d16cb33 100644 --- a/R/registry.R +++ b/R/registry.R @@ -3,6 +3,24 @@ # Strata, as amRml names them, and the genome-level amRdata column holding each. .STRATA <- c(year = "year_bin", country = "country_abbr") +# The MDR target labelling each genome with its combination of resistant classes, as +# amRdata's `resistant_classes` does. +.MDR_COMBINATION <- "resistant_classes" + +# Targets whose phenotype is one label among several, not Resistant or Susceptible. +# Eligibility counts their labels, and fitting uses multinomial regression. Every other +# target is Resistant/Susceptible. +.MULTICLASS_TARGETS <- tibble::tribble( + ~unit, ~target, ~description, + "mdr", .MDR_COMBINATION, "Each genome's combination of resistant classes." +) + +# Whether each `unit` and `target` pair is a multi-class target. +.isMulticlass <- function(unit, target) { + paste(unit, target, sep = "\r") %in% + paste(.MULTICLASS_TARGETS$unit, .MULTICLASS_TARGETS$target, sep = "\r") +} + # The modeling modes of the amRml v1 design. `grouping` is what a mode splits genomes # by: a stratum, drugs, or nothing. .MODES <- tibble::tribble( diff --git a/design/mdr.md b/design/mdr.md index 0af31ea..46061d8 100644 --- a/design/mdr.md +++ b/design/mdr.md @@ -63,6 +63,11 @@ which matters for the cross modes below. ## Design, by stage +The combination target is declared once, in `.MULTICLASS_TARGETS` in +`R/registry.R`, next to `.MODES`. Eligibility and fitting read it through +`.isMulticlass(unit, target)` instead of checking names or guessing from +labels. A new multi-class target only needs a row there and a profile builder. + ### 1. Profile (`profileORB()`) Add two kinds of `unit = "mdr"` rows to `profile$phenotypes`: @@ -129,7 +134,9 @@ training can't be predicted correctly, so they are dropped from the test set The combination label has the same blind spot as the yes/no one: a genome tested in only one class is labelled from that class alone. Old amRml did the -same, so the port keeps it, but the per-scope record should count these genomes. +same, so the port keeps these genomes on purpose, while the yes/no target leaves +them out as `undetermined`. **Open:** should the combination target also leave +out `undetermined` genomes? That would change the results compared with old amRml. ### 3. Matrices (`planMatrices()`, `buildMatrices()`) diff --git a/tests/testthat/test-fit-models.R b/tests/testthat/test-fit-models.R index 4492666..d3f3d45 100644 --- a/tests/testthat/test-fit-models.R +++ b/tests/testthat/test-fit-models.R @@ -38,6 +38,13 @@ fxCombinationFit <- function(n = 20, env = parent.frame()) { fitModels(expandTasks(built), BPPARAM = BiocParallel::SerialParam()) } +test_that("only the declared MDR combination target is multi-class", { + expect_equal( + .isMulticlass(c("mdr", "mdr", "drug", "drug_class"), c("resistant_classes", "mdr", "DRA", "CLX")), + c(TRUE, FALSE, FALSE, FALSE) + ) +}) + test_that("folds are assigned for every label, not only Resistant and Susceptible", { phenotype <- factor(rep(c("CLX", "CLX_CLY", "Susceptible"), each = 4)) fold <- withr::with_seed(1, .stratifiedFolds(phenotype, 2)) From 250d9c4ac1f011f48383fcdcd6a5cb3a61407836 Mon Sep 17 00:00:00 2001 From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com> Date: Wed, 21 Oct 2026 12:03:47 -0600 Subject: [PATCH 09/14] print MDR in eligibility --- R/eligibility.R | 11 +++++++++++ tests/testthat/test-eligibility.R | 6 ++++++ 2 files changed, 17 insertions(+) diff --git a/R/eligibility.R b/R/eligibility.R index fc8e114..c51dfff 100644 --- a/R/eligibility.R +++ b/R/eligibility.R @@ -473,6 +473,17 @@ print.amr_eligibility <- function(x, ...) { s$min_test_genomes, "genomes,", s$min_test_minority, "of the rarer phenotype\n" ) + combination <- x$scopes[.isMulticlass(x$scopes$unit, x$scopes$target), , drop = FALSE] + cat( + " MDR : ", x$profile$settings$mdr_min_classes, "+ classes; combination labels need ", + s$min_label_genomes, " training genomes", + if (nrow(combination)) { + paste0("; ", sum(combination$eligible), " of ", nrow(combination), " combination scopes eligible") + }, + "\n", + sep = "" + ) + labels <- format(x$modes$mode_id) for (i in seq_len(nrow(x$modes))) { mode <- x$modes[i, , drop = FALSE] diff --git a/tests/testthat/test-eligibility.R b/tests/testthat/test-eligibility.R index 7420c82..1ed84ff 100644 --- a/tests/testthat/test-eligibility.R +++ b/tests/testthat/test-eligibility.R @@ -502,3 +502,9 @@ test_that("the print shows each mode's result", { expect_output(print(eligibility), "stratified_year +: not run \\(too_few_groups\\)") expect_output(print(eligibility), "lodo +: not run \\(mode_not_supported\\)") }) + +test_that("the print shows the MDR settings", { + eligibility <- fxEligibility(modes = "vanilla") + + expect_output(print(eligibility), "MDR +: 2\\+ classes; combination labels need 25 training genomes") +}) From 9249ebf5d79cd9e9fa8df1dea1b72e14a1c85621 Mon Sep 17 00:00:00 2001 From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com> Date: Wed, 21 Oct 2026 12:55:15 -0600 Subject: [PATCH 10/14] print model units --- R/fit-models.R | 9 +++++++++ tests/testthat/test-fit-models.R | 6 ++++++ 2 files changed, 15 insertions(+) diff --git a/R/fit-models.R b/R/fit-models.R index ac93195..7b545a7 100644 --- a/R/fit-models.R +++ b/R/fit-models.R @@ -578,6 +578,15 @@ print.amr_fit <- function(x, ...) { cat("", orb$dataset_id, "-", orb$dataset_label, "\n") cat(" in", x$built$out_dir, "\n") cat(" ", nrow(t), " tasks: ", counts(t), "\n", sep = "") + units <- table(factor( + x$built$matrices$unit[match(t$matrix_id, x$built$matrices$matrix_id)], + c("drug", "drug_class", "mdr") + )) + cat( + " by target: ", units[["drug"]], " drug, ", units[["drug_class"]], " drug class, ", + units[["mdr"]], " MDR\n", + sep = "" + ) cat(.modeLines(names(by_mode), vapply(by_mode, counts, character(1))), sep = "\n") invisible(x) } diff --git a/tests/testthat/test-fit-models.R b/tests/testthat/test-fit-models.R index d3f3d45..787ddfb 100644 --- a/tests/testthat/test-fit-models.R +++ b/tests/testthat/test-fit-models.R @@ -86,3 +86,9 @@ test_that("importance is per label, and predictions carry the predicted label's expect_true(all(fit$predictions$prob_predicted >= 1 / 3 & fit$predictions$prob_predicted <= 1)) expect_true(all(fit$predictions$truth %in% c("CLX", "CLX_CLY", "Susceptible"))) }) + +test_that("the print counts tasks by drug, drug class and MDR", { + fit <- fxCombinationFit() + + expect_output(print(fit), "by target: 0 drug, 0 drug class, 1 MDR") +}) From 6ca9bfe808092fcdccddbd997612abe5b0021c0d Mon Sep 17 00:00:00 2001 From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com> Date: Wed, 21 Oct 2026 13:11:19 -0600 Subject: [PATCH 11/14] add mdr parameters to workflow --- R/workflow.R | 8 ++++++-- man/planMatrices.Rd | 19 +++++++++++++++---- man/runModelingWorkflow.Rd | 11 ++++++++++- 3 files changed, 31 insertions(+), 7 deletions(-) diff --git a/R/workflow.R b/R/workflow.R index 3b9b235..ef84fb2 100644 --- a/R/workflow.R +++ b/R/workflow.R @@ -21,8 +21,10 @@ #' @param seeds Passed to [expandTasks()]. #' @param labels `"real"`, `"shuffled"` or both (the default), passed to #' [expandTasks()]. -#' @param min_genomes,n_fold,min_test_genomes,min_test_minority Passed to +#' @param min_genomes,n_fold,min_test_genomes,min_test_minority,min_label_genomes Passed to #' [eligibleScopes()]. +#' @param mdr_min_classes Passed to [profileORB()]: fewest drug classes a genome +#' must be Resistant in to be MDR; at least 2. #' @param overwrite Replace an earlier run's matrices and results in `out_dir`; #' if `FALSE`, an `out_dir` holding earlier output is refused. The default run #' folder is always new. @@ -52,6 +54,7 @@ runModelingWorkflow <- function(orb_dir, min_test_genomes = 5, min_test_minority = 2, min_label_genomes = 25, + mdr_min_classes = 2, overwrite = FALSE, BPPARAM = BiocParallel::SnowParam(workers = 8L, tasks = 128L), reuse_fits = TRUE, @@ -62,12 +65,13 @@ runModelingWorkflow <- function(orb_dir, .checkArgSeeds(seeds, call = call) .checkArgLabels(labels, call = call) .checkArgClass(BPPARAM, "BPPARAM", "BiocParallelParam", "BiocParallel::SnowParam()", call = call) + .checkArgCount(mdr_min_classes, "mdr_min_classes", min = 2, call = call) .checkArgFlag(reuse_fits, "reuse_fits", call = call) orb <- readORB(orb_dir) out_dir <- out_dir %||% .defaultOutDir(orb$directory) eligibility <- eligibleScopes( - profileORB(orb), + profileORB(orb, mdr_min_classes = mdr_min_classes), modes = modes, min_genomes = min_genomes, n_fold = n_fold, min_test_genomes = min_test_genomes, min_test_minority = min_test_minority, min_label_genomes = min_label_genomes diff --git a/man/planMatrices.Rd b/man/planMatrices.Rd index 1ff1720..d031a69 100644 --- a/man/planMatrices.Rd +++ b/man/planMatrices.Rd @@ -4,13 +4,19 @@ \alias{planMatrices} \title{Plan the matrices to build} \usage{ -planMatrices(eligibility, modes = NULL) +planMatrices(eligibility, modes = NULL, feature_types = "standard") } \arguments{ \item{eligibility}{An \code{amr_eligibility} from \code{\link[=eligibleScopes]{eligibleScopes()}}.} \item{modes}{Mode IDs from \code{\link[=modelingModes]{modelingModes()}} to plan matrices for, or \code{NULL} for every mode in \code{eligibility}.} + +\item{feature_types}{A preset (\code{"quick"}, \code{"standard"}, or \code{"full"}) or an +explicit character vector of feature types. Default \code{"standard"}. +Quick uses gene, protein, and Pfam; Standard also includes COG and +AMRFinder; Full additionally includes struct and DefenseCas. Presets use +only tables present in the ORB. Explicit selections require exact names.} } \value{ An \code{amr_matrix_plan} list: @@ -23,13 +29,15 @@ the scope's \code{n_train} and \code{n_test} genomes (\code{NA} without a test s } \description{ Lists every matrix the eligible scopes need, before any is built: one per -eligible scope, feature type and encoding. Each feature type is encoded as -binary, and also as counts when its table has values above 1. +eligible scope, selected feature type and encoding. Each feature type is +encoded as binary, and also as counts when its table has values above 1. +Feature types are selected before their tables are scanned or matrices built. } \section{Errors}{ All errors have class \code{amrml_error}, plus \code{amrml_invalid_argument} if -\code{eligibility} is not from \code{\link[=eligibleScopes]{eligibleScopes()}} or \code{modes} names an unknown mode. +\code{eligibility} is not from \code{\link[=eligibleScopes]{eligibleScopes()}}, \code{modes} names an unknown mode, +or \code{feature_types} is invalid. } \examples{ @@ -39,5 +47,8 @@ plan$matrices # Only some modes planMatrices(eligibility, modes = c("vanilla", "loto")) + +# Quick feature selection +planMatrices(eligibility, feature_types = "quick") } } diff --git a/man/runModelingWorkflow.Rd b/man/runModelingWorkflow.Rd index 99f7284..50fc246 100644 --- a/man/runModelingWorkflow.Rd +++ b/man/runModelingWorkflow.Rd @@ -15,9 +15,11 @@ runModelingWorkflow( min_test_genomes = 5, min_test_minority = 2, min_label_genomes = 25, + mdr_min_classes = 2, overwrite = FALSE, BPPARAM = BiocParallel::SnowParam(workers = 8L, tasks = 128L), - reuse_fits = TRUE + reuse_fits = TRUE, + feature_types = "standard" ) } \arguments{ @@ -36,6 +38,9 @@ runModelingWorkflow( \item{min_genomes, n_fold, min_test_genomes, min_test_minority, min_label_genomes}{Passed to \code{\link[=eligibleScopes]{eligibleScopes()}}.} +\item{mdr_min_classes}{Passed to \code{\link[=profileORB]{profileORB()}}: fewest drug classes a genome +must be Resistant in to be MDR; at least 2.} + \item{overwrite}{Replace an earlier run's matrices and results in \code{out_dir}; if \code{FALSE}, an \code{out_dir} holding earlier output is refused. The default run folder is always new.} @@ -44,6 +49,10 @@ folder is always new.} eight-worker cluster with 128 chunks of work to chew through.} \item{reuse_fits}{Passed to \code{\link[=fitModels]{fitModels()}}; reuse compatible cross-drug fits.} + +\item{feature_types}{Feature types to model: \code{"quick"}, \code{"standard"} (default), +\code{"full"}, or an explicit character vector of ORB feature types. +Selected before matrices are planned or built; see \code{\link[=planMatrices]{planMatrices()}}.} } \value{ The \code{amr_fit} from \code{\link[=fitModels]{fitModels()}}; every earlier stage's result can be From 5c0c8f0bc3dda840b3aed6a760fa1d0d0a26f936 Mon Sep 17 00:00:00 2001 From: Emily Boyer Date: Mon, 26 Oct 2026 13:16:29 -0600 Subject: [PATCH 12/14] refactor: rename min_label_genomes to min_mdr_label_genomes Clearer for users, since MDR combinations are the only labels it applies to (review suggestion from @AbhirupaGhosh). --- R/eligibility.R | 20 ++++++++++---------- R/workflow.R | 8 ++++---- design/mdr.md | 4 ++-- man/eligibleScopes.Rd | 6 +++--- man/runModelingWorkflow.Rd | 4 ++-- tests/testthat/test-eligibility.R | 8 ++++---- tests/testthat/test-fit-models.R | 2 +- 7 files changed, 26 insertions(+), 26 deletions(-) diff --git a/R/eligibility.R b/R/eligibility.R index c51dfff..2dd967a 100644 --- a/R/eligibility.R +++ b/R/eligibility.R @@ -122,10 +122,10 @@ .asMembers(train, test) } -# Multi-class scopes keep the labels with at least `min_label_genomes` training genomes, and +# Multi-class scopes keep the labels with at least `min_mdr_label_genomes` training genomes, and # drop every genome, training or test, with another label: a test genome whose label the # model never trained on can't be predicted. Other scopes keep every genome. -.keepLabels <- function(members, min_label_genomes) { +.keepLabels <- function(members, min_mdr_label_genomes) { members$kept <- TRUE multiclass <- .isMulticlass(members$unit, members$target) if (!any(multiclass)) { @@ -135,7 +135,7 @@ keys <- c(.SCOPE_KEYS, "phenotype") train <- members[multiclass & members$role == "train", keys] n <- dplyr::count(train, dplyr::across(dplyr::all_of(keys))) - labels <- n[n$n >= min_label_genomes, keys] + labels <- n[n$n >= min_mdr_label_genomes, keys] key <- function(x) do.call(paste, c(unname(as.list(x[keys])), sep = "\r")) members$kept[multiclass] <- key(members[multiclass, ]) %in% key(labels) members @@ -331,7 +331,7 @@ #' @details #' Counts are of genomes, each with one phenotype. In the MDR combination target #' (`resistant_classes`), each combination is a label: a scope keeps the labels -#' with at least `min_label_genomes` training genomes and drops the genomes, +#' with at least `min_mdr_label_genomes` training genomes and drops the genomes, #' training or test, with any other (`n_dropped`, `test_n_dropped`). It then #' needs 2 labels, and its rarest label must fill the folds; its test set needs #' only `min_test_genomes`. A drug or class needs 2 of a @@ -352,7 +352,7 @@ #' @param n_fold Number of cross-validation folds. #' @param min_test_genomes Fewest genomes a separate test set needs. #' @param min_test_minority Fewest of the rarer phenotype a test set needs. -#' @param min_label_genomes Fewest training genomes an MDR combination needs to +#' @param min_mdr_label_genomes Fewest training genomes an MDR combination needs to #' stay a label in a scope (old amRml's `min_n`). #' #' @return An `amr_eligibility` list: @@ -376,7 +376,7 @@ eligibleScopes <- function(profile, n_fold = 5, min_test_genomes = 5, min_test_minority = 2, - min_label_genomes = 25) { + min_mdr_label_genomes = 25) { call <- rlang::current_env() # Check the arguments. @@ -385,7 +385,7 @@ eligibleScopes <- function(profile, .checkArgCount(n_fold, "n_fold", min = 2, call = call) .checkArgCount(min_test_genomes, "min_test_genomes", call = call) .checkArgCount(min_test_minority, "min_test_minority", call = call) - .checkArgCount(min_label_genomes, "min_label_genomes", call = call) + .checkArgCount(min_mdr_label_genomes, "min_mdr_label_genomes", call = call) modes <- modes %||% .MODES$mode_id .checkArgModes(modes, call = call) @@ -394,7 +394,7 @@ eligibleScopes <- function(profile, n_fold = n_fold, min_test_genomes = min_test_genomes, min_test_minority = min_test_minority, - min_label_genomes = min_label_genomes, + min_mdr_label_genomes = min_mdr_label_genomes, holdout = .CV_HOLDOUT ) @@ -408,7 +408,7 @@ eligibleScopes <- function(profile, judged <- lapply(seq_len(nrow(runnable)), function(i) { mode <- runnable[i, , drop = FALSE] built <- .buildScopes(mode, profile) - built$members <- .keepLabels(built$members, settings$min_label_genomes) + built$members <- .keepLabels(built$members, settings$min_mdr_label_genomes) decisions <- .judgeScopes(built, mode, profile, settings) eligible <- decisions[decisions$eligible, .SCOPE_KEYS] list( @@ -476,7 +476,7 @@ print.amr_eligibility <- function(x, ...) { combination <- x$scopes[.isMulticlass(x$scopes$unit, x$scopes$target), , drop = FALSE] cat( " MDR : ", x$profile$settings$mdr_min_classes, "+ classes; combination labels need ", - s$min_label_genomes, " training genomes", + s$min_mdr_label_genomes, " training genomes", if (nrow(combination)) { paste0("; ", sum(combination$eligible), " of ", nrow(combination), " combination scopes eligible") }, diff --git a/R/workflow.R b/R/workflow.R index 4149897..90f7ccd 100644 --- a/R/workflow.R +++ b/R/workflow.R @@ -22,7 +22,7 @@ #' @param seeds Passed to [expandTasks()]. #' @param labels `"real"`, `"shuffled"` or both (the default), passed to #' [expandTasks()]. -#' @param min_genomes,n_fold,min_test_genomes,min_test_minority,min_label_genomes Passed to +#' @param min_genomes,n_fold,min_test_genomes,min_test_minority,min_mdr_label_genomes Passed to #' [eligibleScopes()]. #' @param mdr_min_classes Passed to [profileORB()]: fewest drug classes a genome #' must be Resistant in to be MDR; at least 2. @@ -59,7 +59,7 @@ runModelingWorkflow <- function(orb_dir, n_fold = 5, min_test_genomes = 5, min_test_minority = 2, - min_label_genomes = 25, + min_mdr_label_genomes = 25, mdr_min_classes = 2, overwrite = FALSE, BPPARAM = BiocParallel::SnowParam(workers = 8L, tasks = 128L), @@ -84,7 +84,7 @@ runModelingWorkflow <- function(orb_dir, modes = modes, seeds = as.integer(seeds), labels = labels, min_genomes = min_genomes, n_fold = n_fold, min_test_genomes = min_test_genomes, min_test_minority = min_test_minority, - min_label_genomes = min_label_genomes, mdr_min_classes = mdr_min_classes, + min_mdr_label_genomes = min_mdr_label_genomes, mdr_min_classes = mdr_min_classes, reuse_fits = reuse_fits, feature_types = feature_types ) if (resume) { @@ -106,7 +106,7 @@ runModelingWorkflow <- function(orb_dir, profileORB(orb, mdr_min_classes = mdr_min_classes), modes = modes, min_genomes = min_genomes, n_fold = n_fold, min_test_genomes = min_test_genomes, min_test_minority = min_test_minority, - min_label_genomes = min_label_genomes + min_mdr_label_genomes = min_mdr_label_genomes ) built <- buildMatrices(planMatrices(eligibility, feature_types = feature_types), out_dir, overwrite = overwrite) diff --git a/design/mdr.md b/design/mdr.md index 46061d8..6d4d796 100644 --- a/design/mdr.md +++ b/design/mdr.md @@ -107,9 +107,9 @@ MDR uses the existing scope builders, so it gets scopes in every mode with a grouping of `none`, `year` or `country`. The drug-grouped modes (`cross_drug`, `lodo`) only take `unit == "drug"` rows, so they leave MDR out with no change. -New argument `min_label_genomes = 25` (old `min_n`). Per MDR scope: +New argument `min_mdr_label_genomes = 25` (old `min_n`). Per MDR scope: -1. Keep the labels with at least `min_label_genomes` **training** genomes, and +1. Keep the labels with at least `min_mdr_label_genomes` **training** genomes, and drop the scope's genomes (train and test) with any other label. On the whole dataset this gives exactly old amRml's genome set. Doing it per scope keeps a label with 3 genomes in one year from breaking that year's folds. diff --git a/man/eligibleScopes.Rd b/man/eligibleScopes.Rd index 8989873..b541ca0 100644 --- a/man/eligibleScopes.Rd +++ b/man/eligibleScopes.Rd @@ -11,7 +11,7 @@ eligibleScopes( n_fold = 5, min_test_genomes = 5, min_test_minority = 2, - min_label_genomes = 25 + min_mdr_label_genomes = 25 ) } \arguments{ @@ -27,7 +27,7 @@ eligibleScopes( \item{min_test_minority}{Fewest of the rarer phenotype a test set needs.} -\item{min_label_genomes}{Fewest training genomes an MDR combination needs to +\item{min_mdr_label_genomes}{Fewest training genomes an MDR combination needs to stay a label in a scope (old amRml's \code{min_n}).} } \value{ @@ -51,7 +51,7 @@ built or fitted. \details{ Counts are of genomes, each with one phenotype. In the MDR combination target (\code{resistant_classes}), each combination is a label: a scope keeps the labels -with at least \code{min_label_genomes} training genomes and drops the genomes, +with at least \code{min_mdr_label_genomes} training genomes and drops the genomes, training or test, with any other (\code{n_dropped}, \code{test_n_dropped}). It then needs 2 labels, and its rarest label must fill the folds; its test set needs only \code{min_test_genomes}. A drug or class needs 2 of a diff --git a/man/runModelingWorkflow.Rd b/man/runModelingWorkflow.Rd index f3a8ae3..519f5e0 100644 --- a/man/runModelingWorkflow.Rd +++ b/man/runModelingWorkflow.Rd @@ -14,7 +14,7 @@ runModelingWorkflow( n_fold = 5, min_test_genomes = 5, min_test_minority = 2, - min_label_genomes = 25, + min_mdr_label_genomes = 25, mdr_min_classes = 2, overwrite = FALSE, BPPARAM = BiocParallel::SnowParam(workers = 8L, tasks = 128L), @@ -36,7 +36,7 @@ runModelingWorkflow( \item{labels}{\code{"real"}, \code{"shuffled"} or both (the default), passed to \code{\link[=expandTasks]{expandTasks()}}.} -\item{min_genomes, n_fold, min_test_genomes, min_test_minority, min_label_genomes}{Passed to +\item{min_genomes, n_fold, min_test_genomes, min_test_minority, min_mdr_label_genomes}{Passed to \code{\link[=eligibleScopes]{eligibleScopes()}}.} \item{mdr_min_classes}{Passed to \code{\link[=profileORB]{profileORB()}}: fewest drug classes a genome diff --git a/tests/testthat/test-eligibility.R b/tests/testthat/test-eligibility.R index 1ed84ff..8d94ff9 100644 --- a/tests/testthat/test-eligibility.R +++ b/tests/testthat/test-eligibility.R @@ -138,7 +138,7 @@ test_that("vanilla scopes are the profile's whole-dataset counts", { expect_true(all(is.na(scopes$test_n_genomes))) # The MDR combination target counts only the genomes its label filter keeps: here none, - # since its one label has 5 genomes, under min_label_genomes. + # since its one label has 5 genomes, under min_mdr_label_genomes. combination <- scopes$target == "resistant_classes" expect_equal(scopes$n_genomes[!combination], whole$n_genomes[!combination]) expect_equal(scopes$n_resistant[!combination], whole$n_resistant[!combination]) @@ -387,7 +387,7 @@ fxCombinationScopes <- function(eligibility, mode_id) { test_that("a combination scope drops labels with too few training genomes, and their genomes", { eligibility <- fxEligibility( fxCombinations(4, 3, 2), - modes = "vanilla", min_genomes = 4, n_fold = 2, min_label_genomes = 3 + modes = "vanilla", min_genomes = 4, n_fold = 2, min_mdr_label_genomes = 3 ) scope <- fxCombinationScopes(eligibility, "vanilla") @@ -403,7 +403,7 @@ test_that("a combination scope drops labels with too few training genomes, and t test_that("a combination scope with one label left is single_phenotype", { eligibility <- fxEligibility( fxCombinations(4, 1, 1), - modes = "vanilla", min_genomes = 2, n_fold = 2, min_label_genomes = 3 + modes = "vanilla", min_genomes = 2, n_fold = 2, min_mdr_label_genomes = 3 ) scope <- fxCombinationScopes(eligibility, "vanilla") @@ -420,7 +420,7 @@ test_that("a cross test set drops genomes whose label the training set doesn't k eligibility <- fxEligibility( metadata, modes = "cross_year", min_genomes = 4, n_fold = 2, min_test_genomes = 2, - min_label_genomes = 3 + min_mdr_label_genomes = 3 ) scopes <- fxCombinationScopes(eligibility, "cross_year") y1_to_y2 <- scopes[scopes$train_group == "Y1", ] diff --git a/tests/testthat/test-fit-models.R b/tests/testthat/test-fit-models.R index 787ddfb..b83ea80 100644 --- a/tests/testthat/test-fit-models.R +++ b/tests/testthat/test-fit-models.R @@ -29,7 +29,7 @@ fxCombinationFit <- function(n = 20, env = parent.frame()) { data <- fxLabelled(n) eligibility <- fxEligibility( data$metadata, modes = "vanilla", parquets = data$parquets, - min_genomes = n, n_fold = 2, min_label_genomes = n, env = env + min_genomes = n, n_fold = 2, min_mdr_label_genomes = n, env = env ) eligibility$scopes <- eligibility$scopes[eligibility$scopes$target == "resistant_classes", ] eligibility$members <- eligibility$members[eligibility$members$target == "resistant_classes", ] From 3cdb8894e2e830e0cc3512ca6d50c2aa82c8df1e Mon Sep 17 00:00:00 2001 From: Abhirupa Ghosh <100681585+AbhirupaGhosh@users.noreply.github.com> Date: Mon, 26 Oct 2026 16:53:15 -0600 Subject: [PATCH 13/14] Update orb-read.R Add the time and location --- R/orb-read.R | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/R/orb-read.R b/R/orb-read.R index eb3f58f..1f2153f 100644 --- a/R/orb-read.R +++ b/R/orb-read.R @@ -445,7 +445,10 @@ listORBs <- function(root = NULL) { } pick <- utils::menu( - paste0(orbs$dataset_label, " [", orbs$dataset_id, "]"), + paste0( + orbs$dataset_label, " [", orbs$dataset_id, "], finished ", orbs$finished_at, + ", ", orbs$directory + ), title = "Select an ORB to read:" ) From 85a1860c10b6a9d909ea53a44d94fe99645375d8 Mon Sep 17 00:00:00 2001 From: Emily Boyer Date: Tue, 27 Oct 2026 10:53:28 -0600 Subject: [PATCH 14/14] feat: skip MDR combination scopes whose labels span fewer than two classes With one class, the combination labels are that class or Susceptible, so the model repeats the class model. New eligibility rule single_class_labels, which also covers a dataset tested on one class only (where the yes/no MDR target has no model either). Suggested by @AbhirupaGhosh. --- R/eligibility.R | 28 +++++++++++++++++++++++++++- R/errors.R | 2 ++ design/mdr.md | 4 ++++ man/eligibleScopes.Rd | 3 ++- tests/testthat/test-eligibility.R | 21 +++++++++++++++++++++ 5 files changed, 56 insertions(+), 2 deletions(-) diff --git a/R/eligibility.R b/R/eligibility.R index 2dd967a..92819bc 100644 --- a/R/eligibility.R +++ b/R/eligibility.R @@ -242,6 +242,27 @@ ) } +# For each multi-class scope, how many drug classes its kept training genomes are resistant +# to. Fewer than two means its labels are one class or Susceptible: that class's model again. +.labelClassCounts <- function(members, profile) { + phenotypes <- profile$phenotypes + resistant <- phenotypes[ + phenotypes$unit == "drug_class" & phenotypes$phenotype == "Resistant", + c("genome_id", "target") + ] + names(resistant) <- c("genome_id", "class") + kept <- members[ + .isMulticlass(members$unit, members$target) & members$kept & members$role == "train", + c(.SCOPE_KEYS, "genome_id") + ] + classes <- dplyr::inner_join(kept, resistant, by = "genome_id", relationship = "many-to-many") + dplyr::summarise( + dplyr::group_by(classes, dplyr::across(dplyr::all_of(.SCOPE_KEYS))), + n_label_classes = dplyr::n_distinct(.data$class), + .groups = "drop" + ) +} + # Each scope's first failed rule, or NA if it passes them all. .firstFailedRule <- function(scopes, mode, settings, has_test) { # A class with data for only one of its drugs is that drug's model again. @@ -254,6 +275,8 @@ multiclass, scopes$n_rarest, pmin(scopes$n_resistant, scopes$n_susceptible) ) single <- ifelse(multiclass, scopes$n_labels < 2, rarer == 0) + # A scope with no resistant kept genome has no count (NA), so fewer than two too. + one_class <- multiclass & !(!is.na(scopes$n_label_classes) & scopes$n_label_classes >= 2) enough_for_cv <- if (has_test) { rarer >= settings$n_fold } else { @@ -269,6 +292,7 @@ failed <- list( duplicate_class = duplicate, single_phenotype = single, + single_class_labels = one_class, too_few_genomes = scopes$n_genomes < settings$min_genomes, too_few_for_cv = !enough_for_cv, test_too_small = !test_ok @@ -300,6 +324,7 @@ has_test <- mode$LOO || mode$cross_test scopes <- .scopeCounts(built$members, has_test) scopes <- dplyr::left_join(scopes, .classDrugCounts(built$members, profile), by = .SCOPE_KEYS) + scopes <- dplyr::left_join(scopes, .labelClassCounts(built$members, profile), by = .SCOPE_KEYS) scopes$rule_id <- .firstFailedRule(scopes, mode, settings, has_test) decisions <- scopes[c(.SCOPE_KEYS, counts, "rule_id")] @@ -333,7 +358,8 @@ #' (`resistant_classes`), each combination is a label: a scope keeps the labels #' with at least `min_mdr_label_genomes` training genomes and drops the genomes, #' training or test, with any other (`n_dropped`, `test_n_dropped`). It then -#' needs 2 labels, and its rarest label must fill the folds; its test set needs +#' needs 2 labels involving at least 2 drug classes (otherwise it repeats a +#' class model), and its rarest label must fill the folds; its test set needs #' only `min_test_genomes`. A drug or class needs 2 of a #' mode's groups (3 for leave-one-out) to get scopes, and a stratified one needs #' 2 eligible groups. A scope fails, in order, if it's a class with data for only diff --git a/R/errors.R b/R/errors.R index 83b4393..be8656d 100644 --- a/R/errors.R +++ b/R/errors.R @@ -73,6 +73,8 @@ "The class has data for only one of its drugs.", "single_phenotype", "eligibility", "scope", "Training genomes have one phenotype only.", + "single_class_labels", "eligibility", "scope", + "MDR combination labels involve fewer than two drug classes, so it repeats a class model.", "too_few_genomes", "eligibility", "scope", "Too few training genomes.", "too_few_for_cv", "eligibility", "scope", diff --git a/design/mdr.md b/design/mdr.md index 6d4d796..6133d32 100644 --- a/design/mdr.md +++ b/design/mdr.md @@ -115,6 +115,10 @@ New argument `min_mdr_label_genomes = 25` (old `min_n`). Per MDR scope: label with 3 genomes in one year from breaking that year's folds. 2. Then apply the usual rules, with "phenotype" meaning label: - `single_phenotype`: fewer than 2 labels kept. + - `single_class_labels`: the kept labels involve fewer than 2 drug classes + (e.g. only `CLX` and `Susceptible`), so the model repeats that class's model. + This covers a dataset tested on one class only, where the yes/no target + also has no model (every genome is `undetermined`). (Suggested by Abhirupa.) - `too_few_genomes`: fewer than `min_genomes` genomes kept. - `too_few_for_cv`: the rarest kept label can't fill `n_fold` folds after the holdout. With 25 genomes per label this always passes. diff --git a/man/eligibleScopes.Rd b/man/eligibleScopes.Rd index b541ca0..6782224 100644 --- a/man/eligibleScopes.Rd +++ b/man/eligibleScopes.Rd @@ -53,7 +53,8 @@ Counts are of genomes, each with one phenotype. In the MDR combination target (\code{resistant_classes}), each combination is a label: a scope keeps the labels with at least \code{min_mdr_label_genomes} training genomes and drops the genomes, training or test, with any other (\code{n_dropped}, \code{test_n_dropped}). It then -needs 2 labels, and its rarest label must fill the folds; its test set needs +needs 2 labels involving at least 2 drug classes (otherwise it repeats a +class model), and its rarest label must fill the folds; its test set needs only \code{min_test_genomes}. A drug or class needs 2 of a mode's groups (3 for leave-one-out) to get scopes, and a stratified one needs 2 eligible groups. A scope fails, in order, if it's a class with data for only diff --git a/tests/testthat/test-eligibility.R b/tests/testthat/test-eligibility.R index 8d94ff9..e6e398e 100644 --- a/tests/testthat/test-eligibility.R +++ b/tests/testthat/test-eligibility.R @@ -412,6 +412,27 @@ test_that("a combination scope with one label left is single_phenotype", { expect_false(any(eligibility$members$target == "resistant_classes")) }) +test_that("a combination scope whose labels involve one class is skipped", { + # CLX-only and Susceptible genomes: the labels are CLX and Susceptible, the CLX model again. + one_class <- fxEligibility( + fxCombinations(4, 0, 4), + modes = "vanilla", min_genomes = 4, n_fold = 2, min_mdr_label_genomes = 3 + ) + scope <- fxCombinationScopes(one_class, "vanilla") + + expect_equal(scope$n_labels, 2) + expect_equal(scope$rule_id, "single_class_labels") + expect_false(any(one_class$members$target == "resistant_classes")) + + # A whole dataset tested on one class: no combination model, and no yes/no MDR model either. + metadata <- rbind(fxRows(40, 40, drug = "DRA"), fxRows(40, 40, drug = "DRB")) + eligibility <- fxEligibility(metadata, modes = "vanilla", min_mdr_label_genomes = 25) + scopes <- fxScopes(eligibility, "vanilla") + + expect_equal(scopes$rule_id[scopes$target == "resistant_classes"], "single_class_labels") + expect_false("mdr" %in% scopes$target) +}) + test_that("a cross test set drops genomes whose label the training set doesn't keep", { metadata <- rbind( fxCombinations(3, 3, 0, from = 1, year_bin = "Y1"),