diff --git a/.Rbuildignore b/.Rbuildignore
index 1ec22b7..ebfabc7 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -7,4 +7,5 @@
 ^CODE_OF_CONDUCT\.md$
 ^doc$
 ^Meta$
-^results$
\ No newline at end of file
+^results$
+^design$
diff --git a/R/eligibility.R b/R/eligibility.R
index 07ddafc..92819bc 100644
--- a/R/eligibility.R
+++ b/R/eligibility.R
@@ -122,6 +122,25 @@
   .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) {
@@ -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))
@@ -219,12 +242,41 @@
   )
 }
 
+# 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 {
@@ -232,13 +284,15 @@
   }
   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]
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/R/fit-execution.R b/R/fit-execution.R
index 89a7141..b1158f5 100644
--- a/R/fit-execution.R
+++ b/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`.
diff --git a/R/fit-models.R b/R/fit-models.R
index aa02828..59d7b15 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 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")
@@ -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 = as.character(ifelse(unname(coefs) > 0, "POS", "NEG"))
+    )
+  }
 
-  tibble::tibble(
-    Variable = features[as.integer(sub("^x", "", names(coefs)))],
-    Importance = abs(unname(coefs)),
-    Sign = as.character(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)
+  }))
 }
 
 # Cross-test tasks can share a fit when the training inputs exactly match
@@ -239,12 +269,18 @@
   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
+  # 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)
   }
@@ -257,23 +293,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)
-  if (is.null(fitted)) fitted <- .fitModel(train, task$n_fold, sparse = use_sparse)
+                     rows$phenotype[split$test], sparse = use_sparse, labels = labels)
+  if (is.null(fitted)) 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,
@@ -289,11 +329,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]]
@@ -326,6 +368,15 @@ print.amr_fit <- function(x, ...) {
   cat("<amr_fit>", 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/R/orb-profile.R b/R/orb-profile.R
index 1425da7..4631b02 100644
--- a/R/orb-profile.R
+++ b/R/orb-profile.R
@@ -99,6 +99,44 @@
   )
 }
 
+# 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 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"))
@@ -145,7 +183,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"
   )
@@ -158,9 +197,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
@@ -169,6 +208,18 @@
 #' if they disagree, profiling stops. Missing feature rows represent zeros in
 #' sparse feature tables.
 #'
+#' 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.
+#'
+#' 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:
 #'
@@ -181,13 +232,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 +255,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 +303,11 @@ 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]
+  combinations <- .classCombinations(class_phenotypes)
   phenotypes <- dplyr::bind_rows(
     tibble::tibble(
       genome_id = drug_phenotypes$genome_id,
@@ -260,6 +320,19 @@ 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])
+    ),
+    tibble::tibble(
+      genome_id = combinations$genome_id,
+      unit = "mdr",
+      target = .MDR_COMBINATION,
+      phenotype = combinations$phenotype
     )
   )
 
@@ -284,6 +357,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 +392,17 @@ 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 = ""
+  )
+  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/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:"
   )
 
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/R/workflow.R b/R/workflow.R
index 0c7e240..90f7ccd 100644
--- a/R/workflow.R
+++ b/R/workflow.R
@@ -22,8 +22,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_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.
 #' @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.
@@ -57,6 +59,8 @@ runModelingWorkflow <- function(orb_dir,
                                 n_fold = 5,
                                 min_test_genomes = 5,
                                 min_test_minority = 2,
+                                min_mdr_label_genomes = 25,
+                                mdr_min_classes = 2,
                                 overwrite = FALSE,
                                 BPPARAM = BiocParallel::SnowParam(workers = 8L, tasks = 128L),
                                 reuse_fits = TRUE,
@@ -68,6 +72,7 @@ runModelingWorkflow <- function(orb_dir,
   .checkArgLabels(labels, call = call)
   .checkArgFlag(overwrite, "overwrite", call = call)
   .checkArgClass(BPPARAM, "BPPARAM", "BiocParallelParam", "BiocParallel", call = call)
+  .checkArgCount(mdr_min_classes, "mdr_min_classes", min = 2, call = call)
   .checkArgFlag(reuse_fits, "reuse_fits", call = call)
   .checkArgFlag(resume, "resume", call = call)
   if (resume && (is.null(out_dir) || overwrite)) {
@@ -79,6 +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_mdr_label_genomes = min_mdr_label_genomes, mdr_min_classes = mdr_min_classes,
     reuse_fits = reuse_fits, feature_types = feature_types
   )
   if (resume) {
@@ -97,8 +103,10 @@ runModelingWorkflow <- function(orb_dir,
     orb <- readORB(orb_dir)
     out_dir <- out_dir %||% .defaultOutDir(orb$directory)
     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
+      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_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
new file mode 100644
index 0000000..6133d32
--- /dev/null
+++ b/design/mdr.md
@@ -0,0 +1,203 @@
+# MDR: modeling combinations of resistant classes
+
+Status: both MDR targets are implemented. See "Checking against old amRml" for results.
+
+## 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
+
+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`:
+
+- `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_mdr_label_genomes = 25` (old `min_n`). Per MDR scope:
+
+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.
+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.
+   - `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 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()`)
+
+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
+
+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.**
+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.**
diff --git a/man/eligibleScopes.Rd b/man/eligibleScopes.Rd
index ecc53c8..6782224 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_mdr_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_mdr_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,13 @@ 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_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 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
 one of its drugs, has one phenotype, has fewer than \code{min_genomes} genomes,
