Test/add changing variable type with ~ka

nlmixr2 · Dec 3, 2023 · 9eca019 · 9eca019
1 parent 150b8a6
commit 9eca019
Show file tree

Hide file tree

Showing 4 changed files with 193 additions and 4 deletions.
diff --git a/R/piping-ini.R b/R/piping-ini.R
@@ -408,6 +408,70 @@
   assign("iniDf", ret, envir=rxui)
   invisible()
 }
+#' Handle switching theta to eta and vice versa
+#'
+#' This is coded as model |> ini(~par)
+#'
+#' @param expr Expression, this would be the ~par expression
+#' @param rxui rxui uncompressed environment
+#' @param envir Environment for evaluation (if needed)
+#' @return Nothing, called for side effects
+#' @noRd
+#' @author Matthew L. Fidler
+.iniHandleSwitchType <- function(expr, rxui, envir=parent.frame()) {
+  .var <- as.character(expr[[2]])
+  .iniDf <- rxui$iniDf
+  .w <- which(.iniDf$name == .var)
+  if (length(.w) != 1L) stop("cannot switch parameter type for '", .var, "'", call.=FALSE)
+  .theta <- .iniDf[!is.na(.iniDf$ntheta),, drop = FALSE]
+  .eta <- .iniDf[is.na(.iniDf$ntheta),, drop = FALSE]
+  if (is.na(.iniDf$ntheta[.w])) {
+    # switch eta to theta
+    .neta <- .iniDf$neta1[.w]
+    .eta <- .eta[.eta$neta1 != .neta,, drop = FALSE]
+    .eta <- .eta[.eta$neta2 != .neta,, drop = FALSE]
+    .eta$neta1 <- .eta$neta1 - ifelse(.eta$neta1 < .neta, 0L, 1L)
+    .eta$neta2 <- .eta$neta2 - ifelse(.eta$neta2 < .neta, 0L, 1L)
+    .newTheta <- .iniDf[.w, ]
+    .newTheta$neta1 <- NA_integer_
+    .newTheta$neta2 <- NA_integer_
+    if (length(.theta$ntheta) == 0L) {
+      .newTheta$ntheta <- 1L
+    } else {
+      .newTheta$ntheta <- max(.theta$ntheta) + 1L
+    }
+    .minfo(paste0("convert '", .var, "' from between subject variability to population parameter"))
+    .theta <- rbind(.theta, .newTheta)
+  } else {
+    # switch theta to eta
+    if (!is.na(.iniDf$err[.w])){
+      stop("cannot switch error parameter '", .var,
+           "' to a different type", call. = FALSE)
+    }
+    .ntheta <- .iniDf$ntheta[.w]
+    .theta <- .theta[.theta$ntheta != .ntheta,, drop = FALSE]
+    .theta$ntheta <- .theta$ntheta - ifelse(.theta$ntheta < .ntheta, 0L, 1L)
+    .newEta <- .iniDf[.w, ]
+    .newEta$ntheta <- NA_integer_
+    if (length(.eta$neta1) == 0L) {
+      .newEta$neta1 <- .newEta$neta2 <- 1L
+    } else {
+      .newEta$neta1 <- .newEta$neta2 <- max(.eta$neta1) + 1L
+    }
+    .minfo(paste0("convert '", .var, "' from population parameter to between subject variability"))
+    if (.newEta$est == 0) {
+      .minfo("old initial estimate is zero, changing to 1")
+      .newEta$est <- 1
+    } else if (.newEta$est < 0) {
+      .minfo("old initial estimate was negative, changing to positive")
+      .newEta$est <- -.newEta$est
+    }
+    .eta <- rbind(.eta, .newEta)
+  }
+  .ini <- rbind(.theta, .eta)
+  assign("iniDf", .ini, envir=rxui)
+  invisible()
+}
 
 #' Update the iniDf of a model
 #'
@@ -464,6 +528,8 @@
     expr[[3]] <- eval(as.call(list(quote(`lotri`), as.call(list(quote(`{`), expr)))),
                       envir=envir)[1, 1]
     .iniHandleFixOrUnfixEqual(expr=expr, rxui=rxui, envir=envir, maxLen=1L)
+  } else if (.isTildeExpr(expr)) {
+    .iniHandleSwitchType(expr=expr, rxui=rxui, envir=envir)
   } else {
     # Can this error be improved to clarify what is the expression causing the
     # issue?  It needs a single character string representation of something

diff --git a/R/piping-model.R b/R/piping-model.R
@@ -180,6 +180,10 @@ model.rxModelVars <- model.rxode2
     .matchesLangTemplate(expr, str2lang(". = ."))
 }
 
