Merge pull request #352 from nlmixr2/augPred-fix

Aug pred fix
nlmixr2 · May 24, 2023 · 24a372e · 24a372e
2 parents 3887055 + b9024f2
commit 24a372e
Show file tree

Hide file tree

Showing 5 changed files with 70 additions and 17 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -65,15 +65,14 @@ Imports:
     rxode2 (>= 2.0.12),
     stats,
     symengine,
-    ucminf,
     utils
 Suggests: 
     broom.mixed,
     crayon,
     data.table,
     devtools,
     digest,
-    dplyr,
+    dplyr (>= 1.1.0),
     generics,
     nloptr,
     qs,

diff --git a/NEWS.md b/NEWS.md
@@ -12,6 +12,9 @@
 
 - `vpcSim()` works when an eta value is fixed to 0 (#341)
 
+- `augPred()` now consistently uses the simulation model (instead of
+  the inner model used for `CWRES` calculation).
+
 # nlmixr2est 2.1.5
 
 - Add `$fitMergeFull`, `$fitMergInner`, `$fitMergeLeft`,

diff --git a/R/augPred.R b/R/augPred.R
@@ -1,6 +1,11 @@
+.augPredIpredModel <- function(fit) {
+  .ipredModel <- .getSimModel(fit, hideIpred=FALSE,tad=FALSE)
+  eval(as.call(list(quote(`rxModelVars`), .ipredModel[[-1]])))
+}
+
 .augPredExpandData <- function(fit, covsInterpolation = c("locf", "nocb", "linear", "midpoint"),
                                minimum = NULL, maximum = NULL, length.out = 51L) {
-  .origData <- rxode2::etTrans(fit$dataSav, fit$ipredModel, addCmt=TRUE, keepDosingOnly=TRUE, allTimeVar=TRUE)
+  .origData <- rxode2::etTrans(fit$dataSav, .augPredIpredModel(fit), addCmt=TRUE, keepDosingOnly=TRUE, allTimeVar=TRUE)
   .predDf <- fit$ui$predDf
   .range <- range(.origData$TIME)
   .covs <- fit$ui$allCovs
@@ -86,6 +91,7 @@ nlmixr2AugPredSolve <- function(fit, covsInterpolation = c("locf", "nocb", "line
   .events <- .augPredExpandData(fit, covsInterpolation = covsInterpolation,
                                 minimum = minimum, maximum = maximum,
                                 length.out = length.out)
+
   # ipred
   .sim <- rxode2::rxSolve(object=.rx, .params, .events,
                           keep=c("DV", "CMT"), returnType="data.frame")
@@ -109,7 +115,16 @@ nlmixr2AugPredSolve <- function(fit, covsInterpolation = c("locf", "nocb", "line
   .stk$id <- .sim$id
   .stk$time <- .sim$time
   .stk$cmt <- as.integer(.sim$CMT)
-  levels(.stk$cmt) <- c(fit$ipredModel$state, fit$ipredModel$stateExtra)
+  .ipredModel <- .augPredIpredModel(fit)
+  .lvl <- c(.ipredModel$state, .ipredModel$stateExtra)
+  if (length(.lvl) == 1L && .lvl == "rxLinCmt") {
+    if (rxModelVars(fit)$flags["ka"] == c(ka=1L)) {
+      .lvl <- c("depot", "central")
+    } else {
+      .lvl <- "central"
+    }
+  }
+  levels(.stk$cmt) <- .lvl
   class(.stk$cmt) <- "factor"
   .stk <- .stk[!is.na(.stk$values), ]
   class(.stk) <- c("nlmixr2AugPred", "data.frame")

diff --git a/R/focei.R b/R/focei.R
@@ -26,18 +26,6 @@ is.latex <- function() {
   get("is_latex_output", asNamespace("knitr"))()
 }
 
-
-.ucminf <- function(par, fn, gr, lower = -Inf, upper = Inf, control = list(), ...) {
-  rxode2::rxReq("ucminf")
-  .ctl <- control
-  .ctl$stepmax <- control$rhobeg
-  .ctl$maxeval <- control$maxOuterIterations
-  .ctl <- .ctl[names(.ctl) %in% c("stepmax", "maxeval")]
-  .ret <- ucminf::ucminf(par, fn, gr = NULL, ..., control = list(), hessian = 2)
-  .ret$x <- .ret$par
-  return(.ret)
-}
-
 .bobyqa <- function(par, fn, gr, lower = -Inf, upper = Inf, control = list(), ...) {
   .ctl <- control
   if (is.null(.ctl$npt)) .ctl$npt <- length(par) * 2 + 1

diff --git a/tests/testthat/test-augpred.R b/tests/testthat/test-augpred.R
@@ -84,7 +84,7 @@ nmTest({
                       lapply(doses, function(x) {
                         ids <- dat %>%
                           dplyr::filter(DOSE == x) %>%
-                          dplyr::summarize(ids=unique(ID)) %>%
+                          dplyr::reframe(ids=unique(ID)) %>%
                           dplyr::pull()
                         ids <- ids[seq(1, nid)]
                         dat %>%
@@ -143,4 +143,52 @@ nmTest({
 
     expect_error(augPred(fit2), NA)
   })
+
+  test_that("mixed pkpd with effect compartment augpred", {
+
+    dat <- nlmixr2data::warfarin
+
+    mod <- function () {
+      ini({
+        tktr <- -0.0407039444259225
+        tcl <- -1.94598426244892
+        tv <- 2.09185800199064
+        eps.pkprop <- c(0, 0.103324277742915)
+        eps.pkadd <- c(0, 0.41909330805883)
+        tc50 <- 0.121920646717701
+        tkout <- -3.1790123282253
+        te0 <- 4.43865746919035
+        eps.pdadd <- c(0, 5.96046447753906e-07)
+        eta.ktr ~ 0.625781701127507
+        eta.cl ~ 0.102936117458982
+        eta.v ~ 0.0491541297805943
+        eta.c50 ~ 0.41589473728507
+        eta.kout ~ 0.108340169259417
+        eta.e0 ~ 0.0632084214464694
+      })
+      model({
+        ktr <- exp(tktr + eta.ktr)
+        cl <- exp(tcl + eta.cl)
+        v <- exp(tv + eta.v)
+        c50 <- exp(tc50 + eta.c50)
+        kout <- exp(tkout + eta.kout)
+        e0 <- exp(te0 + eta.e0)
+        emax <- 1
+        cp <- central/v
+        d/dt(depot) <- -ktr * depot
+        d/dt(central) <- ktr * trans - cl * cp
+        d/dt(trans) <- ktr * depot - ktr * trans
+        d/dt(ce) = kout * (cp - ce)
+        effect <- e0 * (1 - emax * ce/(c50 + ce))
+        cp ~ prop(eps.pkprop) + add(eps.pkadd) | cp
+        effect ~ add(eps.pdadd) | pca
+      })
+    }
+
+    fit <- nlmixr2(mod, dat, "posthoc")
+
+    expect_error(augPred(fit), NA)
+
+
+  })
 })