Merge pull request #478 from nlmixr2/478-mu3

`mu` referencing based on derived variables
nlmixr2 · Sep 17, 2024 · 020daf4 · 020daf4
2 parents 1ffa015 + 34bbc11
commit 020daf4
Show file tree

Hide file tree

Showing 4 changed files with 162 additions and 7 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -1,9 +1,13 @@
-# nlmixr2est (development version)
+# nlmixr2est 3.0.0
 
 - No binary linking to `rxode2`, `lbfgsb3c` and `n1q1`, which means
   that updating these will not make `nlmixr2est` crash without
   recompiling.
 
+- New `mu`3 referencing will take context from the model to see if the
+  algebraic expression can be completed from defined model variables;
+  These variable would have to be unique.
+
 # nlmixr2est 2.2.2
 
 ## Breaking changes

diff --git a/R/mu2.R b/R/mu2.R
@@ -27,6 +27,16 @@ mu2env$expit <- rxode2::expit
     expr
   }
 }
+.uiGetMu3 <- function(data, .datEnv, .tmp) {
+  .tmp <- eval(str2lang(paste0("rxode2::rxToSE(", .tmp, ", NULL)")))
+  .tmp <- str2lang(paste0("with(.datEnv$symengine, ", .tmp, ")"))
+  .tmp <- eval(.tmp)
+  .tmp <- as.character(.tmp)
+  .tmp <- str2lang(paste0("rxode2::rxFromSE(", .tmp, ")"))
+  .tmp <- eval(.tmp)
+  .tmp <- str2lang(paste0("with(.datEnv, with(data,",  .tmp, "))"))
+  eval(.tmp)
+}
 #' This function handles mu2 covariates
 #'
 #' In general the dataset is modified with nlmixrMuDerCov# and the mu2
@@ -42,24 +52,51 @@ mu2env$expit <- rxode2::expit
   .datEnv$data <- data
   .datEnv$model <- rxode2::as.model(ui)
   .datEnv$ui  <- ui
+  .datEnv$symengine <- NULL
+  if (use.utf()) {
+    .mu2 <- "\u03BC\u2082"
+    .mu3 <- "\u03BC\u2083"
+  } else {
+    .mu2 <- "mu2"
+    .mu3 <- "mu3"
+  }
+
   lapply(seq_along(ui$mu2RefCovariateReplaceDataFrame$covariate),
          function(i) {
            .datEnv$i <- i
            .tmp <- try(with(.datEnv,
                         with(data,
                              eval(str2lang(ui$mu2RefCovariateReplaceDataFrame$covariate[i])))),
                        silent=TRUE)
+           if (inherits(.tmp, "try-error")) {
+             if (is.null(.datEnv$symengine)) {
+               .minfo(paste0("loading model to look for ", .mu3, "references"))
+               .datEnv$symengine <- ui$loadPruneSaem
+               .minfo("done")
+             }
+             .tmp <- try(.uiGetMu3(data, .datEnv,
+                                   ui$mu2RefCovariateReplaceDataFrame$covariate[i]), silent=TRUE)
+             if (!inherits(.tmp, "try-error")) {
+               .txt <- paste0(.mu3, " item: ", ui$mu2RefCovariateReplaceDataFrame$covariate[i])
+               .minfo(.txt)
+               # Will put into the fit information
+               warning(.txt, call.=FALSE)
+             }
+           } else {
+             .txt <- paste0(.mu2, " item: ", ui$mu2RefCovariateReplaceDataFrame$covariate[i])
+             .minfo(.txt)
+             warning(.txt, call.=FALSE)
+           }
            if (!inherits(.tmp, "try-error")) {
              .datEnv$data[[paste0("nlmixrMuDerCov", i)]] <- .tmp
              .new <- str2lang(paste0("nlmixrMuDerCov", i, "*",
                                      ui$mu2RefCovariateReplaceDataFrame$covariateParameter[i]))
              .old <- str2lang(ui$mu2RefCovariateReplaceDataFrame$modelExpression[i])
              .datEnv$model <- .uiModifyForCovsRep(.datEnv$model, .old, .new)
            } else {
-             warning(paste0("algebraic mu expression failed for '",
-                            ui$mu2RefCovariateReplaceDataFrame$modelExpression[i],
-                            "'"),
-                     call.=FALSE)
+             .txt <- paste0("not ",.mu2," or ", .mu3, " item: ", ui$mu2RefCovariateReplaceDataFrame$covariate[i])
+             .minfo(.txt)
+             warning(.txt, call.=FALSE)
            }
            invisible()
          })
@@ -120,5 +157,8 @@ mu2env$expit <- rxode2::expit
       }
     }
   }
+  # Reset symengine environments
+  .saemModelEnv$symengine <- NULL
+  .saemModelEnv$predSymengine <- NULL
   ret
 }