+.isTildeExpr <- function(expr) {
+  .matchesLangTemplate(expr, str2lang("~ ."))
+}
+
 # get the left hand side of an assignment or endpoint; returns NULL if the input
 # is not an assignment or endpoint
 .getLhs <- function(expr) {

diff --git a/R/piping.R b/R/piping.R
@@ -295,10 +295,12 @@
     } else if (identical(.quoted[[1]], quote(`as.formula`))) {
       .quoted <- .quoted[[2]]
     } else if (identical(.quoted[[1]], quote(`~`))) {
-      .quoted[[3]] <- .iniSimplifyFixUnfix(.quoted[[3]])
-      if (identical(.quoted[[3]], quote(`fix`)) ||
-            identical(.quoted[[3]], quote(`unfix`))) {
-        .quoted <- as.call(list(quote(`<-`), .quoted[[2]], .quoted[[3]]))
+      if (length(.quoted) == 3L) {
+        .quoted[[3]] <- .iniSimplifyFixUnfix(.quoted[[3]])
+        if (identical(.quoted[[3]], quote(`fix`)) ||
+              identical(.quoted[[3]], quote(`unfix`))) {
+          .quoted <- as.call(list(quote(`<-`), .quoted[[2]], .quoted[[3]]))
+        }
       }
     } else if (identical(.quoted[[1]], quote(`$`))) {
       .tmp <- try(eval(.quoted), silent=TRUE)

diff --git a/tests/testthat/test-piping-ini.R b/tests/testthat/test-piping-ini.R
@@ -569,3 +569,120 @@ test_that("append allows promoting from covariate (#472)", {
   )
   expect_equal(newmod$iniDf$name, c("lka", "ka_dose", "lcl", "lvc", "propSd"))
 })
+
+test_that("change ini type with ~", {
+
+  mod <- function() {
+    ini({
+      lka <- 0.45
+      lcl <- 1
+      lvc  <- 3.45
+      propSd <- 0.5
+    })
+    model({
+      ka <- exp(lka)
+      cl <- exp(lcl)
+      vc  <- exp(lvc)
+      kel <- cl / vc
+      d/dt(depot) <- -ka*depot
+      d/dt(central) <- ka*depot-kel*central
+      cp <- central / vc
+      cp ~ prop(propSd)
+    })
+  }
+
+  mod1 <- mod |> ini( ~ lka)
+  expect_equal(mod1$omega, lotri(lka ~ 0.45))
+
+  mod2 <- mod1 |> ini( ~ lka)
+  expect_equal(mod2$omega, NULL)
+
+  expect_error(mod1 |> ini( ~ propSd))
+
+  expect_error(mod1 |> ini( ~ matt))
+
+  ## all etas
+
+  mod <- function() {
+    ini({
+      lka ~ 0.45
+      lcl ~ 1
+      lvc ~ 3.45
+    })
+    model({
+      ka <- exp(lka)
+      cl <- exp(lcl)
+      vc  <- exp(lvc)
+      kel <- cl / vc
+      d/dt(depot) <- -ka*depot
+      d/dt(central) <- ka*depot-kel*central
+      cp <- central / vc
+    })
+  }
+
+  mod2 <- mod |> ini( ~ lka)
+
+  expect_equal(mod2$omega, lotri(lcl ~ 1, lvc ~ 3.45))
+
+  # remove correlated eta
+
+  mod <- function() {
+    ini({
+      lka + lcl + lvc ~
+        c(0.45,
+          0.01, 1,
+          0.01, -0.01, 3.45)
+    })
+    model({
+      ka <- exp(lka)
+      cl <- exp(lcl)
+      vc  <- exp(lvc)
+      kel <- cl / vc
+      d/dt(depot) <- -ka*depot
+      d/dt(central) <- ka*depot-kel*central
+      cp <- central / vc
+    })
+  }
+
+  mod2 <- mod |> ini( ~ lka)
+
+  expect_equal(mod2$omega, lotri(lcl + lvc ~ c(1,
+                                               -0.01, 3.45)))
+
+
+  # negative and zero
+
+  mod <- function() {
+    ini({
+      lka <- 0.45
+      lcl <- -1
+      lvc <- 0
+    })
+    model({
+      ka <- exp(lka)
+      cl <- exp(lcl)
+      vc  <- exp(lvc)
+      kel <- cl / vc
+      d/dt(depot) <- -ka*depot
+      d/dt(central) <- ka*depot-kel*central
+      cp <- central / vc
+    })
+  }
+
+  mod2 <- mod |> ini( ~ lcl)
+
+  expect_equal(mod2$omega, lotri(lcl ~ 1))
+
+  mod2 <- mod |> ini( ~ lvc)
+
+  expect_equal(mod2$omega, lotri(lvc ~ 1))
+
+  mod3 <- mod2 |> ini( ~ lvc)
+
+  expect_equal(mod3$omega, lotri(lcl ~ 1, lvc ~ 1))
+
+  mod4 <- mod3 |> ini( ~ lvc)
+
+  expect_equal(mod4$omega, lotri(lcl ~ 1))
+
+})