Merge pull request #136 from nlmixr2/135-ex11

Add babelBpopIdx
nlmixr2 · Oct 19, 2024 · 503f39a · 503f39a
2 parents e2beb1b + f8caa3d
commit 503f39a
Show file tree

Hide file tree

Showing 6 changed files with 171 additions and 19 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -144,6 +144,7 @@ export(as.nlmixr)
 export(as.nlmixr2)
 export(as.nonmem2rx)
 export(babel.poped.database)
+export(babelBpopIdx)
 export(bblDatToMonolix)
 export(bblDatToMrgsolve)
 export(bblDatToNonmem)

diff --git a/NEWS.md b/NEWS.md
@@ -13,6 +13,10 @@
   `PopED` code between two different `PopED` databases, like in issue
   #131.
 
+* Added a new function `babelBpopIdx(poped.db, "par")` which will get
+  the poped index for a model generated from `babelmixr2`, which is
+  useful when calculating the power (as in example 11).
+
 # babelmixr2 0.1.4
 
 * Added experimental `PopED` integration

diff --git a/R/poped.R b/R/poped.R
@@ -3042,7 +3042,13 @@ babel.poped.database <- function(popedInput, ..., optTime=NA) {
   if (is.environment(popedInput$babelmixr2)) {
     .babelmixr2 <- popedInput$babelmixr2
     .db <- PopED::create.poped.database(popedInput=popedInput, ...)
-    .db$babelmixr2 <- .babelmixr2
+    .b2 <- new.env(parent=emptyenv())
+    for (v in ls(envir=.babelmixr2, all=TRUE)) {
+      assign(v, get(v, envir=.babelmixr2), envir=.b2)
+    }
+    .b2$modelNumber <- .poped$modelNumber
+    .poped$modelNumber <- .poped$modelNumber + 1
+    .db$babelmixr2 <- .b2
     if (is.logical(optTime) && length(optTime) == 1L &&
           !is.na(optTime)) {
       popedMultipleEndpointResetTimeIndex()
@@ -3054,3 +3060,68 @@ babel.poped.database <- function(popedInput, ..., optTime=NA) {
          call.=FALSE)
   }
 }
+
+#' Get the bpop_idx by variable name for a poped database created by `babelmixr2`
+#'
+#' This may work for other poped databases if the population parameters are named.
+#'
+#' @param popedInput The babaelmixr2 created database
+#' @param var variable to query
+#' @return index of the variable
+#' @export
+#' @author Matthew L. Fidler
+#'
+#' @examples
+#'
+#' if (requireNamespace("PopED", quietly=TRUE)) {
+#'
+#' f <- function() {
+#'   ini({
+#'     tV <- 72.8
+#'     tKa <- 0.25
+#'     tCl <- 3.75
+#'     tF <- fix(0.9)
+#'     pedCL <- 0.8
+#'
+#'     eta.v ~ 0.09
+#'     eta.ka ~ 0.09
+#'     eta.cl ~0.25^2
+#'
+#'     prop.sd <- fix(sqrt(0.04))
+#'     add.sd <- fix(sqrt(5e-6))
+#'
+#'   })
+#'   model({
+#'     V<-tV*exp(eta.v)
+#'     KA<-tKa*exp(eta.ka) * (pedCL**isPediatric) # add covariate for pediatrics
+#'     CL<-tCl*exp(eta.cl)
+#'     Favail <- tF
+#'
+#'     N <-  floor(t/TAU)+1
+#'     y <- (DOSE*Favail/V)*(KA/(KA - CL/V)) *
+#'       (exp(-CL/V * (t - (N - 1) * TAU)) *
+#'          (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) -
+#'          exp(-KA * (t - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU)))
+#'
+#'     y ~ prop(prop.sd) + add(add.sd)
+#'   })
+#' }
+#'
+#' e <- et(c( 1,8,10,240,245))
+#'
+#' babel.db <- nlmixr2(f, e, "poped",
+#'                     popedControl(m = 2,
+#'                                  groupsize=20,
+#'                                  bUseGrouped_xt=TRUE,
+#'                                  a=list(c(DOSE=20,TAU=24,isPediatric = 0),
+#'                                         c(DOSE=40, TAU=24,isPediatric = 0))))
+#'
+#' babelBpopIdx(babel.db, "pedCL")
+#'
+#' }
+babelBpopIdx <- function(popedInput, var)  {
+  .w <- which(rownames(popedInput$parameters$bpop) == var)
+  if (length(.w) == 1L) return(.w)
+  stop("cannot find '", var, "' in the baelmixr2 poped model",
+       call.=FALSE)
+}
diff --git a/man/babelBpopIdx.Rd b/man/babelBpopIdx.Rd
diff --git a/src/poped.cpp b/src/poped.cpp
@@ -218,6 +218,9 @@ Rcpp::NumericVector popedMultipleEndpointParam(Rcpp::NumericVector p,
                                                Rcpp::IntegerVector modelSwitch,
                                                int maxMT,
                                                bool optTime=true) {
+  if (optTime && globalTimeIndexer.isInitialized()) {
+    globalTimeIndexer.reset();
+  }
   globalTimeIndexer.initialize(modelSwitch, times, optTime);
   Rcpp::NumericVector ret(p.size()-1+maxMT);
   std::fill(ret.begin(), ret.end(), globalTimeIndexer.getMaxTime());
@@ -467,7 +470,9 @@ void popedSolveFidMat(arma::mat &matMT, NumericVector &theta, int id, int nrow,
     setIndIdx(ind, j);
     kk = getIndIx(ind, j);
     curT = getTime(kk, ind);
-    isMT = getIndEvid(ind, kk) >= 10 && getIndEvid(ind, kk) <= 99;
+    int evid = getIndEvid(ind, kk);
+    isMT = evid >= 10 && evid <= 99;
+    // REprintf("curT: %f; evid: %d; isMT: %d; kk: %d\n", curT, evid, isMT, kk);
     if (isDose(getIndEvid(ind, kk))) {
       rxInner.calc_lhs(id, curT, getOpIndSolve(op, ind, j), lhs);
       continue;

diff --git a/tests/testthat/test-poped.R b/tests/testthat/test-poped.R
@@ -207,16 +207,16 @@ if (requireNamespace("PopED", quietly=TRUE)) {
 
       library(PopED)
 
-      expect_equal(list(ofv = 20.947993269923,
-                        fim = lotri({
-                          tcl ~ 39.1196937156543
-                          tv ~ c(-0.386317381405573, 40.0037481563051)
-                          sig_var_add.err ~ 800097.977314259
-                        }),
-                        rse = c(tcl = 3.31152092974365,
-                                tv = 30.9526397537534,
-                                sig_var_add.err = 11.1796553130823)),
-                   evaluate_design(db), tolerance = 1e-4)
+      expect_equal(
+        list(ofv = 20.9455595871912,
+             fim = lotri({
+               tcl ~ 39.0235064308582
+               tv ~ c(-0.385137678397696, 40.0037346438455)
+               sig_var_add.err ~ 800120.483866746
+             }),
+             rse = c(tcl = 3.31559905061048, tv = 30.9526395967922,
+                     sig_var_add.err = 11.1794980759702)),
+        evaluate_design(db), tolerance = 1e-4)
 
       ## v <- poped_optim(db, opt_xt=TRUE)
 
@@ -379,27 +379,22 @@ if (requireNamespace("PopED", quietly=TRUE)) {
         tCl <- 3.75
         tF <- fix(0.9)
         pedCL <- 0.8
-
         eta.v ~ 0.09
         eta.ka ~ 0.09
         eta.cl ~0.25^2
-
         prop.sd <- fix(sqrt(0.04))
         add.sd <- fix(sqrt(5e-6))
-
       })
       model({
         V<-tV*exp(eta.v)
-        KA<-tKa*exp(eta.ka) * (pedCL**isPediatric) # add covariate for pediatrics
-        CL<-tCl*exp(eta.cl)
+        KA<-tKa*exp(eta.ka)
+        CL<-tCl*exp(eta.cl)  * (pedCL^isPediatric) # add covariate for pediatrics
         Favail <- tF
-
         N <-  floor(t/TAU)+1
         y <- (DOSE*Favail/V)*(KA/(KA - CL/V)) *
           (exp(-CL/V * (t - (N - 1) * TAU)) *
              (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) -
              exp(-KA * (t - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU)))
-
         y ~ prop(prop.sd) + add(add.sd)
       })
     }
@@ -436,6 +431,11 @@ if (requireNamespace("PopED", quietly=TRUE)) {
     vdiffr::expect_doppelganger("poped-ped-next",
                                 plot_model_prediction(babel.db.ped,
                                                       model_num_points = 300))
+
+    expect_equal(babelBpopIdx(babel.db.ped, "pedCL"), 4L)
+
+    expect_error(babelBpopIdx(babel.db, "matt"))
+
   })
 
   ##  The tests run interactively runs OK