From 1cf4dd094a5a8e4bc94f7c5e1a54bcd5b1d9b287 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juli=C3=A1n=20D=2E=20Ot=C3=A1lvaro?= Date: Tue, 15 Sep 2026 17:55:40 +0100 Subject: [PATCH] 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]) + ) + } +})