From 53bd164676bc3de98c374a22ddbee867aac7bc76 Mon Sep 17 00:00:00 2001 From: Matthew Fidler Date: Tue, 17 Sep 2024 15:16:55 -0500 Subject: [PATCH] Add and test mu3 referencing --- R/mu2.R | 48 +++++++++++++-- R/saemRxUiGetModel.R | 9 ++- tests/testthat/test-saem-mu.R | 106 ++++++++++++++++++++++++++++++++++ 3 files changed, 157 insertions(+), 6 deletions(-) create mode 100644 tests/testthat/test-saem-mu.R diff --git a/R/mu2.R b/R/mu2.R index 0bc467e5..02911783 100644 --- 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,6 +52,15 @@ 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 @@ -49,6 +68,25 @@ mu2env$expit <- rxode2::expit 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, "*", @@ -56,10 +94,9 @@ mu2env$expit <- rxode2::expit .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 index f8463359..1115d7d2 100644 --- 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 new file mode 100644 index 00000000..39f9d2d8 --- /dev/null +++ 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])) + +})