diff --git a/man/fitModels.Rd b/man/fitModels.Rd
index 334d6d9..2358c87 100644
--- a/man/fitModels.Rd
+++ b/man/fitModels.Rd
@@ -34,7 +34,12 @@ 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 tasks from \code{\link[=expandTasks]{expandTasks()}} in parallel with old amRml's method:
diff --git a/man/profileORB.Rd b/man/profileORB.Rd
index fa39240..5205fe2 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}
@@ -24,9 +30,9 @@ for the whole dataset.
 }
 }
 \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
@@ -34,6 +40,18 @@ Susceptible row. 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. Missing feature rows represent zeros in
 sparse feature tables.
+
+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.
+
+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 06d17aa..519f5e0 100644
--- a/man/runModelingWorkflow.Rd
+++ b/man/runModelingWorkflow.Rd
@@ -14,6 +14,8 @@ runModelingWorkflow(
   n_fold = 5,
   min_test_genomes = 5,
   min_test_minority = 2,
+  min_mdr_label_genomes = 25,
+  mdr_min_classes = 2,
   overwrite = FALSE,
   BPPARAM = BiocParallel::SnowParam(workers = 8L, tasks = 128L),
   reuse_fits = TRUE,
@@ -34,9 +36,12 @@ 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}{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
+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.}
diff --git a/tests/testthat/test-eligibility.R b/tests/testthat/test-eligibility.R
index 5702672..e6e398e 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_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])
+  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,105 @@ 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_mdr_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_mdr_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 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"),
+    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_mdr_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 +519,13 @@ 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\\)")
 })
+
+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")
+})
diff --git a/tests/testthat/test-fit-models.R b/tests/testthat/test-fit-models.R
new file mode 100644
index 0000000..b83ea80
--- /dev/null
+++ b/tests/testthat/test-fit-models.R
@@ -0,0 +1,94 @@
+# 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_mdr_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("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))
+
+  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")))
+})
+
+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")
+})
diff --git a/tests/testthat/test-orb-profile.R b/tests/testthat/test-orb-profile.R
index f6bf95d..28b0456 100644
--- a/tests/testthat/test-orb-profile.R
+++ b/tests/testthat/test-orb-profile.R
@@ -207,6 +207,89 @@ 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$target == "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("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())
+
+  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", {
@@ -215,12 +298,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", {