feat: MDR targets (yes/no and combination) and design spec

.Rbuildignore

@@ -7,4 +7,5 @@
 ^CODE_OF_CONDUCT\.md$
 ^doc$
 ^Meta$
-^results$
+^results$
+^design$

R/eligibility.R

@@ -122,14 +122,33 @@
   .asMembers(train, test)
 }
 
+# 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_mdr_label_genomes) {
+  members$kept <- TRUE
+  multiclass <- .isMulticlass(members$unit, members$target)
+  if (!any(multiclass)) {
+    return(members)
+  }
+
+  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_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
+}
+
 # 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) {
   counts <- dplyr::summarise(
     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))
@@ -219,26 +242,57 @@
   )
 }
 
+# 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.
   duplicate <- scopes$unit == "drug_class" & scopes$n_drugs %in% 1
 
-  rarer <- pmin(scopes$n_resistant, scopes$n_susceptible)
+  # 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(
+    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 {
     rarer * (1 - settings$holdout) >= settings$n_fold
   }
   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
+      (multiclass |
+        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,
+    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
@@ -263,18 +317,22 @@
 # 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)
   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")]
 
   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 +354,13 @@
 #' 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_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 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
 #' one of its drugs, has one phenotype, has fewer than `min_genomes` genomes,
@@ -314,12 +378,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_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:
 #'   * `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 +401,8 @@ eligibleScopes <- function(profile,
                            min_genomes = 40,
                            n_fold = 5,
                            min_test_genomes = 5,
-                           min_test_minority = 2) {
+                           min_test_minority = 2,
+                           min_mdr_label_genomes = 25) {
   call <- rlang::current_env()
 
   # Check the arguments.
@@ -343,6 +411,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_mdr_label_genomes, "min_mdr_label_genomes", call = call)
   modes <- modes %||% .MODES$mode_id
   .checkArgModes(modes, call = call)
 
@@ -351,6 +420,7 @@ eligibleScopes <- function(profile,
     n_fold = n_fold,
     min_test_genomes = min_test_genomes,
     min_test_minority = min_test_minority,
+    min_mdr_label_genomes = min_mdr_label_genomes,
     holdout = .CV_HOLDOUT
   )
 
@@ -364,13 +434,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_mdr_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 +457,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(),
@@ -423,6 +499,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_mdr_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]

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",

R/fit-execution.R

@@ -14,7 +14,7 @@
           task = task, path = job$paths[[j]],
           rows = job$rows[[j]], columns = job$columns[[j]],
           encoding = job$encoding[[j]], cache = cache,
-          fitted = shared_fit, return_fitted = keep_fit
+          fitted = shared_fit, return_fitted = keep_fit, multi = job$multi[[j]]
         ),
         .rng_kind = "Mersenne-Twister", .rng_normal_kind = "Inversion",
         .rng_sample_kind = "Rejection"
@@ -40,6 +40,7 @@
     rows = index$rows[index$matrix_scopes[positions]],
     columns = index$columns[positions],
     encoding = matrices$encoding[positions],
+    multi = .isMulticlass(matrices$unit[positions], matrices$target[positions]),
     checkpoint = checkpoint
   )
 }
@@ -109,19 +110,20 @@
     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()
     ),
     importance = tibble::tibble(
-      task_id = character(), Variable = character(),
+      task_id = character(), label = character(), Variable = character(),
       Importance = numeric(), Sign = character()
     ),
     predictions = tibble::tibble(
       task_id = character(), genome_id = character(), truth = character(),
-      predicted = character(), prob_resistant = numeric()
+      predicted = character(), prob_resistant = numeric(), prob_predicted = numeric()
     )
   )
 }
@@ -228,7 +230,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:
 #' Invalid arguments or incompatible resume inputs raise `amrml_invalid_argument`.