fix: associate PTA results with the requested regimen when ids are text

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

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
         }
 

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

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)
+}

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])
+    )
+  }
+})