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..a1d7cdb8b 100755
--- a/R/PM_pta.R
+++ b/R/PM_pta.R
@@ -23,6 +23,13 @@
 #' 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. 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
 PM_pta <- R6::R6Class(
@@ -49,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)`.
@@ -85,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".
@@ -102,8 +111,9 @@ 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`.
-    #' * **label** Annotation of the simulation, supplied by the `simlabels` argument.
+    #' * **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, 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
@@ -227,15 +237,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 = dplyr::dense_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 = dplyr::dense_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 +298,17 @@ PM_pta <- R6::R6Class(
         )
         private$populate(pta)
       } else { # try simdata as a filename
-        pta <- 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)
       }
     },
@@ -333,6 +350,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)
@@ -346,13 +375,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")
@@ -369,21 +394,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
@@ -443,7 +456,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)
@@ -466,10 +484,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,
-          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(
@@ -482,50 +508,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, label, target, type,
+          reg_num, id, label, target, type,
           success_ratio, prop_success, success, pdi,
           start, end
-        ) |>
-        ungroup() # remove rowwise
+        )
 
 
       # add intersection if multiple target types
@@ -538,7 +584,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 +600,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 +632,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),
@@ -851,7 +941,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,
@@ -1209,9 +1300,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 +1326,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 +1361,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 +1372,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/R/PM_sim.R b/R/PM_sim.R
index dee1d6a8a..4818b0d32 100755
--- a/R/PM_sim.R
+++ b/R/PM_sim.R
@@ -1019,19 +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 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
         }
         
         
@@ -1402,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
@@ -1413,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,
@@ -1459,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) {
@@ -1467,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,
@@ -1537,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
         
         
@@ -1568,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
             }
           }
@@ -1599,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."
@@ -1639,7 +1640,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 +1855,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 +1879,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
@@ -2713,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/R/id_order.R b/R/id_order.R
new file mode 100644
index 000000000..6956d76fe
--- /dev/null
+++ b/R/id_order.R
@@ -0,0 +1,30 @@
+# 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 {
+    # 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 b2b6a64d3..0907c5483 100755
--- a/man/PM_pta.Rd
+++ b/man/PM_pta.Rd
@@ -20,6 +20,13 @@ 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. 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{
 
@@ -128,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)}.
@@ -165,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".}
@@ -193,8 +202,9 @@ 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{label} Annotation of the simulation, supplied by the \code{simlabels} argument.
+\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, 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/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
new file mode 100644
index 000000000..7fe325eab
--- /dev/null
+++ b/tests/testthat/test-pm-pta-regimen-order.R
@@ -0,0 +1,406 @@
+testthat::skip_on_cran()
+
+# Regression: PTA rows must be associated with the regimen that produced them.
+#
+# 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,
+                        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),
+    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)
+  )
+
+  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.
+  # `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(
+      dose = dose * scales[as.integer(id)],
+      id = ids[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 = 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
+  )
+}
+
+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 = target_type, success = success,
+    outeq = 1, free_fraction = 1, start = start, end = end, ...
+  ))$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 = "|")
+
+# 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)
+    observed <- as.character(unique(fx$sim_all$data$obs$id))
+
+    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("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)
+  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(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("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("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"),
+    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 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)) {
+      x$label <- paste("Regimen", x$reg_num)
+      legacy$data[[i]] <- x[, setdiff(names(x), "id"), drop = FALSE]
+    }
+  }
+  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))
+  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)))
+})
+
+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)))
+})
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)
+    )
+  }
+})