From 1cf4dd094a5a8e4bc94f7c5e1a54bcd5b1d9b287 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Juli=C3=A1n=20D=2E=20Ot=C3=A1lvaro?=
 <juliandavid347@gmail.com>
Date: Tue, 15 Sep 2026 17:55:40 +0100
Subject: [PATCH 1/6] fix: associate PTA results with the requested regimen
 when ids are text

Subject ids are text whenever a data file labels subjects with anything other
than plain integers, and Pmetrics ordered and grouped data by that text. Text
sorts by bytes, so a 12-regimen simulation came out as 1, 10, 11, 12, 2, 3, ...
PM_pta numbers its rows by regimen position, so row 4 described the regimen with
id 12 rather than id 4. The PDI calculation and `prop_success == mean(success)`
were correct throughout; only the association was wrong.

`pm_id_rank()` (new, internal) returns an integer rank per id: numeric ids sort
numerically, anything else keeps byte order, and equal ids share a rank. It
replaces `id` as the sort key at every site that decides regimen order:

- PM_data standardization, where `group_by(id, ...)` orders text ids as text
- PM_sim, where `.id` was ranked with `dense_rank(id)` and two `arrange(.id, ...)`
  calls sorted the id string
- PM_pta, where `dense_rank(id)` reported the same ranking

`dense_rank(id)` in PM_pta came from #341 ("fix: mis-sorting during PTA"), which
converted ids to ranks to fix a sorting problem; the rank of a text id is
lexicographic, so ids of different lengths were not covered. This is an
incomplete earlier fix rather than a later regression.

Ids that are already numeric are unaffected, and ids such as "P01"/"P02" keep
the order they had.

The regression test builds 12 regimens with text ids and log-spaced doses
spanning the target boundary, then compares each multi-regimen PM_pta row with
PM_pta run on that regimen alone, so it tests the association without
reimplementing the PDI calculation. Without this fix it fails with the reported
pattern ("1" "10" "11" "12" "2" ...).

Full testthat suite: no failures.
---
 R/PM_data.R                                | 10 ++-
 R/PM_pta.R                                 |  4 +-
 R/PM_sim.R                                 |  7 +-
 R/id_order.R                               | 28 ++++++
 tests/testthat/test-pm-pta-regimen-order.R | 99 ++++++++++++++++++++++
 5 files changed, 140 insertions(+), 8 deletions(-)
 create mode 100644 R/id_order.R
 create mode 100644 tests/testthat/test-pm-pta-regimen-order.R

diff --git a/R/PM_data.R b/R/PM_data.R
index 10fc4fa6a..55019b47f 100755
--- a/R/PM_data.R
+++ b/R/PM_data.R
@@ -430,20 +430,24 @@ PM_data <- R6::R6Class("PM_data",
       }
 
       # Assign an observation occasion number for each id, incremented at each EVID = 4.
+      # `group_by(id)` orders text ids as text, putting subject 10 before subject
+      # 2; the final arrange restores the order the ids were written in.
       dataObj <- makePMdataOccasion(dataObj) |>
         group_by(id, occasion) |>
         arrange(time, desc(evid), .by_group = TRUE) |>
         ungroup() |>
-        relocate(occasion, .after = last_col())
+        relocate(occasion, .after = last_col()) |>
+        arrange(pm_id_rank(id), occasion, time)
 
       if ("evid" %in% names(dataObj_orig)) {
         dataObj_orig <- makePMdataOccasion(dataObj_orig) |>
           group_by(id, occasion) |>
           arrange(time, desc(evid), .by_group = TRUE) |>
           ungroup() |>
+          arrange(pm_id_rank(id), occasion, time) |>
           select(-occasion)
       } else {
-        dataObj_orig <- dataObj_orig |> arrange(id, time, out)
+        dataObj_orig <- dataObj_orig |> arrange(pm_id_rank(id), time, out)
       }
 
 
@@ -1868,7 +1872,7 @@ plot.PM_data <- function(
   # time after dose
   if (tad) {
     dat$standard_data$time <- calcTAD(dat$standard_data)
-    dat$standard_data <- dat$standard_data |> arrange(id, time)
+    dat$standard_data <- dat$standard_data |> arrange(pm_id_rank(id), time)
   }
 
   # filter
diff --git a/R/PM_pta.R b/R/PM_pta.R
index 4e7d1d3d8..75c52c24d 100755
--- a/R/PM_pta.R
+++ b/R/PM_pta.R
@@ -228,13 +228,13 @@ PM_pta <- R6::R6Class(
             simdata <- simdata$data
           } else { # just one sim
 
-            simdata$data$obs <- simdata$data$obs |> mutate(id = dplyr::dense_rank(id)) # ensure numeric
+            simdata$data$obs <- simdata$data$obs |> mutate(id = pm_id_rank(id)) # ensure numeric
             simdata <- split(simdata$data$obs, as.factor(simdata$data$obs$id))
           }
         }
 
         if (dataType == 1) { # PM_sim_data object
-          simdata$obs <- simdata$obs |> mutate(id = dplyr::dense_rank(id)) # ensure numeric
+          simdata$obs <- simdata$obs |> mutate(id = pm_id_rank(id)) # ensure numeric
           simdata <- split(simdata, as.factor(simdata$id)) # split by id
         }
 
diff --git a/R/PM_sim.R b/R/PM_sim.R
index dee1d6a8a..5a4e42a97 100755
--- a/R/PM_sim.R
+++ b/R/PM_sim.R
@@ -1639,7 +1639,7 @@ PM_sim <- R6::R6Class(
           )
         
         sim_res <- sim_res |>
-        arrange(.id, comp, nsim, time, outeq) |>
+        arrange(pm_id_rank(.id), comp, nsim, time, outeq) |>
         select(-.id)
         
         obs <- sim_res |> filter(comp == min(comp, na.rm = TRUE)) |> # obs are duplicated in every compartment
@@ -1854,8 +1854,9 @@ PM_sim <- R6::R6Class(
         }
         
         # first, add temporary index to ensure id order remains the same
+        # (natural id order: "2" follows "1", not "10")
         dat2 <- template |>
-        mutate(.id = dplyr::dense_rank(id))
+        mutate(.id = pm_id_rank(id))
         
         # second, add predInt if necessary
         if (!is.na(predTimes[1])) {
@@ -1877,7 +1878,7 @@ PM_sim <- R6::R6Class(
           }) |>
           bind_rows()
           new_dat <- bind_rows(dat2, dat3) |>