diff --git a/R/saemRxUiGetModel.R b/R/saemRxUiGetModel.R
@@ -300,6 +300,7 @@ attr(rxUiGet.saemParams, "desc") <- "Get the params() for a saem model"
 #' @export
 rxUiGet.saemModel <- function(x, ...) {
   .s <- rxUiGet.loadPruneSaem(x, ...)
+
   .prd <- get("rx_pred_", envir = .s)
   .prd <- paste0("rx_pred_=", rxode2::rxFromSE(.prd))
   ## .lhs0 <- .s$..lhs0
@@ -381,7 +382,10 @@ rxUiGet.saemModelPredReplaceLst <- function(x, ...) {
 }
 #attr(rxUiGet.saemModelPredReplaceLst, "desc") <- "Replace the mu referenced thetas with these values"
 
-.saemModelPredSymengineEnvironment <- NULL
+.saemModelEnv <- NULL
+.saemModelEnv$symengine <- NULL
+.saemModelEnv$predSymengine <- NULL
+
 
 #' @export
 rxUiGet.interpLinesStr <- function(x, ...) {
@@ -408,8 +412,9 @@ rxUiGet.saemModelPred <- function(x, ...) {
     .levels <- paste(.levels, collapse="\n")
   }
   .s <- rxUiGet.loadPruneSaemPred(x, ...)
+  .saemModelEnv$symengine <- .s
   .replaceLst <- rxUiGet.saemModelPredReplaceLst(x, ...)
-  assignInMyNamespace(".saemModelPredSymengineEnvironment", .s)
+  .saemModelEnv$predSymengine <- .s
   .prd <- get("rx_pred_", envir = .s)
   .prd <- paste0("rx_pred_=", rxode2::rxFromSE(.prd))
   .r <- get("rx_r_", envir = .s)

diff --git a/tests/testthat/test-saem-mu.R b/tests/testthat/test-saem-mu.R
@@ -0,0 +1,106 @@
+test_that("saem mu reference 1", {
+
+  .nlmixr2 <- function(...) { suppressMessages(nlmixr2est::nlmixr2(...)) }
+
+  theo_sd2 <- nlmixr2data::theo_sd
+
+  theo_sd2$lwt<-log(theo_sd2$WT/70)
+
+  # The basic model consists of an ini block that has initial estimates
+  # Original mu-referencing
+  one.compartment <- function() {
+    ini({
+      tka <- log(1.57); label("Ka")
+      tcl <- log(2.72); label("Cl")
+      tv <- log(31.5); label("V")
+      covwt<- 0.01
+      eta.ka ~ 0.6
+      eta.cl ~ 0.3
+      eta.v ~ 0.1
+      add.sd <- 0.7
+    })
+    model({
+      ka <- exp(tka + eta.ka)
+      cl <- exp(tcl + eta.cl + covwt*lwt)
+      v <- exp(tv + eta.v)
+      d/dt(depot) <- -ka * depot
+      d/dt(center) <- ka * depot - cl / v * center
+      cp <- center / v
+      cp ~ add(add.sd)
+    })
+  }
+
+  fit1 <- .nlmixr2(one.compartment, theo_sd2,  est="saem",
+                  control=saemControl(print=0,seed = 1234, nBurn = 1, nEm = 1,
+                                      calcTables = FALSE))
+
+  # true mu expression should not have information in $runInfo
+  expect_true(is.null(fit1$runInfo))
+
+  # mu2-referencing
+  theo_sd2 <- nlmixr2data::theo_sd
+  one.compartment <- function() {
+    ini({
+      tka <- log(1.57); label("Ka")
+      tcl <- log(2.72); label("Cl")
+      tv <- log(31.5); label("V")
+      covwt <- 0.01
+      eta.ka ~ 0.6
+      eta.cl ~ 0.3
+      eta.v ~ 0.1
+      add.sd <- 0.7
+    })
+    model({
+      ka <- exp(tka + eta.ka)
+      cl <- exp(tcl + eta.cl + log(WT/70)*covwt)
+      v <- exp(tv + eta.v)
+      d/dt(depot) <- -ka * depot
+      d/dt(center) <- ka * depot - cl / v * center
+      cp <- center / v
+      cp ~ add(add.sd)
+    })
+  }
+
+
+  fit2 <-
+    withr::with_options(list(cli.unicode=FALSE),{
+      .nlmixr2(one.compartment, theo_sd2,  est="saem",
+              control=saemControl(print=0,seed = 1234, nBurn = 1, nEm = 1,
+                                  calcTables = FALSE))
+    })
+
+  expect_true(grepl("mu2 item:", fit2$runInfo[1]))
+
+  # mu3 referencing
+  one.compartment <- function() {
+    ini({
+      tka <- log(1.57); label("Ka")
+      tcl <- log(2.72); label("Cl")
+      tv <- log(31.5); label("V")
+      covwt <- 0.01
+      eta.ka ~ 0.6
+      eta.cl ~ 0.3
+      eta.v ~ 0.1
+      add.sd <- 0.7
+    })
+    model({
+      ka <- exp(tka + eta.ka)
+      wt70 <- WT/70
+      cl <- exp(tcl + eta.cl + log(wt70)*covwt)
+      v <- exp(tv + eta.v)
+      d/dt(depot) <- -ka * depot
+      d/dt(center) <- ka * depot - cl / v * center
+      cp <- center / v
+      cp ~ add(add.sd)
+    })
+  }
+
+  fit3 <-   withr::with_options(list(cli.unicode=FALSE), {
+    .nlmixr2(one.compartment, theo_sd2,  est="saem",
+            control=saemControl(print=0,seed = 1234, nBurn = 1, nEm = 1,
+                                calcTables = FALSE))
+  })
+
+  expect_true(grepl("mu3 item", fit3$runInfo[1]))
+
+})