UUPharmacometrics / xpose

Graphical diagnostics for pharmacometric models
https://uupharmacometrics.github.io/xpose
GNU Lesser General Public License v3.0
55 stars 28 forks source link

Correct sensitivity of irep() to sort order #226

Open jpryby opened 4 months ago

jpryby commented 4 months ago

The current irep() depends upon the input x already being sorted within each repeat, but NONMEM does not have this requirement. The proposed change works whether or not x is sorted this way.

library(xpose)
#> Loading required package: ggplot2
#> Warning: package 'ggplot2' was built under R version 4.2.3
#> 
#> Attaching package: 'xpose'
#> The following object is masked from 'package:stats':
#> 
#>     filter

irep2 <- function(x, quiet = FALSE) {
  if (missing(x)) stop('argument "x" is missing, with no default', call. = FALSE)
  if (is.factor(x)) x <- as.numeric(as.character(x))
  lagcheck <- dplyr::lag(x, default = x[1]) != x
  dupcheck <- duplicated(x)
  check <- dplyr::if_else(lagcheck & dupcheck, 1, 0)
  ilen <- which(check==1)[1] - 1
  if (is.na(ilen)) ilen <- length(x)
  x <- rep(1:(length(x)/ilen), each=ilen)
  msg(c('irep: ', max(x), ' simulations found.'), quiet)
  x
}

reps <- 15

sorted <- rep(c(1,1,1,2,2,2,2,3,3,3),reps)
unsorted <- rep(c(10,10,10,2,2,2,2,3,3,3),reps)

max(irep(sorted))==reps
#> irep: 15 simulations found.
#> [1] TRUE
max(irep(unsorted))==reps
#> irep: 16 simulations found.
#> [1] FALSE

max(irep2(sorted))==reps
#> irep: 15 simulations found.
#> [1] TRUE
max(irep2(unsorted))==reps
#> irep: 15 simulations found.
#> [1] TRUE

Created on 2024-06-24 with reprex v2.0.2

jpryby commented 3 months ago

To add to this, there is no protection in the current version of the function nor this fix that allows the reuse of an ID number across two individuals, even though that is how NONMEM would interpret separate blocks of the same ID number. See below for the simtab output of what I mean (with TABLE breaks for each subproblem); this is two simulations from a dataset containing 3 subjects, as NONMEM would see it.

TABLE 1

ID 1 1 1 2 2 2 1 1 1 | TABLE 2 | ID 1 1 1 2 2 2 1 1 1

To me, this is a rarer edge-case than the one the fix is addressing, but may be worth mentioning in the documentation.