-          arrange(.id, time, outeq) |>
+          arrange(pm_id_rank(.id), time, outeq) |>
           select(-.id)
         } else { # predInt was not specified
           new_dat <- template # the original data without .id
diff --git a/R/id_order.R b/R/id_order.R
new file mode 100644
index 000000000..f0f87c630
--- /dev/null
+++ b/R/id_order.R
@@ -0,0 +1,28 @@
+# Ordering of subject identifiers
+#
+# Pmetrics identifiers are text whenever a data file labels subjects with
+# anything other than plain integers, and sorting text is byte order: "10"
+# precedes "2". Any ordering that follows id order - simulated regimens, PTA
+# rows - then follows that order instead of the one the data was written in,
+# which silently associates results with the wrong subject.
+#
+# `pm_id_rank()` returns an integer rank per element: numeric ids sort
+# numerically, anything else keeps byte order, and equal ids share a rank so
+# they stay together. Use it wherever data are ordered or grouped by id.
+
+pm_id_rank <- function(id) {
+  id <- as.character(id)
+  unique_ids <- unique(id)
+  if (length(unique_ids) == 0) {
+    return(integer(0))
+  }
+
+  numeric_ids <- suppressWarnings(as.numeric(unique_ids))
+  ordered <- if (all(!is.na(numeric_ids))) {
+    unique_ids[order(numeric_ids)]
+  } else {
+    sort(unique_ids, na.last = TRUE)
+  }
+
+  match(id, ordered)
+}
diff --git a/tests/testthat/test-pm-pta-regimen-order.R b/tests/testthat/test-pm-pta-regimen-order.R
new file mode 100644
index 000000000..5768aa491
--- /dev/null
+++ b/tests/testthat/test-pm-pta-regimen-order.R
@@ -0,0 +1,99 @@
+testthat::skip_on_cran()
+
+# Regression: PTA rows must follow the simulation ids a user supplied.
+#
+# When ids are text ("1", "2", ... "12"), everything that ordered data by id
+# sorted it as text, so the regimens came out as 1, 10, 11, 12, 2, 3, ...
+# PM_pta numbers its rows by regimen position, so row 4 described the regimen
+# with id 12 rather than id 4. Nothing in the calculation was wrong; the results
+# were attached to the wrong regimen.
+
+skip_if_not_pmcore <- function() {
+  testthat::skip_if_not_installed("Pmetrics")
+}
+
+pta_fixture <- function(n_ids = 12, nsim = 4) {
+  manual_poppar <- list(
+    wt = 1,
+    mean = as.list(NPex$final$popMean),
+    cov = diag(diag(NPex$final$popCov))
+  )
+  manual_limits <- tibble::tibble(
+    par = c("ka", "ke", "v", "tlag1"),
+    min = c(0.2, 0.01, 40, 0),
+    max = c(0.8, 0.09, 110, 3)
+  )
+
+  # Log-spaced dose scales put the regimens across the target boundary, so each
+  # one has a different response rather than all succeeding or all failing.
+  scales <- exp(seq(log(0.15), log(0.45), length.out = n_ids))
+  template_data <- dataEx$standard_data |>
+    dplyr::filter(id %in% seq_len(n_ids)) |>
+    dplyr::mutate(
+      id = as.character(id),
+      dose = dose * scales[as.integer(id)]
+    )
+  template <- PM_data$new(as.data.frame(template_data), quiet = TRUE)
+
+  sim_for <- function(include) {
+    suppressMessages(PM_sim$new(
+      poppar = manual_poppar,
+      data = template,
+      model = NPex$model,
+      include = include,
+      nsim = nsim,
+      predInt = 5,
+      limits = manual_limits
+    ))
+  }
+
+  list(
+    ids = as.character(seq_len(n_ids)),
+    sim_all = sim_for(seq_len(n_ids)),
+    sim_for = sim_for
+  )
+}
+
+pta_pdi <- function(sim, mic = 2) {
+  suppressMessages(PM_pta$new(
+    simdata = sim, target = mic, target_type = "time", success = 0.5,
+    outeq = 1, free_fraction = 1, start = 0, end = Inf
+  ))$data$data
+}
+
+testthat::test_that("simulated regimens follow numeric id order, not text order", {
+  fx <- pta_fixture()
+  observed <- unique(fx$sim_all$data$obs$id)
+
+  testthat::expect_equal(observed, fx$ids)
+  # The text order that caused the mis-association.
+  testthat::expect_false(identical(observed, sort(fx$ids)))
+})
+
+testthat::test_that("each PTA row describes the regimen with the same id", {
+  fx <- pta_fixture()
+  rows <- pta_pdi(fx$sim_all)
+  testthat::expect_equal(rows$reg_num, seq_along(fx$ids))
+
+  # PM_pta on a single regimen gives that regimen's own result, so comparing the
+  # multi-regimen rows against single-regimen runs tests the association without
+  # reimplementing the PDI calculation.
+  single_pdi <- lapply(fx$ids, function(id) {
+    unlist(pta_pdi(fx$sim_for(id))$pdi[[1]])
+  })
+
+  fingerprints <- vapply(seq_along(fx$ids), function(i) {
+    paste(round(unlist(rows$pdi[[i]]), 8), collapse = "|")
+  }, character(1))
+  # Guard against a vacuous pass: if every regimen returned the same vector the
+  # comparison below would hold for any ordering.
+  testthat::expect_gte(length(unique(fingerprints)), 4)
+
+  for (i in seq_along(fx$ids)) {
+    testthat::expect_equal(
+      round(unlist(rows$pdi[[i]]), 8),
+      round(single_pdi[[i]], 8),
+      info = sprintf("reg_num %d should describe id %s", i, fx$ids[i])
+    )
+  }
+})

From e64770688d04d3f35d9b108393dc914d27474ca6 Mon Sep 17 00:00:00 2001
From: Markus <66058642+mhovd@users.noreply.github.com>
Date: Wed, 16 Sep 2026 12:27:08 +0200
Subject: [PATCH 2/6] fix: keep the regimen id with each PTA result

PM_pta identified a regimen by its position in the list of simulations: the id
was replaced by a rank in `initialize()` and the rows were numbered by the loop
index in `make()`. Nothing downstream could tell which regimen a result
belonged to, so whenever the simulations came back in a different order than
the data the rows described the wrong regimen. That happens for any id that is
not a plain integer, because sorting text is byte order and "10" precedes "2".

The ranking added in #391 orders the results so that this is less likely, but
it does not restore the association: ids such as P1..P12 keep byte order, and
a rank is not an identifier. This keeps the identifier with the result instead:

- `initialize()` splits the simulated observations by id while preserving the
  order the ids first appear, and keeps the id itself as the name of each
  element rather than replacing it with a rank. This also fixes the
  `PM_sim_data` path (`PM_pta$new(simEx$data)`), which split the whole object
  instead of its observations and errored.
- `make()` carries the id of each regimen into every result tibble, including
  `intersect`, and `summary()` reports it. `reg_num` stays the position of the
  regimen in `simdata` - the order in which its id first appears - which is
  what `start`, `end`, and `simlabels` are keyed on, and for user-supplied
  lists that repeat an id it stays positional.
- Results saved before the id column existed fall back to the regimen number,
  which is the only identifier those rows were ever keyed on, so old objects
  still summarise and plot.

The ranking is left in place for ordering: the order of the regimens is still
user-visible through `reg_num` and `label`, and the order of the data is what
PM_data and PM_sim already report. Nothing depends on it for correctness now.

The regression test is kept and extended. It asserts the id column for
numeric, character, and prefixed ids; compares each row with a single-regimen
run for the id it reports, so the association is checked without depending on
the order; and covers a shuffled dose assignment, a subset of regimens, a
user-supplied list of regimens, and a saved object without the id column. The
new assertions fail on both the base commit and on #391.

Not changed, but worth knowing: the `postToUse[i]`/`toInclude[i]` pairing in
PM_sim's posterior branch joins posteriors to subjects by position, so a
posterior file listing ids in a different order than the data still attaches
the wrong posterior. PM_pta's new id column does not fix that, but it makes
such a mix-up visible in the results. The single-regimen paths for PM_post and
PM_data (dataType 5-8) are broken before this change - `simdata` is nested one
level too deep for `make()` to read it - and are left as they are.

Full testthat suite: no failures.
---
 R/PM_pta.R                                 |  90 +++++++--
 man/PM_pta.Rd                              |   9 +-
 man/summary.PM_pta.Rd                      |   3 +-
 tests/testthat/test-pm-pta-regimen-order.R | 205 ++++++++++++++++-----
 4 files changed, 249 insertions(+), 58 deletions(-)

diff --git a/R/PM_pta.R b/R/PM_pta.R
index 75c52c24d..5fdbbe39e 100755
--- a/R/PM_pta.R
+++ b/R/PM_pta.R
@@ -23,6 +23,12 @@
 #' Both methods require the prior creation of a simulation of
 #' appropriate regimens.
 #'
+#' Every result is labeled with the `id` of the regimen it was calculated from,
+#' so results can be associated with the requested regimens without depending on
+#' the order in which the regimens were simulated. `reg_num` and the `simlabels`
+#' annotation are positional: they follow the order of the regimens in
+#' `simdata`, which for a [PM_sim] is the order of the ids in the data.
+#'
 #' @author Julian Otalvaro and Michael Neely
 #' @export
 PM_pta <- R6::R6Class(
@@ -102,7 +108,8 @@ PM_pta <- R6::R6Class(
     #' or `NA` if only one `target_type` was specified.
     #' The individual elements are tibbles with all possible combinations
     #' of `target`s and simulated regimens for a given `target_type`. The tibbles have the following columns:
-    #' * **reg_num** The simulation number in `simdata`.
+    #' * **reg_num** The position of the simulation in `simdata`, i.e. the order in which its id first appears.
+    #' * **id** The identifier of the simulated regimen, as it appears in the `id` column of the simulation.
     #' * **label** Annotation of the simulation, supplied by the `simlabels` argument.
     #' * **target** is the specified `target` for the results row. If a distribution created by [makePTAtarget],
     #' this will be a tibble with  the simulated targets
@@ -227,15 +234,12 @@ PM_pta <- R6::R6Class(
           if (inherits(simdata$data, "PM_simlist")) { # multiple sims
             simdata <- simdata$data
           } else { # just one sim
-
-            simdata$data$obs <- simdata$data$obs |> mutate(id = pm_id_rank(id)) # ensure numeric
-            simdata <- split(simdata$data$obs, as.factor(simdata$data$obs$id))
+            simdata <- split_by_regimen(simdata$data$obs) # split by id
           }
         }
 
         if (dataType == 1) { # PM_sim_data object
-          simdata$obs <- simdata$obs |> mutate(id = pm_id_rank(id)) # ensure numeric
-          simdata <- split(simdata, as.factor(simdata$id)) # split by id
+          simdata <- split_by_regimen(simdata$obs) # split by id
         }
 
         if (dataType == 2) { # PM_simlist
@@ -291,7 +295,7 @@ PM_pta <- R6::R6Class(
         )
         private$populate(pta)
       } else { # try simdata as a filename
-        pta <- PM_upgrade(readRDS(simdata))
+        pta <- backfill_pta_id(PM_upgrade(readRDS(simdata)))
         private$populate(pta)
       }
     },
@@ -333,6 +337,18 @@ PM_pta <- R6::R6Class(
 
       # define some global variables
       n_reg <- length(simdata) # number of regimens
+      # The identifier of each regimen travels with its simulated observations,
+      # so a result can be associated with the regimen it was calculated from
+      # regardless of the order the simulations were returned in. `reg_num`
+      # remains the position of the regimen in `simdata`, i.e. the order in
+      # which its id first appears, and is what `start` and `end` are recycled
+      # by; `id` is the identifier itself and is authoritative.
+      reg_ids <- purrr::map_chr(simdata, regimen_id)
+      if (anyNA(reg_ids) && length(names(simdata)) == n_reg) {
+        # a user-supplied list may name its elements with the regimen ids
+        reg_ids[is.na(reg_ids)] <- names(simdata)[is.na(reg_ids)]
+      }
+      reg_ids <- unname(reg_ids)
       n_sim <- sapply(simdata, function(x) length(unique(x$nsim))) # number of sims per regimen
       n_type <- length(target_type)
       n_success <- length(success)
@@ -469,6 +485,7 @@ PM_pta <- R6::R6Class(
       master_pta <- purrr::map(1:n_reg, \(x){
         purrr::map(1:n_type, \(y) tidyr::expand_grid(
           reg_num = x,
+          id = reg_ids[x],
           target = if (simTarg) { # simulated targets
             if (y == 1) {
               list(
@@ -521,7 +538,7 @@ PM_pta <- R6::R6Class(
         mutate(prop_success = sum(success) / length(success)) |>
         mutate(label = sim_labels[reg_num]) |>
         dplyr::relocate(
-          reg_num, label, target, type,
+          reg_num, id, label, target, type,
           success_ratio, prop_success, success, pdi,
           start, end
         ) |>
@@ -538,7 +555,7 @@ PM_pta <- R6::R6Class(
           x
         })
         names(master_pta) <- NULL
-        master_pta$intersect <- master_pta[[1]] |> select(reg_num, target, success_1 = success, label) # get primary success
+        master_pta$intersect <- master_pta[[1]] |> select(reg_num, id, target, success_1 = success, label) # get primary success
         for (i in 2:n_type) { # add additional success
           master_pta$intersect[[paste0("success_", i)]] <- master_pta[[i]]$success[match(master_pta[[1]]$reg_num, master_pta[[i]]$reg_num)]
         }
@@ -554,7 +571,7 @@ PM_pta <- R6::R6Class(
           mutate(success_ratio = paste0("(", success, ")", collapse = "")) |>
           mutate(success = list(total_success)) |>
           select(
-            reg_num, label, target, type,
+            reg_num, id, label, target, type,
             success_ratio, prop_success, success
           ) |>
           ungroup()
@@ -586,6 +603,50 @@ PM_pta$load <- function(file_name = "PMpta.rds") {
 
 # ACCESSORY INTERNAL FUNCTIONS --------------------------------------------
 
+# Split simulated observations into one element per regimen.
+#
+# `split()` sorts the levels it is given, and for text ids that is byte order
+# ("10" before "2"), which silently reorders the regimens. Building the levels
+# from the order the ids first appear keeps the regimens in the order they were
+# simulated. The id itself - not a rank - becomes the name of the element, so
+# each result can be traced back to the regimen that produced it.
+split_by_regimen <- function(obs) {
+  ids <- as.character(obs$id)
+  split(obs, factor(ids, levels = unique(ids)))
+}
+
+# The identifier of one regimen in `simdata`, i.e. the first id of the
+# observations it contains. `split_by_regimen()` names each element with its id,
+# while user-supplied lists carry the id in the element itself.
+regimen_id <- function(x) {
+  ids <- unique(as.character(x[["id"]]))
+  if (length(ids) == 0) NA_character_ else unname(ids[1])
+}
+
+# PTA results saved before every row carried an `id` column cannot recover the
+# regimen identifiers - they were replaced by a rank when the object was made -
+# so the regimen number, which is what those rows were keyed on, is the best
+# identifier available. Accepts a `PM_pta_data` list or a single result tibble.
+backfill_pta_id <- function(pta) {
+  add_id <- function(x) {
+    if (is.data.frame(x) && "reg_num" %in% names(x) && !"id" %in% names(x)) {
+      x$id <- as.character(x$reg_num)
+      x <- dplyr::relocate(x, id, .after = reg_num)
+    }
+    x
+  }
+
+  if (inherits(pta, "PM_pta_data")) {
+    # the elements are not always named, e.g. one unnamed tibble per target type
+    for (i in seq_along(pta)) {
+      pta[[i]] <- add_id(pta[[i]])
+    }
+    return(pta)
+  }
+
+  add_id(pta)
+}
+
 pta_auc <- function(sims, .target, .simTarg, .start, .end, .pb) {
   .pb$tick()
   auc <- rlang::try_fetch(make_AUC(sims, out ~ time | nsim, start = .start, end = .end),
@@ -1209,9 +1270,10 @@ plot.PM_pta <- function(
 #' intersection of all the PTAs.
 #' @param ci Width of the interval for pharmacodynamic index reporting.  Default is 0.95, i.e. 2.5th to 97.5th percentile.
 #' @param ... Not used.
-#' @return A tibble with the following columns (only the first five if `at = "intersect"`):
+#' @return A tibble with the following columns (only the first six if `at = "intersect"`):
 #'
 #' * **reg_num** is the number of the simulation regimen
+#' * **id** is the identifier of the simulation regimen
 #' * **label** is the simulation label, for reference
 #' * **target** is the target for the row, if targets are discrete, not used for simulated targets
 #' * **type** is the target type for the row, e.g. "auc", "time", "-min", etc.
@@ -1234,7 +1296,7 @@ plot.PM_pta <- function(
 #' @export
 
 summary.PM_pta <- function(object, at = "intersect", ci = 0.95, ...) {
-  pta <- object$clone()$data
+  pta <- backfill_pta_id(object$clone()$data)
 
   if (at == "intersect") {
     if (length(pta$intersect) > 1) {
@@ -1269,7 +1331,7 @@ summary.PM_pta <- function(object, at = "intersect", ci = 0.95, ...) {
 
   if (at != "intersect") {
     pdi <- pdi |> mutate(
-      label = label, type = type, prop_success = prop_success,
+      label = label, id = id, type = type, prop_success = prop_success,
       median = median(unlist(pdi), na.rm = TRUE),
       lower = quantile(unlist(pdi), probs = 0.5 - ci / 2, na.rm = TRUE),
       upper = quantile(unlist(pdi), probs = 0.5 + ci / 2, na.rm = TRUE),
@@ -1280,7 +1342,7 @@ summary.PM_pta <- function(object, at = "intersect", ci = 0.95, ...) {
       .keep = "none"
     )
   } else {
-    pdi <- pdi |> mutate(type = type, label = label, prop_success = prop_success, .keep = "none")
+    pdi <- pdi |> mutate(id = id, type = type, label = label, prop_success = prop_success, .keep = "none")
   }
 
   pdi <- pdi |> ungroup()
diff --git a/man/PM_pta.Rd b/man/PM_pta.Rd
index b2b6a64d3..16b27341d 100755
--- a/man/PM_pta.Rd
+++ b/man/PM_pta.Rd
@@ -20,6 +20,12 @@ and creates the PTA. It is described here.
 
 Both methods require the prior creation of a simulation of
 appropriate regimens.
+
+Every result is labeled with the \code{id} of the regimen it was calculated from,
+so results can be associated with the requested regimens without depending on
+the order in which the regimens were simulated. \code{reg_num} and the \code{simlabels}
+annotation are positional: they follow the order of the regimens in
+\code{simdata}, which for a \link{PM_sim} is the order of the ids in the data.
 }
 \examples{
 
@@ -193,7 +199,8 @@ or \code{NA} if only one \code{target_type} was specified.
 The individual elements are tibbles with all possible combinations
 of \code{target}s and simulated regimens for a given \code{target_type}. The tibbles have the following columns:
 \itemize{
-\item \strong{reg_num} The simulation number in \code{simdata}.
+\item \strong{reg_num} The position of the simulation in \code{simdata}, i.e. the order in which its id first appears.
+\item \strong{id} The identifier of the simulated regimen, as it appears in the \code{id} column of the simulation.
 \item \strong{label} Annotation of the simulation, supplied by the \code{simlabels} argument.
 \item \strong{target} is the specified \code{target} for the results row. If a distribution created by \link{makePTAtarget},
 this will be a tibble with  the simulated targets
diff --git a/man/summary.PM_pta.Rd b/man/summary.PM_pta.Rd
index 05d95b62e..34406ef52 100755
--- a/man/summary.PM_pta.Rd
+++ b/man/summary.PM_pta.Rd
@@ -21,9 +21,10 @@ intersection of all the PTAs.}
 \item{...}{Not used.}
 }
 \value{
-A tibble with the following columns (only the first five if \code{at = "intersect"}):
+A tibble with the following columns (only the first six if \code{at = "intersect"}):
 \itemize{
 \item \strong{reg_num} is the number of the simulation regimen
+\item \strong{id} is the identifier of the simulation regimen
 \item \strong{label} is the simulation label, for reference
 \item \strong{target} is the target for the row, if targets are discrete, not used for simulated targets
 \item \strong{type} is the target type for the row, e.g. "auc", "time", "-min", etc.
diff --git a/tests/testthat/test-pm-pta-regimen-order.R b/tests/testthat/test-pm-pta-regimen-order.R
index 5768aa491..f8392ea90 100644
--- a/tests/testthat/test-pm-pta-regimen-order.R
+++ b/tests/testthat/test-pm-pta-regimen-order.R
@@ -1,18 +1,23 @@
 testthat::skip_on_cran()
 
-# Regression: PTA rows must follow the simulation ids a user supplied.
+# Regression: PTA rows must be associated with the regimen that produced them.
 #
-# When ids are text ("1", "2", ... "12"), everything that ordered data by id
-# sorted it as text, so the regimens came out as 1, 10, 11, 12, 2, 3, ...
-# PM_pta numbers its rows by regimen position, so row 4 described the regimen
-# with id 12 rather than id 4. Nothing in the calculation was wrong; the results
-# were attached to the wrong regimen.
-
-skip_if_not_pmcore <- function() {
-  testthat::skip_if_not_installed("Pmetrics")
-}
+# PM_pta used to discard the regimen id and keep only the position of the
+# regimen in the list of simulations, so a result could be matched to a regimen
+# only by its position. Whenever the simulations came back in a different order
+# than the data - which happens for any id that is not a plain integer, because
+# sorting text is byte order and "10" precedes "2" - the rows described the
+# wrong regimen. Every row now carries the `id` of the regimen it was calculated
+# from, so the association survives any order.
+#
+# Note: the posterior branch of PM_sim cannot be exercised from the shipped
+# data, because `NPex$postMean` is empty, so all regimens here are simulated
+# from the population parameters.
 
-pta_fixture <- function(n_ids = 12, nsim = 4) {
+pta_fixture <- function(n_ids = 12, nsim = 4,
+                        shape = c("numeric", "character", "prefixed"),
+                        scale_order = seq_len(n_ids)) {
+  shape <- match.arg(shape)
   manual_poppar <- list(
     wt = 1,
     mean = as.list(NPex$final$popMean),
@@ -24,14 +29,24 @@ pta_fixture <- function(n_ids = 12, nsim = 4) {
     max = c(0.8, 0.09, 110, 3)
   )
 
+  ids <- switch(shape,
+    numeric = seq_len(n_ids),
+    character = as.character(seq_len(n_ids)),
+    prefixed = paste0("P", seq_len(n_ids))
+  )
+
   # Log-spaced dose scales put the regimens across the target boundary, so each
   # one has a different response rather than all succeeding or all failing.
-  scales <- exp(seq(log(0.15), log(0.45), length.out = n_ids))
+  # `scale_order` hands out the scales out of order, decoupling the response of
+  # a regimen from its position so that only its id identifies it.
+  scales <- exp(seq(log(0.15), log(0.45), length.out = n_ids))[scale_order]
   template_data <- dataEx$standard_data |>
     dplyr::filter(id %in% seq_len(n_ids)) |>
+    # `dose` is scaled before `id` is relabeled, since it is indexed by the
+    # subject number.
     dplyr::mutate(
-      id = as.character(id),
-      dose = dose * scales[as.integer(id)]
+      dose = dose * scales[as.integer(id)],
+      id = ids[as.integer(id)]
     )
   template <- PM_data$new(as.data.frame(template_data), quiet = TRUE)
 
@@ -48,8 +63,10 @@ pta_fixture <- function(n_ids = 12, nsim = 4) {
   }
 
   list(
-    ids = as.character(seq_len(n_ids)),
-    sim_all = sim_for(seq_len(n_ids)),
+    ids = ids, # the ids as the user supplied them
+    data_ids = as.character(unique(template$standard_data$id)), # as PM_data orders them
+    template = template,
+    sim_all = sim_for(ids),
     sim_for = sim_for
   )
 }
@@ -61,39 +78,143 @@ pta_pdi <- function(sim, mic = 2) {
   ))$data$data
 }
 
+# PDI of a single regimen, keyed by the regimen id.
+single_pdi <- function(fx, id) {
+  unlist(pta_pdi(fx$sim_for(id))$pdi[[1]])
+}
+
+pdi_fingerprint <- function(x) paste(round(unlist(x), 8), collapse = "|")
+
 testthat::test_that("simulated regimens follow numeric id order, not text order", {
-  fx <- pta_fixture()
-  observed <- unique(fx$sim_all$data$obs$id)
+  for (shape in c("numeric", "character")) {
+    fx <- pta_fixture(shape = shape)
+    observed <- as.character(unique(fx$sim_all$data$obs$id))
 
-  testthat::expect_equal(observed, fx$ids)
-  # The text order that caused the mis-association.
-  testthat::expect_false(identical(observed, sort(fx$ids)))
+    testthat::expect_equal(observed, as.character(fx$ids), info = shape)
+    # The text order that caused the mis-association.
+    testthat::expect_false(identical(observed, sort(observed)), info = shape)
+  }
 })
 
-testthat::test_that("each PTA row describes the regimen with the same id", {
-  fx <- pta_fixture()
+testthat::test_that("simulated regimens follow the id order of the data", {
+  for (shape in c("numeric", "character", "prefixed")) {
+    fx <- pta_fixture(shape = shape)
+
+    testthat::expect_equal(
+      as.character(unique(fx$sim_all$data$obs$id)),
+      fx$data_ids,
+      info = shape
+    )
+  }
+})
+
+testthat::test_that("each PTA row carries the id of its own regimen", {
+  for (shape in c("numeric", "character", "prefixed")) {
+    fx <- pta_fixture(shape = shape)
+    rows <- pta_pdi(fx$sim_all)
+    sim_ids <- as.character(unique(fx$sim_all$data$obs$id))
+
+    # The id is the regimen identifier itself, not a position and not a rank, so
+    # it names the same regimen however the ids are written.
+    testthat::expect_false(anyNA(rows$id), info = shape)
+    testthat::expect_equal(rows$id, sim_ids[rows$reg_num], info = shape)
+    testthat::expect_equal(rows$id, fx$data_ids[rows$reg_num], info = shape)
+
+    # PM_pta on a single regimen gives that regimen's own result, so comparing
+    # each multi-regimen row with the single-regimen run for the id it reports
+    # tests the association without reimplementing the PDI calculation.
+    fingerprints <- vapply(rows$pdi, pdi_fingerprint, character(1))
+    # Guard against a vacuous pass: if every regimen returned the same vector
+    # the comparison below would hold for any association.
+    testthat::expect_gte(length(unique(fingerprints)), 4)
+
+    for (i in seq_len(nrow(rows))) {
+      testthat::expect_equal(
+        round(unlist(rows$pdi[[i]]), 8),
+        round(single_pdi(fx, rows$id[i]), 8),
+        info = sprintf("%s: row %d reports id %s", shape, i, rows$id[i])
+      )
+    }
+  }
+})
+
+testthat::test_that("the association does not depend on the order of the regimens", {
+  # The dose scales are shuffled, so the response of each regimen is not in the
+  # same order as its id: a row can only match its own regimen by id.
+  fx <- pta_fixture(shape = "character", scale_order = c(5, 11, 1, 8, 12, 3, 9, 2, 7, 4, 10, 6))
   rows <- pta_pdi(fx$sim_all)
-  testthat::expect_equal(rows$reg_num, seq_along(fx$ids))
-
-  # PM_pta on a single regimen gives that regimen's own result, so comparing the
-  # multi-regimen rows against single-regimen runs tests the association without
-  # reimplementing the PDI calculation.
-  single_pdi <- lapply(fx$ids, function(id) {
-    unlist(pta_pdi(fx$sim_for(id))$pdi[[1]])
-  })
-
-  fingerprints <- vapply(seq_along(fx$ids), function(i) {
-    paste(round(unlist(rows$pdi[[i]]), 8), collapse = "|")
-  }, character(1))
-  # Guard against a vacuous pass: if every regimen returned the same vector the
-  # comparison below would hold for any ordering.
-  testthat::expect_gte(length(unique(fingerprints)), 4)
-
-  for (i in seq_along(fx$ids)) {
+  sim_ids <- as.character(unique(fx$sim_all$data$obs$id))
+
+  testthat::expect_equal(rows$id, sim_ids[rows$reg_num])
+  # Guard against a vacuous pass, as above.
+  testthat::expect_gte(length(unique(vapply(rows$pdi, pdi_fingerprint, character(1)))), 4)
+
+  for (i in seq_len(nrow(rows))) {
     testthat::expect_equal(
       round(unlist(rows$pdi[[i]]), 8),
-      round(single_pdi[[i]], 8),
-      info = sprintf("reg_num %d should describe id %s", i, fx$ids[i])
+      round(single_pdi(fx, rows$id[i]), 8),
+      info = sprintf("row %d reports id %s", i, rows$id[i])
     )
   }
 })
+
+testthat::test_that("a subset of the regimens keeps each id with its own result", {
+  fx <- pta_fixture(shape = "prefixed")
+  keep <- c("P7", "P2", "P11", "P4") # out of id order on purpose
+  rows <- pta_pdi(fx$sim_for(keep))
+
+  testthat::expect_setequal(rows$id, keep)
+  testthat::expect_equal(rows$id, fx$data_ids[fx$data_ids %in% keep][rows$reg_num])
+
+  for (i in seq_len(nrow(rows))) {
+    testthat::expect_equal(
+      round(unlist(rows$pdi[[i]]), 8),
+      round(single_pdi(fx, rows$id[i]), 8),
+      info = sprintf("row %d reports id %s", i, rows$id[i])
+    )
+  }
+})
+
+testthat::test_that("a user-supplied list of regimens is associated by its id", {
+  fx <- pta_fixture(shape = "prefixed")
+  obs <- fx$sim_all$data$obs
+  # Deliberately not in id order, and unnamed so that the ids can only come from
+  # the regimens themselves.
+  order <- c("P5", "P1", "P12", "P7")
+  regimens <- unname(split(obs, as.factor(obs$id))[order])
+  rows <- pta_pdi(regimens)
+
+  testthat::expect_equal(rows$id, order)
+  testthat::expect_equal(rows$reg_num, seq_along(order))
+
+  for (i in seq_along(order)) {
+    testthat::expect_equal(
+      round(unlist(rows$pdi[[i]]), 8),
+      round(single_pdi(fx, order[i]), 8),
+      info = sprintf("row %d reports id %s", i, order[i])
+    )
+  }
+})
+
+testthat::test_that("results saved before the id column still summarise and plot", {
+  pta <- suppressMessages(PM_pta$new(
+    simdata = simEx, target = list(5, 10), target_type = c("min", "max"),
+    success = c(1, 1), start = 120, end = 144
+  ))
+  # Emulate a PTA made before results carried the regimen id. Its identifiers
+  # were replaced by a rank, so the regimen number is the only one left.
+  legacy <- pta$clone(deep = TRUE)
+  for (i in seq_along(legacy$data)) {
+    x <- legacy$data[[i]]
+    if (is.data.frame(x) && "id" %in% names(x)) {
+      legacy$data[[i]] <- x[, setdiff(names(x), "id"), drop = FALSE]
+    }
+  }
+  testthat::expect_false("id" %in% names(legacy$data$data))
+
+  summary <- suppressMessages(legacy$summary())
+  testthat::expect_equal(summary$id, as.character(summary$reg_num))
+  single <- suppressMessages(legacy$summary(at = 1))
+  testthat::expect_equal(single$id, as.character(single$reg_num))
+  testthat::expect_no_error(suppressMessages(legacy$plot(print = FALSE)))
+})

From ec36ca7cfe0d93095f26c3241926252a7e64db51 Mon Sep 17 00:00:00 2001
From: Markus <66058642+mhovd@users.noreply.github.com>
Date: Wed, 16 Sep 2026 13:14:06 +0200
Subject: [PATCH 3/6] fix: label a PTA regimen with its id by default

A PTA result is now labeled with the id of the regimen it was calculated from,
instead of the generic "Regimen n". The regimen is already reported in the `id`
column, so the default label repeats the identifier rather than a position, and
a plot legend or summary reads as the regimens the user asked for. Explicit
`simlabels` are unchanged: they are still assigned positionally in `reg_num`
order, and both length warnings are untouched.

A regimen whose id cannot be recovered - an unnamed list of objects that carry
no id column - keeps the generic "Regimen n" text rather than becoming NA,
because `plot.PM_pta` builds its legend and axis labels from `unique(label)` and
a missing label would break plotting.

Docs for `simlabels`, the `label` field, and the plot legend are updated, and
`man/` is regenerated.

The regression test gains the label policy in both directions: the default
equals the id for numeric, character, and prefixed ids and carries through
`summary()`; supplied labels follow `reg_num` while the id still identifies the
regimen; and an unrecoverable id falls back to the generic text and still plots.
The legacy-object test now also resets the labels, so it emulates an object saved
before either change.

Full testthat suite: no failures.
---
 R/PM_pta.R                                 | 21 ++++---
 man/PM_pta.Rd                              | 11 ++--
 man/plot.PM_pta.Rd                         |  3 +-
 tests/testthat/test-pm-pta-regimen-order.R | 71 ++++++++++++++++++++--
 4 files changed, 89 insertions(+), 17 deletions(-)

diff --git a/R/PM_pta.R b/R/PM_pta.R
index 5fdbbe39e..bb97762d4 100755
--- a/R/PM_pta.R
+++ b/R/PM_pta.R
@@ -25,9 +25,10 @@
 #'
 #' Every result is labeled with the `id` of the regimen it was calculated from,
 #' so results can be associated with the requested regimens without depending on
-#' the order in which the regimens were simulated. `reg_num` and the `simlabels`
-#' annotation are positional: they follow the order of the regimens in
-#' `simdata`, which for a [PM_sim] is the order of the ids in the data.
+#' the order in which the regimens were simulated. By default that id is also the
+#' `label`. `reg_num` is positional: it is the order in which the id first appears
+#' in `simdata`, which for a [PM_sim] is the order of the ids in the data, and
+#' explicit `simlabels` are assigned in that order.
 #'
 #' @author Julian Otalvaro and Michael Neely
 #' @export
@@ -55,7 +56,7 @@ PM_pta <- R6::R6Class(
     #' * `NPex$data$standard_data` *PM_data_data*
     #' * `"simout.txt"` matches files with wildcard ability, see [PM_sim]
     #'
-    #' @param simlabels Optional character vector of labels for each simulation.  Default is `c('Regimen 1', 'Regimen 2',...)`.
+    #' @param simlabels Optional character vector of labels for each simulation.  Default is the `id` of each regimen, falling back to `c('Regimen 1', 'Regimen 2',...)` for a regimen whose id is not recoverable. Labels are assigned in the order of `reg_num`, so the first label is for the first regimen.
     #' @param target One of several options.
     #'
     #' * A vector of pharmacodynamic targets, such as Minimum Inhibitory Concentrations (MICs), e.g. `c(0.25, 0.5, 1, 2, 4, 8, 16, 32)`.
@@ -110,7 +111,7 @@ PM_pta <- R6::R6Class(
     #' of `target`s and simulated regimens for a given `target_type`. The tibbles have the following columns:
     #' * **reg_num** The position of the simulation in `simdata`, i.e. the order in which its id first appears.
     #' * **id** The identifier of the simulated regimen, as it appears in the `id` column of the simulation.
-    #' * **label** Annotation of the simulation, supplied by the `simlabels` argument.
+    #' * **label** Annotation of the simulation, supplied by the `simlabels` argument, or the regimen `id` by default.
     #' * **target** is the specified `target` for the results row. If a distribution created by [makePTAtarget],
     #' this will be a tibble with  the simulated targets
     #' * **type** is the specified `target_type` for the results row
@@ -459,7 +460,12 @@ PM_pta <- R6::R6Class(
 
 
       # Check the simulation labels
-      sim_labels <- paste("Regimen", 1:n_reg)
+      # The default label is the regimen id, so a result is labeled with the
+      # regimen it belongs to. A regimen whose id could not be recovered keeps
+      # the generic text, since a missing label would break plotting.
+      sim_labels <- reg_ids
+      missing_ids <- is.na(sim_labels)
+      sim_labels[missing_ids] <- paste("Regimen", 1:n_reg)[missing_ids]
 
       if (length(simlabels) > 0) { # replace generic labels with user labels
         n_reglabels <- length(simlabels)
@@ -912,7 +918,8 @@ makePTAtarget <- function(x) {
 #' @param grid `r template("grid")`
 #' @param legend `r template("legend")` Default will be the labeled regimen names as an argument
 #' when creating a [PM_pta] object,
-#' or if missing, "Regimen 1, Regimen 2,...Regimen n", where *n* is the number of
+#' or if missing, the `id` of each regimen, or "Regimen 1, Regimen 2,...Regimen n"
+#' for a regimen whose id is not recoverable, where *n* is the number of
 #' regimens in the PM_pta object.
 #' @param ci Confidence interval around curves on `type = "pdi"` plot, on scale of 0 to 1. Default is 0.9.
 #' @param xlab `r template("xlab")`  Default is "Target" when targets are discrete,
diff --git a/man/PM_pta.Rd b/man/PM_pta.Rd
index 16b27341d..a462de6a4 100755
--- a/man/PM_pta.Rd
+++ b/man/PM_pta.Rd
@@ -23,9 +23,10 @@ appropriate regimens.
 
 Every result is labeled with the \code{id} of the regimen it was calculated from,
 so results can be associated with the requested regimens without depending on
-the order in which the regimens were simulated. \code{reg_num} and the \code{simlabels}
-annotation are positional: they follow the order of the regimens in
-\code{simdata}, which for a \link{PM_sim} is the order of the ids in the data.
+the order in which the regimens were simulated. By default that id is also the
+\code{label}. \code{reg_num} is positional: it is the order in which the id first appears
+in \code{simdata}, which for a \link{PM_sim} is the order of the ids in the data, and
+explicit \code{simlabels} are assigned in that order.
 }
 \examples{
 
@@ -134,7 +135,7 @@ Michael Neely and Jan Strojil
 \item \code{NPex$data$standard_data} \emph{PM_data_data}
 \item \code{"simout.txt"} matches files with wildcard ability, see \link{PM_sim}
 }}
-      \item{\code{simlabels}}{Optional character vector of labels for each simulation.  Default is \code{c('Regimen 1', 'Regimen 2',...)}.}
+      \item{\code{simlabels}}{Optional character vector of labels for each simulation.  Default is the \code{id} of each regimen, falling back to \code{c('Regimen 1', 'Regimen 2',...)} for a regimen whose id is not recoverable. Labels are assigned in the order of \code{reg_num}, so the first label is for the first regimen.}
       \item{\code{target}}{One of several options.
 \itemize{
 \item A vector of pharmacodynamic targets, such as Minimum Inhibitory Concentrations (MICs), e.g. \code{c(0.25, 0.5, 1, 2, 4, 8, 16, 32)}.
@@ -201,7 +202,7 @@ of \code{target}s and simulated regimens for a given \code{target_type}. The tib
 \itemize{
 \item \strong{reg_num} The position of the simulation in \code{simdata}, i.e. the order in which its id first appears.
 \item \strong{id} The identifier of the simulated regimen, as it appears in the \code{id} column of the simulation.
-\item \strong{label} Annotation of the simulation, supplied by the \code{simlabels} argument.
+\item \strong{label} Annotation of the simulation, supplied by the \code{simlabels} argument, or the regimen \code{id} by default.
 \item \strong{target} is the specified \code{target} for the results row. If a distribution created by \link{makePTAtarget},
 this will be a tibble with  the simulated targets
 \item \strong{type} is the specified \code{target_type} for the results row
diff --git a/man/plot.PM_pta.Rd b/man/plot.PM_pta.Rd
index 03ce271b2..5252fbaae 100755
--- a/man/plot.PM_pta.Rd
+++ b/man/plot.PM_pta.Rd
@@ -114,7 +114,8 @@ The default value is 12.
 Examples: \itemize{ \item{\code{legend = T}} \item{\code{legend = list(orientation = "h", font = list(color = "blue"))}} }
 Default will be the labeled regimen names as an argument
 when creating a \link{PM_pta} object,
-or if missing, "Regimen 1, Regimen 2,...Regimen n", where \emph{n} is the number of
+or if missing, the \code{id} of each regimen, or "Regimen 1, Regimen 2,...Regimen n"
+for a regimen whose id is not recoverable, where \emph{n} is the number of
 regimens in the PM_pta object.}
 
 \item{log}{Boolean operator to plot the x axis in log base 10. This argument maps to
diff --git a/tests/testthat/test-pm-pta-regimen-order.R b/tests/testthat/test-pm-pta-regimen-order.R
index f8392ea90..91d5c1012 100644
--- a/tests/testthat/test-pm-pta-regimen-order.R
+++ b/tests/testthat/test-pm-pta-regimen-order.R
@@ -71,10 +71,10 @@ pta_fixture <- function(n_ids = 12, nsim = 4,
   )
 }
 
-pta_pdi <- function(sim, mic = 2) {
+pta_pdi <- function(sim, mic = 2, ...) {
   suppressMessages(PM_pta$new(
     simdata = sim, target = mic, target_type = "time", success = 0.5,
-    outeq = 1, free_fraction = 1, start = 0, end = Inf
+    outeq = 1, free_fraction = 1, start = 0, end = Inf, ...
   ))$data$data
 }
 
@@ -196,17 +196,79 @@ testthat::test_that("a user-supplied list of regimens is associated by its id",
   }
 })
 
+testthat::test_that("a PTA without simlabels is labeled with the regimen id", {
+  for (shape in c("numeric", "character", "prefixed")) {
+    fx <- pta_fixture(shape = shape)
+    rows <- pta_pdi(fx$sim_all)
+
+    testthat::expect_equal(rows$label, as.character(rows$id), info = shape)
+    testthat::expect_equal(rows$label, as.character(fx$data_ids), info = shape)
+    # Guard against a vacuous pass: the previous default was the generic text.
+    testthat::expect_false(any(rows$label == paste("Regimen", rows$reg_num)), info = shape)
+
+    summary <- suppressMessages(PM_pta$new(
+      simdata = fx$sim_all, target = 2, target_type = "time", success = 0.5,
+      outeq = 1, free_fraction = 1, start = 0, end = Inf
+    )$summary())
+    testthat::expect_equal(summary$label, summary$id, info = shape)
+  }
+})
+
+testthat::test_that("id-based default labels plot", {
+  # 4 regimens, since the default color palette holds at most 9.
+  fx <- pta_fixture(shape = "prefixed", n_ids = 4)
+  pta <- suppressMessages(PM_pta$new(
+    simdata = fx$sim_all, target = c(1, 2), target_type = "time", success = 0.5,
+    outeq = 1, free_fraction = 1, start = 0, end = Inf
+  ))
+
+  testthat::expect_no_error(suppressMessages(pta$plot(print = FALSE)))
+})
+
+testthat::test_that("supplied simlabels stay positional, in reg_num order", {
+  fx <- pta_fixture(shape = "prefixed")
+  labels <- paste0("L", seq_along(fx$ids))
+  rows <- pta_pdi(fx$sim_all, simlabels = labels)
+
+  testthat::expect_equal(rows$label, labels[rows$reg_num])
+  # The label is the supplied one, while the id still identifies the regimen.
+  testthat::expect_equal(rows$id, as.character(fx$data_ids)[rows$reg_num])
+})
+
+testthat::test_that("a regimen whose id is not recoverable keeps a generic label", {
+  # 4 regimens, since the default color palette holds at most 9.
+  fx <- pta_fixture(shape = "numeric", n_ids = 4)
+  obs <- fx$sim_all$data$obs
+  # Drop the id from each regimen, as an unnamed list of objects without ids.
+  regimens <- unname(lapply(split(obs, as.factor(obs$id)), function(x) {
+    x[, setdiff(names(x), "id"), drop = FALSE]
+  }))
+  # "min" is used rather than "time" because the latter arranges by id.
+  pta <- suppressMessages(PM_pta$new(
+    simdata = regimens, target = 2, target_type = "min", success = 1,
+    outeq = 1, free_fraction = 1, start = 0, end = Inf
+  ))
+  rows <- pta$data$data
+
+  testthat::expect_equal(rows$label, paste("Regimen", seq_along(regimens)))
+  testthat::expect_true(all(is.na(rows$id)))
+  # A missing label would break the legend, which is built from the labels.
+  testthat::expect_no_error(suppressMessages(pta$plot(print = FALSE)))
+})
+
 testthat::test_that("results saved before the id column still summarise and plot", {
   pta <- suppressMessages(PM_pta$new(
     simdata = simEx, target = list(5, 10), target_type = c("min", "max"),
     success = c(1, 1), start = 120, end = 144
   ))
   # Emulate a PTA made before results carried the regimen id. Its identifiers
-  # were replaced by a rank, so the regimen number is the only one left.
+  # were replaced by a rank and its labels were the generic text, so the regimen
+  # number is the only identifier left.
   legacy <- pta$clone(deep = TRUE)
   for (i in seq_along(legacy$data)) {
     x <- legacy$data[[i]]
-    if (is.data.frame(x) && "id" %in% names(x)) {
+    if (is.data.frame(x)) {
+      x$label <- paste("Regimen", x$reg_num)
       legacy$data[[i]] <- x[, setdiff(names(x), "id"), drop = FALSE]
     }
   }
@@ -214,6 +276,7 @@ testthat::test_that("results saved before the id column still summarise and plot
 
   summary <- suppressMessages(legacy$summary())
   testthat::expect_equal(summary$id, as.character(summary$reg_num))
+  testthat::expect_equal(summary$label, paste("Regimen", summary$reg_num))
   single <- suppressMessages(legacy$summary(at = 1))
   testthat::expect_equal(single$id, as.character(single$reg_num))
   testthat::expect_no_error(suppressMessages(legacy$plot(print = FALSE)))

From 630b61d9f5c14604da96da4509aaa43127cc16d7 Mon Sep 17 00:00:00 2001
From: Markus <66058642+mhovd@users.noreply.github.com>
Date: Wed, 16 Sep 2026 14:19:22 +0200
Subject: [PATCH 4/6] refactor: build PTA results from a regimen x target table

`private$make()` identified a regimen by index arithmetic: it nested a map over
regimens and target types, then looked up `reg_ids[x]`, `simdata[[x]]`,
`target[[y]]`, `success[[y]]` and `start[y]`, and rebuilt the name of the `pta_*`
helper as a string in order to call it. Nothing about a result carried the
regimen it came from.

It now builds one row per regimen and target, with the regimen's data, id, start
and end in the row, and maps over those rows with `purrr::pmap`, choosing the
helper from a lookup keyed by the target type. `expand_grid()` is used rather
than `crossing()` because it preserves the order of its inputs, so the rows stay
regimen-major and then follow the order of the target types and values, which is
the order results have always been reported in. The `-` sign handling, the
progress bar, the sampled-target branch, the `this_type` split and the
`intersect` assembly are unchanged.

This exposed a bug in `start` and `end`: they were recycled to one value per
regimen, then remapped with `start[x]` where `x` is the *target type* index, so
the per-regimen values were discarded and every regimen silently got the first
regimen's window. A vector `start` or `end` is now honoured per regimen, and a
numeric `target_type` still overrides both with its own absolute value.

`pm_id_rank()`'s fallback sorted with `sort()`, which follows the session
collation, while dplyr's `arrange()` and `dense_rank()` use the C locale, so the
order of the regimens and of `standard_data` depended on the user's locale. It
now sorts with `method = "radix"` (byte order) and keeps `na.last = TRUE`.

Simulating from posteriors paired `postToUse[i]` with `toInclude[i]` by
position, so a posterior file listing its ids in a different order than the data
attached the wrong posterior to each subject. The two are paired by id now, and
a template with no posterior aborts with a clear message. `NPex$postMean` is
empty, so this cannot be exercised from the shipped data.

Verified with a parity harness against the previous commit: 26 of 31 fixtures
produce identical `PM_pta$data`, and the 5 that differ are the per-regimen
`start`/`end` cases, differing only in the intended way.

Full testthat suite: no failures.
---
 R/PM_pta.R                                 | 141 +++++++++++----------
 R/PM_sim.R                                 |  12 ++
 R/id_order.R                               |   4 +-
 man/PM_pta.Rd                              |   8 +-
 tests/testthat/test-pm-pta-regimen-order.R |  75 ++++++++++-
 5 files changed, 169 insertions(+), 71 deletions(-)

diff --git a/R/PM_pta.R b/R/PM_pta.R
index bb97762d4..8f2401b0c 100755
--- a/R/PM_pta.R
+++ b/R/PM_pta.R
@@ -92,11 +92,13 @@ PM_pta <- R6::R6Class(
     #' @param free_fraction Proportion of free, active drug, expressed as a numeric value >=0 and <=1.  Default is 1, i.e.,
     #' 100% free drug or 0% protein binding.
     #' @param start Specify the time to begin PTA calculations. Default is a vector with the first observation time for subjects
-    #' in each element of `simdata`, e.g. dose regimen. If specified as a vector, values will be recycled as necessary.
+    #' in each element of `simdata`, e.g. dose regimen. If specified as a vector, one value per regimen is used, recycled as
+    #' necessary, so each regimen keeps its own start. A numerical `target_type` overrides `start` for that target with its own absolute value.
     #' @param end Specify the time to end PTA calculations so that PTA is calculated
     #' from `start` to `end`.  Default for end is the maximum observation
-    #' time for subjects in each element of `simdata`, e.g. dose regimen.  If specified as a vector, values will be recycled
-    #' as necessary. Subjects with insufficient data (fewer than 5 simulated observations) for a specified interval will trigger a warning.
+    #' time for subjects in each element of `simdata`, e.g. dose regimen.  If specified as a vector, one value per regimen is used,
+    #' recycled as necessary, so each regimen keeps its own end, and a numerical `target_type` overrides it as for `start`.
+    #' Subjects with insufficient data (fewer than 5 simulated observations) for a specified interval will trigger a warning.
     #' Ideally then, the simulated datset should contain sufficient observations within the interval specified by `start` and `end`.
     #' @param icen Can be either "median" for the predictions based on medians of `pred.type` parameter value
     #' distributions, or "mean".  Default is "median".
@@ -363,13 +365,9 @@ PM_pta <- R6::R6Class(
         simTarg <- FALSE
       }
 
-      # fill in start and end times for each regimen
-      if (length(start) < n_reg) {
-        start <- rep(start, n_reg)[1:n_reg]
-      }
-      if (length(end) < n_reg) {
-        end <- rep(end, n_reg)[1:n_reg]
-      }
+      # start and end are per regimen, and a vector is recycled as necessary.
+      reg_start <- rep_len(start, n_reg)
+      reg_end <- rep_len(end, n_reg)
 
       # check for valid arguments
       target_type <- stringr::str_replace_all(target_type, "peak", "max")
@@ -386,21 +384,9 @@ PM_pta <- R6::R6Class(
         ))
       }
 
-      # adjust start and end for any specific times
-      start <- unlist(map(1:n_type, \(x) {
-        if (suppressWarnings(!is.na(as.numeric(target_type[x])))) {
-          abs(as.numeric(target_type[x]))
-        } else {
-          start[x]
-        }
-      }))
-      end <- unlist(map(1:n_type, \(x) {
-        if (suppressWarnings(!is.na(as.numeric(target_type[x])))) {
-          abs(as.numeric(target_type[x]))
-        } else {
-          end[x]
-        }
-      }))
+      # A numeric `target_type` is a specific time, and both start and end become
+      # that time. It is applied to each (regimen, target) row below, so the
+      # regimen's own start and end are kept for every other target type.
 
 
       # check to make sure secondary targets are only length 1
@@ -488,11 +474,18 @@ PM_pta <- R6::R6Class(
       on.exit(pb$terminate(), add = TRUE)
 
       ###### MAKE THE PTA OBJECT
-      master_pta <- purrr::map(1:n_reg, \(x){
-        purrr::map(1:n_type, \(y) tidyr::expand_grid(
-          reg_num = x,
-          id = reg_ids[x],
-          target = if (simTarg) { # simulated targets
+      # One row per regimen and target, carrying the regimen's own data, id,
+      # start and end rather than an index into them. `expand_grid()` keeps the
+      # order of its inputs - unlike `crossing()`, which sorts - so the rows stay
+      # regimen-major and then follow the order of the target types and values,
+      # which is the order results have always been reported in.
+      regimens <- tibble::tibble(reg_num = seq_along(simdata), id = reg_ids, sims = simdata)
+
+      # A regimen's targets are the same for every regimen, except that sampled
+      # targets draw one target per simulated subject of that regimen.
+      targets_for <- function(x) {
+        purrr::map(seq_len(n_type), \(y) {
+          this_target <- if (simTarg) { # simulated targets
             if (y == 1) {
               list(
                 tidyr::tibble(
@@ -505,50 +498,70 @@ PM_pta <- R6::R6Class(
                 )
               )
             } else {
-              target[y]
+              list(target[[y]])
             }
           } else {
-            target[[y]]
-          }, # discrete targets
-          this_type = y,
-          type = target_type[[y]],
-          success_ratio = success[[y]],
-          start = start[y],
-          end = end[y]
-          # sim_data = list(simdata[[x]])
-        )) |> # end inner map
-          purrr::list_rbind()
-      }) |> # end outer map
-        purrr::list_rbind() |>
-        mutate(type = stringr::str_replace_all(type, "\\d+", "specific")) |>
-        rowwise() |>
-        mutate(pdi = list(do.call(
-          paste0("pta_", stringr::str_replace_all(type, "-", "")), # call the appropriate pta function below
-          list( # arguments to the pta function
-            sims = simdata[[reg_num]],
-            # sims = sim_data,
-            .target = target,
-            .simTarg = simTarg,
-            .start = start,
-            .end = end,
-            .pb = pb
+            target[[y]] # discrete targets keep their own type
+          }
+          # a numeric target type is a specific time, and both start and end
+          # become that time; otherwise the regimen's own values are used
+          specific <- suppressWarnings(as.numeric(target_type[y]))
+          tibble::tibble(
+            target = this_target,
+            type = target_type[[y]],
+            success_ratio = success[[y]],
+            this_type = y,
+            start = if (is.na(specific)) reg_start[x] else abs(specific),
+            end = if (is.na(specific)) reg_end[x] else abs(specific)
           )
-        ))) |>
-        mutate(success = list(purrr::map_dbl(pdi, \(x) {
+        }) |> purrr::list_rbind()
+      }
+
+      grid <- purrr::map(seq_len(n_reg), \(x) {
+        tidyr::expand_grid(regimens[x, ], targets_for(x))
+      }) |>
+        purrr::list_rbind() |>
+        mutate(type = stringr::str_replace_all(type, "\\d+", "specific"))
+
+      # choose the pta function by name rather than rebuilding its name
+      pta_helpers <- list(
+        time = pta_time, auc = pta_auc, max = pta_max,
+        min = pta_min, specific = pta_specific
+      )
+
+      # the helper lookup is a named list, so drop the names it would pass on
+      grid$pdi <- unname(purrr::pmap(
+        list(
+          pta_helpers[stringr::str_replace_all(grid$type, "-", "")],
+          grid$sims, grid$target, grid$start, grid$end
+        ),
+        \(pta_fn, sims, .target, .start, .end) {
+          pta_fn(sims, .target, simTarg, .start, .end, pb)
+        }
+      ))
+
+      grid$success <- unname(purrr::pmap(
+        list(grid$pdi, grid$type, grid$success_ratio),
+        \(pdi, type, success_ratio) {
           if (stringr::str_detect(type, "-")) {
-            x <= success_ratio # will return NA is x is NA
+            as.numeric(pdi <= success_ratio) # will return NA if pdi is NA
           } else {
-            x >= success_ratio
+            as.numeric(pdi >= success_ratio)
           }
-        }))) |>
-        mutate(prop_success = sum(success) / length(success)) |>
-        mutate(label = sim_labels[reg_num]) |>
+        }
+      ))
+
+      master_pta <- grid |>
+        dplyr::select(-sims) |>
+        mutate(
+          prop_success = purrr::map_dbl(success, \(x) sum(x) / length(x)),
+          label = sim_labels[reg_num]
+        ) |>
         dplyr::relocate(
           reg_num, id, label, target, type,
           success_ratio, prop_success, success, pdi,
           start, end
-        ) |>
-        ungroup() # remove rowwise
+        )
 
 
       # add intersection if multiple target types
diff --git a/R/PM_sim.R b/R/PM_sim.R
index 5a4e42a97..e5793594e 100755
--- a/R/PM_sim.R
+++ b/R/PM_sim.R
@@ -1029,6 +1029,18 @@ PM_sim <- R6::R6Class(
             if (length(postToUse) > 0 && length(postToUse) != nsub) {
               cli::cli_abort(c("x" = "You have {length(postToUse)} posteriors and {nsub} selected subjects in the data file.  These must be equal."))
             }
+            # Pair each template with its own posterior by id rather than by
+            # position, so a posterior file that lists its ids in a different
+            # order than the data still attaches the right posterior to each
+            # subject. `NPex$postMean` is empty, so this cannot be exercised from
+            # the shipped data.
+            postToUse <- postToUse[match(toInclude, postToUse)]
+            if (anyNA(postToUse)) {
+              cli::cli_abort(c(
+                "x" = "No posterior mean was found for subject{?s} {.val {toInclude[is.na(postToUse)]}}.",
+                "i" = "Simulating from posteriors needs one posterior for every included subject."
+              ))
+            }
           }
         } else {
           postToUse <- NULL
diff --git a/R/id_order.R b/R/id_order.R
index f0f87c630..6956d76fe 100644
--- a/R/id_order.R
+++ b/R/id_order.R
@@ -21,7 +21,9 @@ pm_id_rank <- function(id) {
   ordered <- if (all(!is.na(numeric_ids))) {
     unique_ids[order(numeric_ids)]
   } else {
-    sort(unique_ids, na.last = TRUE)
+    # method = "radix" sorts in the C locale, the byte order dplyr's arrange()
+    # and dense_rank() use, so the order does not follow the session collation.
+    sort(unique_ids, na.last = TRUE, method = "radix")
   }
 
   match(id, ordered)
diff --git a/man/PM_pta.Rd b/man/PM_pta.Rd
index a462de6a4..0907c5483 100755
--- a/man/PM_pta.Rd
+++ b/man/PM_pta.Rd
@@ -172,11 +172,13 @@ For example \code{success = 0.4} or \code{success = c(0.4, 100)}. The length mus
       \item{\code{free_fraction}}{Proportion of free, active drug, expressed as a numeric value >=0 and <=1.  Default is 1, i.e.,
 100\% free drug or 0\% protein binding.}
       \item{\code{start}}{Specify the time to begin PTA calculations. Default is a vector with the first observation time for subjects
-in each element of \code{simdata}, e.g. dose regimen. If specified as a vector, values will be recycled as necessary.}
+in each element of \code{simdata}, e.g. dose regimen. If specified as a vector, one value per regimen is used, recycled as
+necessary, so each regimen keeps its own start. A numerical \code{target_type} overrides \code{start} for that target with its own absolute value.}
       \item{\code{end}}{Specify the time to end PTA calculations so that PTA is calculated
 from \code{start} to \code{end}.  Default for end is the maximum observation
-time for subjects in each element of \code{simdata}, e.g. dose regimen.  If specified as a vector, values will be recycled
-as necessary. Subjects with insufficient data (fewer than 5 simulated observations) for a specified interval will trigger a warning.
+time for subjects in each element of \code{simdata}, e.g. dose regimen.  If specified as a vector, one value per regimen is used,
+recycled as necessary, so each regimen keeps its own end, and a numerical \code{target_type} overrides it as for \code{start}.
+Subjects with insufficient data (fewer than 5 simulated observations) for a specified interval will trigger a warning.
 Ideally then, the simulated datset should contain sufficient observations within the interval specified by \code{start} and \code{end}.}
       \item{\code{icen}}{Can be either "median" for the predictions based on medians of \code{pred.type} parameter value
 distributions, or "mean".  Default is "median".}
diff --git a/tests/testthat/test-pm-pta-regimen-order.R b/tests/testthat/test-pm-pta-regimen-order.R
index 91d5c1012..b50f407f0 100644
--- a/tests/testthat/test-pm-pta-regimen-order.R
+++ b/tests/testthat/test-pm-pta-regimen-order.R
@@ -71,10 +71,11 @@ pta_fixture <- function(n_ids = 12, nsim = 4,
   )
 }
 
-pta_pdi <- function(sim, mic = 2, ...) {
+pta_pdi <- function(sim, mic = 2, target_type = "time", success = 0.5,
+                    start = 0, end = Inf, ...) {
   suppressMessages(PM_pta$new(
-    simdata = sim, target = mic, target_type = "time", success = 0.5,
-    outeq = 1, free_fraction = 1, start = 0, end = Inf, ...
+    simdata = sim, target = mic, target_type = target_type, success = success,
+    outeq = 1, free_fraction = 1, start = start, end = end, ...
   ))$data$data
 }
 
@@ -256,6 +257,74 @@ testthat::test_that("a regimen whose id is not recoverable keeps a generic label
   testthat::expect_no_error(suppressMessages(pta$plot(print = FALSE)))
 })
 
+testthat::test_that("start and end are honoured per regimen", {
+  fx <- pta_fixture(shape = "numeric", n_ids = 4)
+  starts <- c(0, 24, 48, 72)
+  ends <- c(96, 120, 144, 168)
+
+  # "min" is the minimum concentration within the window, so it depends on the
+  # window: a regimen given the wrong start or end gets a different PDI.
+  rows <- pta_pdi(fx$sim_all, mic = 2, target_type = "min", success = 1,
+                  start = starts, end = ends)
+
+  testthat::expect_equal(rows$start, starts[rows$reg_num])
+  testthat::expect_equal(rows$end, ends[rows$reg_num])
+
+  # Guard against a vacuous pass: the windows must give distinguishable results.
+  testthat::expect_gte(length(unique(vapply(rows$pdi, pdi_fingerprint, character(1)))), 3)
+
+  for (i in seq_len(nrow(rows))) {
+    one <- pta_pdi(fx$sim_for(rows$id[i]), mic = 2, target_type = "min", success = 1,
+                   start = starts[rows$reg_num[i]], end = ends[rows$reg_num[i]])
+    testthat::expect_equal(
+      round(unlist(rows$pdi[[i]]), 8),
+      round(unlist(one$pdi[[1]]), 8),
+      info = sprintf("row %d (id %s) should use start %s, end %s",
+                     i, rows$id[i], starts[rows$reg_num[i]], ends[rows$reg_num[i]])
+    )
+  }
+})
+
+testthat::test_that("a numeric target_type overrides start and end for every regimen", {
+  fx <- pta_fixture(shape = "numeric", n_ids = 4)
+  rows <- pta_pdi(fx$sim_all, mic = 2, target_type = 144, success = 0.5,
+                  start = c(0, 24, 48, 72))
+
+  # A numeric target_type is a specific time, and start and end both become it.
+  testthat::expect_equal(unique(rows$start), 144)
+  testthat::expect_equal(unique(rows$end), 144)
+  testthat::expect_equal(unique(rows$type), "specific")
+})
+
+testthat::test_that("id ranking uses byte order, not the session collation", {
+  pm_id_rank <- getFromNamespace("pm_id_rank", "Pmetrics")
+
+  # Byte order of these is A2 < B < a < a1, because uppercase letters come
+  # first. A collation that orders lowercase before uppercase gives a different
+  # answer, so the ranks are pinned to byte order rather than to whatever the
+  # session happens to use.
+  ids <- c("a", "B", "a1", "A2")
+  byte_order_ranks <- c(3, 2, 4, 1)
+  testthat::expect_equal(pm_id_rank(ids), byte_order_ranks)
+  testthat::expect_equal(order(pm_id_rank(ids)), order(ids, method = "radix"))
+
+  # The ranks must not change when the collation does. testthat pins
+  # LC_COLLATE to C, which is byte order, so a collation is set explicitly here
+  # to make the difference observable; skip if the locale is not installed, as
+  # on a minimal Linux runner.
+  original <- Sys.getlocale("LC_COLLATE")
+  if (!identical(Sys.setlocale("LC_COLLATE", "en_US.UTF-8"), "")) {
+    ranks <- pm_id_rank(ids)
+    Sys.setlocale("LC_COLLATE", original)
+    testthat::expect_equal(ranks, byte_order_ranks)
+  }
+
+  # ids that are numbers keep numeric order, and missing ids sort last
+  testthat::expect_equal(pm_id_rank(c("10", "2", "1")), c(3, 2, 1))
+  testthat::expect_equal(pm_id_rank(c("b", NA, "a")), c(2, 3, 1))
+  testthat::expect_equal(pm_id_rank(character(0)), integer(0))
+})
+
 testthat::test_that("results saved before the id column still summarise and plot", {
   pta <- suppressMessages(PM_pta$new(
     simdata = simEx, target = list(5, 10), target_type = c("min", "max"),

From c83f5be4e00dcae1189df31a599dc67a9136c38c Mon Sep 17 00:00:00 2001
From: Markus <66058642+mhovd@users.noreply.github.com>
Date: Wed, 16 Sep 2026 14:40:50 +0200
Subject: [PATCH 5/6] fix: pair each subject with its posterior by position,
 not by id

`postToUse` was matched to the templates by id, but `getSimPrior()` consumes that
value as a row index into `postMean` and as a list index into `postCov`, so the
pairing was still positional. With posterior ids `c(10, 20)` and templates
`c(20, 10)` the indices became `c(20, 10)` for a two-row frame: the wrong
posterior, or out of bounds. Ids and positions coincide only when the posteriors
happen to be numbered 1..n in the order they are stored, which is why the
existing tests passed. The conflation is pre-existing; matching by id made it
reachable for ids that are not 1..n.

`match_posteriors()` returns the position of each template's posterior and
aborts on duplicate posterior ids, which is what makes a single `match()` valid
for both structures, or on a template with no posterior, naming the subjects.
The parameter is carried as `post_pos` through `getSimPrior()` and both call
sites. The count check it replaces reported the same situation less clearly.

`postCov` is built by splitting the posterior on id, so it carries the ids as
names, while `postMean` is ordered by a dplyr group; those orderings agree for
ids that are numbers but need not for others. `posterior_cov()` uses the name
when it is available and falls back to the position.

Tested against the shipped `NPex$final`, which has a 20-row `postMean` and 20
`postCov` entries. With `nsim = 1` the posterior mean is drawn exactly, so the
simulated parameters reveal which posterior each subject was given. Reversing
the stored order, relabelling the ids to `10k`, and the helper's own shape
(`c(20, 10)` against `c(10, 20)`) all pass here and all fail before this change:
3 failures and 2 errors.

Full testthat suite: no failures.
---
 R/PM_sim.R                          | 85 ++++++++++++++++++---------
 tests/testthat/test-pm-sim-inputs.R | 91 +++++++++++++++++++++++++++++
 2 files changed, 149 insertions(+), 27 deletions(-)

diff --git a/R/PM_sim.R b/R/PM_sim.R
index e5793594e..4818b0d32 100755
--- a/R/PM_sim.R
+++ b/R/PM_sim.R
@@ -1019,31 +1019,20 @@ PM_sim <- R6::R6Class(
               "i" = "Please remake your model."
             ))
           } else {
-            postToUse <- unique(poppar$postMean$id) # get the id for each posterior mean
             if (split) {
               split <- FALSE
               cli::cli_inform(c("i" = "{.arg split} set to {.val FALSE} for simulations from posteriors."))
             }
-            # now, filter by included IDs
-            postToUse <- postToUse[postToUse %in% toInclude] # subset the posteriors if applicable
-            if (length(postToUse) > 0 && length(postToUse) != nsub) {
-              cli::cli_abort(c("x" = "You have {length(postToUse)} posteriors and {nsub} selected subjects in the data file.  These must be equal."))
-            }
             # Pair each template with its own posterior by id rather than by
-            # position, so a posterior file that lists its ids in a different
-            # order than the data still attaches the right posterior to each
-            # subject. `NPex$postMean` is empty, so this cannot be exercised from
-            # the shipped data.
-            postToUse <- postToUse[match(toInclude, postToUse)]
-            if (anyNA(postToUse)) {
-              cli::cli_abort(c(
-                "x" = "No posterior mean was found for subject{?s} {.val {toInclude[is.na(postToUse)]}}.",
-                "i" = "Simulating from posteriors needs one posterior for every included subject."
-              ))
-            }
+            # position, so posteriors listed in a different order than the data
+            # still attach to the right subject. Positions are kept, because
+            # `getSimPrior()` indexes the posterior mean and covariance with
+            # them; ids and positions only coincide when the posteriors happen to
+            # be numbered 1..n in the order they are stored.
+            post_pos <- match_posteriors(toInclude, poppar$postMean$id)
           }
         } else {
-          postToUse <- NULL
+          post_pos <- NULL
         }
         
         
@@ -1414,7 +1403,7 @@ PM_sim <- R6::R6Class(
           mod <- model
         }
         
-        if (length(postToUse) > 0) {
+        if (length(post_pos) > 0) {
           # simulating from posteriors, each posterior matched to a subject
           # need to set theta as the posterior mean or median for each subject
           ans <- NULL
@@ -1425,7 +1414,7 @@ PM_sim <- R6::R6Class(
               i = i,
               poppar = poppar,
               split = split,
-              postToUse = postToUse[i],
+              post_pos = post_pos[i],
               limits = limits,
               seed = seed[1],
               nsim = nsim,
@@ -1471,7 +1460,7 @@ PM_sim <- R6::R6Class(
           
           class(ret) <- c("PM_sim_data", "list")
           self$data <- ret
-        } else { # postToUse is false
+        } else { # not simulating from posteriors
           
           # set theta as nsim rows drawn from prior
           if (!useTheta) {
@@ -1479,7 +1468,7 @@ PM_sim <- R6::R6Class(
               i = 1,
               poppar = poppar,
               split = split,
-              postToUse = NULL,
+              post_pos = NULL,
               limits = limits,
               seed = seed[1],
               nsim = nsim,
@@ -1549,7 +1538,7 @@ PM_sim <- R6::R6Class(
       }, # end of SIMrun
       
       # get prior density
-      getSimPrior = function(i, poppar, split, postToUse, limits, seed, nsim, toInclude, msg = NULL) {
+      getSimPrior = function(i, poppar, split, post_pos, limits, seed, nsim, toInclude, msg = NULL) {
         # get prior density
         
         
@@ -1580,15 +1569,15 @@ PM_sim <- R6::R6Class(
             ndist <- nrow(popPoints)
             pop_cov <- poppar$popCov / ndist
           } else { # not split
-            if (is.null(postToUse)) { # not simulating from posteriors
+            if (is.null(post_pos)) { # not simulating from posteriors
               pop_weight <- 1
               pop_mean <- poppar$popMean
               pop_cov <- poppar$popCov
               ndist <- 1
             } else { # simulating from posteriors
               pop_weight <- 1
-              pop_mean <- poppar$postMean[postToUse, ] |> select(-id)
-              pop_cov <- poppar$postCov[[postToUse]]
+              pop_mean <- poppar$postMean[post_pos, ] |> select(-id)
+              pop_cov <- posterior_cov(poppar$postCov, post_pos, toInclude[i])
               ndist <- 1
             }
           }
@@ -1611,7 +1600,7 @@ PM_sim <- R6::R6Class(
         if (length(pop_cov) == 1 && pop_cov == 1) {
           return(invisible(NULL)) # quietly abort simulation
         } else if (length(pop_cov) == 1 && pop_cov == -1) {
-          msg <- if (!is.null(postToUse)) {
+          msg <- if (!is.null(post_pos)) {
             glue::glue("Unable to fix covariance for template {.code id = {toInclude[i]}}.")
           } else {
             "Unable to make population covariance positive definite."
@@ -2726,3 +2715,45 @@ generate_multimodal_samples <- function(num_samples, weights, means, cov_matrix,
     total_nsim = total_nsim
   ))
 }
+
+# POSTERIOR PAIRING -------------------------------------------------------
+
+# Pair each template with its own posterior.
+#
+# Returns the position of each template's posterior within `post_ids`, which is
+# the row of `postMean` and the entry of `postCov` belonging to that subject.
+# Positions are returned rather than ids because `getSimPrior()` indexes both
+# structures with the result, and ids coincide with positions only when the
+# posteriors happen to be numbered 1..n in the order they are stored.
+match_posteriors <- function(toInclude, post_ids) {
+  if (anyDuplicated(post_ids)) {
+    cli::cli_abort(c(
+      "x" = "Posterior means contain duplicate subject ids.",
+      "i" = "Each subject needs exactly one posterior."
+    ))
+  }
+
+  post_pos <- match(toInclude, post_ids)
+  if (anyNA(post_pos)) {
+    cli::cli_abort(c(
+      "x" = "No posterior mean was found for subject{?s} {.val {toInclude[is.na(post_pos)]}}.",
+      "i" = "Simulating from posteriors needs one posterior for every included subject."
+    ))
+  }
+
+  post_pos
+}
+
+# The posterior covariance matrix for one subject.
+#
+# `postCov` is built by splitting the posterior on id, so it carries the ids as
+# names, while `postMean` is ordered by a dplyr group. Those two orderings agree
+# for ids that are numbers but need not for other ids, so the name is used when
+# it is available and the position is the fallback for an object without names.
+posterior_cov <- function(postCov, post_pos, id) {
+  key <- as.character(id)
+  if (!is.null(names(postCov)) && key %in% names(postCov)) {
+    return(postCov[[key]])
+  }
+  postCov[[post_pos]]
+}
diff --git a/tests/testthat/test-pm-sim-inputs.R b/tests/testthat/test-pm-sim-inputs.R
index b67a7b945..1f54dbe03 100644
--- a/tests/testthat/test-pm-sim-inputs.R
+++ b/tests/testthat/test-pm-sim-inputs.R
@@ -220,3 +220,94 @@ test_that("PM_sim handles data.frame poppar and ignores nsim", {
   expect_equal(nrow(sim_points$data$parValues), nrow(pop_points_df))
   expect_equal(max(sim_points$data$parValues$nsim), nrow(pop_points_df))
 })
+
+# Each template must be simulated from its own posterior, matched by id rather
+# than by position. With `nsim = 1` the posterior mean is drawn exactly, so the
+# simulated parameters show which posterior a subject was given.
+
+posterior_params <- function(sim, id, pars) {
+  p <- sim$data$parValues
+  unlist(p[p$id == id, pars])
+}
+
+posterior_mean <- function(poppar, id, pars) {
+  m <- poppar$postMean
+  unlist(m[m$id == id, pars])
+}
+
+sim_from_posteriors <- function(poppar, data, include) {
+  suppressMessages(
+    PM_sim$new(
+      poppar = poppar, data = data, model = NPex$model,
+      include = include, nsim = 1, usePost = TRUE, quiet = TRUE
+    )
+  )
+}
+
+test_that("posteriors are paired by id, not by position", {
+  match_posteriors <- getFromNamespace("match_posteriors", "Pmetrics")
+
+  # posterior ids 10, 20 with the templates ordered 20, 10
+  expect_equal(match_posteriors(c(20, 10), c(10, 20)), c(2, 1))
+  # ids that are not 1..n, in the same order
+  expect_equal(match_posteriors(c(10, 30), 10 * (1:3)), c(1, 3))
+
+  expect_error(match_posteriors(c(1, 2), c(1, 1, 2)), "duplicate")
+  expect_error(match_posteriors(c(1, 5), c(1, 2)), "No posterior mean")
+})
+
+test_that("each subject uses its own posterior when they are stored out of order", {
+  pars <- names(NPex$final$popMean)
+  reversed <- NPex$final$clone(deep = TRUE)
+  n <- nrow(reversed$data$postMean)
+  reversed$data$postMean <- reversed$data$postMean[n:1, ]
+  reversed$data$postCov <- reversed$data$postCov[n:1]
+
+  # the posteriors must differ, or the comparison below would prove nothing
+  expect_false(isTRUE(all.equal(
+    posterior_mean(reversed, 1, pars),
+    posterior_mean(reversed, 3, pars)
+  )))
+
+  sim <- sim_from_posteriors(reversed, NPex$data, include = c(1, 3))
+
+  for (id in c(1, 3)) {
+    expect_equal(
+      posterior_params(sim, id, pars),
+      posterior_mean(reversed, id, pars),
+      info = sprintf("subject %s should use its own posterior", id)
+    )
+  }
+
+  # with no names on postCov the position is the only key available, and the two
+  # structures are aligned because both were reversed
+  unnamed <- reversed$clone(deep = TRUE)
+  names(unnamed$data$postCov) <- NULL
+  sim_unnamed <- sim_from_posteriors(unnamed, NPex$data, include = c(1, 3))
+  expect_equal(
+    posterior_params(sim_unnamed, 3, pars),
+    posterior_mean(unnamed, 3, pars)
+  )
+})
+
+test_that("subjects whose ids are not 1..n are paired by id", {
+  pars <- names(NPex$final$popMean)
+
+  relabelled <- NPex$data$standard_data
+  relabelled$id <- 10 * relabelled$id
+  data_10k <- PM_data$new(as.data.frame(relabelled), quiet = TRUE)
+
+  tenths <- NPex$final$clone(deep = TRUE)
+  tenths$data$postMean$id <- 10 * tenths$data$postMean$id
+  names(tenths$data$postCov) <- 10 * as.numeric(names(tenths$data$postCov))
+
+  sim <- sim_from_posteriors(tenths, data_10k, include = c(10, 30))
+
+  for (id in c(10, 30)) {
+    expect_equal(
+      posterior_params(sim, id, pars),
+      posterior_mean(tenths, id, pars),
+      info = sprintf("subject %s should use its own posterior", id)
+    )
+  }
+})

From 64161a23ea36a508d59248ac6d37b837c512f251 Mon Sep 17 00:00:00 2001
From: Markus <66058642+mhovd@users.noreply.github.com>
Date: Wed, 16 Sep 2026 14:47:37 +0200
Subject: [PATCH 6/6] fix: unwrap a saved PTA object when loading it from a
 file

`$save()` writes the whole `PM_pta` object, so `readRDS()` hands back an R6
object. `backfill_pta_id()` understands a `PM_pta_data` list or a data frame, so
it returned that object untouched, and `private$populate()` then stored the whole
loaded PTA as `self$data`. `PM_pta$new(<file>)` produced an object whose
`$summary()` and `$plot()` failed with "wrong arguments for subsetting an
environment", and the `id` backfill added for legacy objects never ran on that
path, so the claim in the description only held for the `readRDS(f)$summary()`
flow.

The loaded object is now unwrapped to its `$data` before backfilling and
populating, with a path kept for a bare `PM_pta_data` list. `PM_upgrade()` is
applied once to whatever was loaded: it walks the data frames of a bare list and
for an R6 object upgrades the list held in `$data`, so both shapes are covered
and nothing is upgraded twice.

Two tests cover the documented round trip and the legacy case, both loading
through `PM_pta$new(<file>)`: the result list is populated, every table has an
`id` column, `summary()` reports it, and `plot()` works. Both fail before this
change with 2 failures and 2 errors. The legacy test's check that the `id`
column really was stripped now looks at every table rather than at `$data`, which
only exists for a single target type and made that assertion vacuous.

Full testthat suite: no failures.
---
 R/PM_pta.R                                 | 12 ++++-
 tests/testthat/test-pm-pta-regimen-order.R | 56 +++++++++++++++++++++-
 2 files changed, 66 insertions(+), 2 deletions(-)

diff --git a/R/PM_pta.R b/R/PM_pta.R
index 8f2401b0c..a1d7cdb8b 100755
--- a/R/PM_pta.R
+++ b/R/PM_pta.R
@@ -298,7 +298,17 @@ PM_pta <- R6::R6Class(
         )
         private$populate(pta)
       } else { # try simdata as a filename
-        pta <- backfill_pta_id(PM_upgrade(readRDS(simdata)))
+        # `$save()` writes the whole `PM_pta` object, so unwrap the result list
+        # before populating. `PM_upgrade()` handles both shapes: it walks the
+        # data frames of a bare `PM_pta_data` list, and for an R6 object it
+        # upgrades the list held in `$data`. It is applied once, to the object
+        # as loaded, so nothing is upgraded twice.
+        loaded <- PM_upgrade(readRDS(simdata))
+        pta <- if (inherits(loaded, "PM_pta")) {
+          backfill_pta_id(loaded$data)
+        } else {
+          backfill_pta_id(loaded)
+        }
         private$populate(pta)
       }
     },
diff --git a/tests/testthat/test-pm-pta-regimen-order.R b/tests/testthat/test-pm-pta-regimen-order.R
index b50f407f0..7fe325eab 100644
--- a/tests/testthat/test-pm-pta-regimen-order.R
+++ b/tests/testthat/test-pm-pta-regimen-order.R
@@ -86,6 +86,19 @@ single_pdi <- function(fx, id) {
 
 pdi_fingerprint <- function(x) paste(round(unlist(x), 8), collapse = "|")
 
+# The result tibbles inside a `PM_pta_data` list, skipping the `NA` intersect
+# placeholder. `$data` only exists when a single target type was requested, so
+# the elements are used directly rather than going through that name.
+pta_tables <- function(pta_data) Filter(is.data.frame, pta_data)
+
+# A two-target-type PTA, so the `intersect` element is exercised too.
+pta_from <- function(simdata) {
+  suppressMessages(PM_pta$new(
+    simdata = simdata, target = list(5, 10), target_type = c("min", "max"),
+    success = c(1, 1), start = 120, end = 144
+  ))
+}
+
 testthat::test_that("simulated regimens follow numeric id order, not text order", {
   for (shape in c("numeric", "character")) {
     fx <- pta_fixture(shape = shape)
@@ -341,7 +354,7 @@ testthat::test_that("results saved before the id column still summarise and plot
       legacy$data[[i]] <- x[, setdiff(names(x), "id"), drop = FALSE]
     }
   }
-  testthat::expect_false("id" %in% names(legacy$data$data))
+  testthat::expect_false(any(vapply(pta_tables(legacy$data), \(x) "id" %in% names(x), logical(1))))
 
   summary <- suppressMessages(legacy$summary())
   testthat::expect_equal(summary$id, as.character(summary$reg_num))
@@ -350,3 +363,44 @@ testthat::test_that("results saved before the id column still summarise and plot
   testthat::expect_equal(single$id, as.character(single$reg_num))
   testthat::expect_no_error(suppressMessages(legacy$plot(print = FALSE)))
 })
+
+testthat::test_that("a saved PTA round-trips through PM_pta$new()", {
+  pta <- pta_from(simEx)
+  file <- withr::local_tempfile(fileext = ".rds")
+  pta$save(file)
+
+  loaded <- pta_from(file)
+
+  # the result list is populated, not the whole saved object
+  testthat::expect_s3_class(loaded$data, "PM_pta_data")
+  # one table per target type, plus the intersect
+  testthat::expect_length(pta_tables(loaded$data), 3)
+  testthat::expect_true(all(vapply(pta_tables(loaded$data), \(x) "id" %in% names(x), logical(1))))
+  testthat::expect_true("id" %in% names(loaded$data$intersect))
+  testthat::expect_no_error(suppressMessages(loaded$summary()))
+  testthat::expect_no_error(suppressMessages(loaded$plot(print = FALSE)))
+})
+
+testthat::test_that("a saved object without ids is backfilled when loaded from a file", {
+  legacy <- pta_from(simEx)$clone(deep = TRUE)
+  for (i in seq_along(legacy$data)) {
+    x <- legacy$data[[i]]
+    if (is.data.frame(x)) {
+      legacy$data[[i]] <- x[, setdiff(names(x), "id"), drop = FALSE]
+    }
+  }
+  # the saved object really did lack the id column
+  testthat::expect_false(any(vapply(pta_tables(legacy$data), \(x) "id" %in% names(x), logical(1))))
+
+  file <- withr::local_tempfile(fileext = ".rds")
+  legacy$save(file)
+
+  loaded <- pta_from(file)
+
+  testthat::expect_s3_class(loaded$data, "PM_pta_data")
+  testthat::expect_true(all(vapply(pta_tables(loaded$data), \(x) "id" %in% names(x), logical(1))))
+  # the identifiers those rows were keyed on are the regimen numbers
+  summary <- suppressMessages(loaded$summary())
+  testthat::expect_equal(summary$id, as.character(summary$reg_num))
+  testthat::expect_no_error(suppressMessages(loaded$plot(print = FALSE)))
+})