diff --git a/NAMESPACE b/NAMESPACE index 2a66d8860..77cfd01a3 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -355,6 +355,7 @@ export(rxExpandGrid_) export(rxExpandIfElse) export(rxExpandSens2_) export(rxExpandSens_) +export(rxFixPop) export(rxForget) export(rxFromSE) export(rxFun) diff --git a/NEWS.md b/NEWS.md index 545036d89..f600f968a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -13,6 +13,9 @@ - Create a function to see if a rxode2 solve is loaded in memory (`rxode2::rxSolveSetup()`) +- Create a new function that fixes the rxode2 population values in the + model (and drops them in the initial estimates); `rxFixPop()` + # rxode2 2.1.2 ## Other changes diff --git a/R/ui-fix.R b/R/ui-fix.R new file mode 100644 index 000000000..aac9198a1 --- /dev/null +++ b/R/ui-fix.R @@ -0,0 +1,137 @@ +#' Hard fix population variables +#' +#' @param item expression to consider +#' @param var variable list to fix +#' @param isLhs is the expression a lhs expression? +#' @return expression with variables replaced with constants +#' @noRd +#' @author Matthew L. Fidler +.rxFixPopVar <- function(item, var, isLhs=FALSE) { + if (is.atomic(item)) { + return(item) + } + if (is.name(item)) { + .n <- as.character(item) + .n <- var[.n] + if (!is.na(.n)) { + return(setNames(.n, NULL)) + } + return(item) + } else if (is.call(item)) { + if (isLhs && identical(item[[1]], quote(`/`))) { + # handle d/dt() differently so that d doesn't get renamed + .num <- item[[2]] + .denom <- item[[3]] + if (is.call(.num)) .num <- as.call(lapply(.num, .rxFixPopVar, var=var, isLhs=TRUE)) + if (is.call(.denom)) .denom <- as.call(lapply(.denom, .rxFixPopVar, var=var, isLhs=TRUE)) + return(as.call(c(list(item[[1]]), .num, .denom))) + } else if (isLhs && length(item) == 2L && + is.numeric(item[[2]])) { + .env <- new.env(parent=emptyenv()) + .env$new <- NULL + lapply(seq_along(var), + function(i) { + if (!is.null(.env$new)) return(NULL) + .curVal <- setNames(var[i], NULL) + .old <- str2lang(names(var[i])) + if (identical(item[[1]], .old)) { + .env$new <- .curVal + } + return(NULL) + }) + if (!is.null(.env$new)) { + # handle x(0) = items + return(as.call(c(.env$new, lapply(item[-1], .rxFixPopVar, var=var, isLhs=isLhs)))) + } + } + if (identical(item[[1]], quote(`=`)) || + identical(item[[1]], quote(`<-`)) || + identical(item[[1]], quote(`~`))) { + .elhs <- lapply(item[c(-1, -3)], .rxFixPopVar, var=var, isLhs=TRUE) + .erhs <- lapply(item[c(-1, -2)], .rxFixPopVar, var=var, isLhs=FALSE) + return(as.call(c(item[[1]], .elhs, .erhs))) + } else { + return(as.call(c(list(item[[1]]), lapply(item[-1], .rxFixPopVar, var=var, isLhs=isLhs)))) + } + } else { + stop("unknown expression", call.=FALSE) + } +} + +#' Apply the fixed population estimated parameters +#' +#' @param ui rxode2 ui function +#' @param returnNull boolean for if unchanged values should return the +#' original ui (`FALSE`) or null (`TRUE`) +#' @return when `returnNull` is TRUE, NULL if nothing was changed, or +#' the changed model ui. When `returnNull` is FALSE, return a ui no +#' matter if it is changed or not. +#' @export +#' @author Matthew L. Fidler +#' @examples +#' +#' One.comp.transit.allo <- function() { +#' ini({ +#' # Where initial conditions/variables are specified +#' lktr <- log(1.15) #log k transit (/h) +#' lcl <- log(0.15) #log Cl (L/hr) +#' lv <- log(7) #log V (L) +#' ALLC <- fix(0.75) #allometric exponent cl +#' ALLV <- fix(1.00) #allometric exponent v +#' prop.err <- 0.15 #proportional error (SD/mean) +#' add.err <- 0.6 #additive error (mg/L) +#' eta.ktr ~ 0.5 +#' eta.cl ~ 0.1 +#' eta.v ~ 0.1 +#' }) +#' model({ +#' #Allometric scaling on weight +#' cl <- exp(lcl + eta.cl + ALLC * logWT70) +#' v <- exp(lv + eta.v + ALLV * logWT70) +#' ktr <- exp(lktr + eta.ktr) +#' # RxODE-style differential equations are supported +#' d/dt(depot) = -ktr * depot +#' d/dt(central) = ktr * trans - (cl/v) * central +#' d/dt(trans) = ktr * depot - ktr * trans +#' ## Concentration is calculated +#' cp = central/v +#' # And is assumed to follow proportional and additive error +#' cp ~ prop(prop.err) + add(add.err) +#' }) +#' } +#' +#' m <- rxFixPop(One.comp.transit.allo) +#' m +#' +#' # now everything is already fixed, so calling again will do nothing +#' +#' rxFixPop(m) +#' +#' # if you call it with returnNull=TRUE when no changes have been +#' # performed, the function will return NULL +#' +#' rxFixPop(m, returnNull=TRUE) +#' +rxFixPop <- function(ui, returnNull=FALSE) { + checkmate::assertLogical(returnNull, any.missing = FALSE, len=1, null.ok=FALSE) + .model <- rxUiDecompress(assertRxUi(ui)) + .model <- .copyUi(.model) + .iniDf <- .model$iniDf + .w <- which(!is.na(.iniDf$ntheta) & is.na(.iniDf$err) & .iniDf$fix) + if (length(.w) == 0L) { + if (returnNull) return(NULL) + return(.model) + } + .v <- setNames(.iniDf$est[.w], .iniDf$name[.w]) + .lst <- lapply(.model$lstExpr, + function(e) { + .rxFixPopVar(e, .v) + }) + .iniDf <- .iniDf[-.w, ] + .iniDf$ntheta <- ifelse(is.na(.iniDf$ntheta), NA_integer_, seq_along(.iniDf$ntheta)) + assign("iniDf", .iniDf, envir=.model) + suppressMessages({ + model(.model) <- .lst + .model + }) +} diff --git a/_pkgdown.yml b/_pkgdown.yml index da8be4eda..198af0332 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -38,6 +38,7 @@ reference: - zeroRe - assertRxUi - rxAppendModel + - rxFixPop - rxRename - update.rxUi - as.rxUi diff --git a/data/rxReservedKeywords.rda b/data/rxReservedKeywords.rda index 80fc13cd1..ddeded968 100644 Binary files a/data/rxReservedKeywords.rda and b/data/rxReservedKeywords.rda differ diff --git a/data/rxResidualError.rda b/data/rxResidualError.rda index df63ee3a5..b5245bbf5 100644 Binary files a/data/rxResidualError.rda and b/data/rxResidualError.rda differ diff --git a/data/rxSyntaxFunctions.rda b/data/rxSyntaxFunctions.rda index bb46aebcc..f2c016e59 100644 Binary files a/data/rxSyntaxFunctions.rda and b/data/rxSyntaxFunctions.rda differ diff --git a/man/reexports.Rd b/man/reexports.Rd index d6ab11a0c..1086b491a 100644 --- a/man/reexports.Rd +++ b/man/reexports.Rd @@ -83,4 +83,3 @@ below to see their documentation. \item{rxode2random}{\code{\link[rxode2random:dot-cbindOme]{.cbindOme}}, \code{\link[rxode2random:dot-expandPars]{.expandPars}}, \code{\link[rxode2random:dot-vecDf]{.vecDf}}, \code{\link[rxode2random]{cvPost}}, \code{\link[rxode2random]{invWR1d}}, \code{\link[rxode2random]{phi}}, \code{\link[rxode2random]{rinvchisq}}, \code{\link[rxode2random]{rLKJ1}}, \code{\link[rxode2random]{rxGetSeed}}, \code{\link[rxode2random]{rxGetSeed}}, \code{\link[rxode2random]{rxRmvn}}, \code{\link[rxode2random]{rxSeedEng}}, \code{\link[rxode2random]{rxSetSeed}}, \code{\link[rxode2random]{rxSetSeed}}, \code{\link[rxode2random]{rxSetSeed}}, \code{\link[rxode2random:rxWithSeed]{rxWithPreserveSeed}}, \code{\link[rxode2random]{rxWithSeed}}, \code{\link[rxode2random]{rxWithSeed}}} }} -\value{ Inherited from parent routine } diff --git a/man/rxFixPop.Rd b/man/rxFixPop.Rd new file mode 100644 index 000000000..e0791249e --- /dev/null +++ b/man/rxFixPop.Rd @@ -0,0 +1,70 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ui-fix.R +\name{rxFixPop} +\alias{rxFixPop} +\title{Apply the fixed population estimated parameters} +\usage{ +rxFixPop(ui, returnNull = FALSE) +} +\arguments{ +\item{ui}{rxode2 ui function} + +\item{returnNull}{boolean for if unchanged values should return the +original ui (\code{FALSE}) or null (\code{TRUE})} +} +\value{ +when \code{returnNull} is TRUE, NULL if nothing was changed, or +the changed model ui. When \code{returnNull} is FALSE, return a ui no +matter if it is changed or not. +} +\description{ +Apply the fixed population estimated parameters +} +\examples{ + +One.comp.transit.allo <- function() { + ini({ + # Where initial conditions/variables are specified + lktr <- log(1.15) #log k transit (/h) + lcl <- log(0.15) #log Cl (L/hr) + lv <- log(7) #log V (L) + ALLC <- fix(0.75) #allometric exponent cl + ALLV <- fix(1.00) #allometric exponent v + prop.err <- 0.15 #proportional error (SD/mean) + add.err <- 0.6 #additive error (mg/L) + eta.ktr ~ 0.5 + eta.cl ~ 0.1 + eta.v ~ 0.1 + }) + model({ + #Allometric scaling on weight + cl <- exp(lcl + eta.cl + ALLC * logWT70) + v <- exp(lv + eta.v + ALLV * logWT70) + ktr <- exp(lktr + eta.ktr) + # RxODE-style differential equations are supported + d/dt(depot) = -ktr * depot + d/dt(central) = ktr * trans - (cl/v) * central + d/dt(trans) = ktr * depot - ktr * trans + ## Concentration is calculated + cp = central/v + # And is assumed to follow proportional and additive error + cp ~ prop(prop.err) + add(add.err) + }) +} + +m <- rxFixPop(One.comp.transit.allo) +m + +# now everything is already fixed, so calling again will do nothing + +rxFixPop(m) + +# if you call it with returnNull=TRUE when no changes have been +# performed, the function will return NULL + +rxFixPop(m, returnNull=TRUE) + +} +\author{ +Matthew L. Fidler +} diff --git a/man/rxode2.Rd b/man/rxode2.Rd index d9b9ea79e..41c77ad83 100644 --- a/man/rxode2.Rd +++ b/man/rxode2.Rd @@ -238,7 +238,7 @@ mod <- rxode2(\{ ## from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2parse.h:33, ## from /home/matt/src/rxode2/inst/include/rxode2.h:9, ## from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2_model_shared.h:3, -## from rx_80ab028288eddd16733200578a7fac4b_.c:117: +## from rx_36fd6c4312c660ddf5390102e91a4bb3_.c:117: ## /usr/share/R/include/R_ext/Complex.h:80:6: warning: ISO C99 doesn’t support unnamed structs/unions [-Wpedantic] ## 80 | \}; ## | ^ @@ -264,7 +264,7 @@ mod <- rxode2(\{ ## from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2parse.h:33, ## from /home/matt/src/rxode2/inst/include/rxode2.h:9, ## from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2_model_shared.h:3, -## from rx_16fc28e9f2ad308e65c8fb7a9b53fdc1_.c:117: +## from rx_4fdd80e5dd5fd33d1951b83e06f7ad84_.c:117: ## /usr/share/R/include/R_ext/Complex.h:80:6: warning: ISO C99 doesn’t support unnamed structs/unions [-Wpedantic] ## 80 | \}; ## | ^ @@ -352,12 +352,12 @@ compilation model. ## from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2parse.h:33, ## from /home/matt/src/rxode2/inst/include/rxode2.h:9, ## from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2_model_shared.h:3, -## from rx_e3306ec84cb8151ac51b0c27ef3dbbe7_.c:117: +## from rx_bbee652505e8bddf28c7e688baef12d9_.c:117: ## /usr/share/R/include/R_ext/Complex.h:80:6: warning: ISO C99 doesn’t support unnamed structs/unions [-Wpedantic] ## 80 | \}; ## | ^ -## rxode2 2.1.2 model named rx_e3306ec84cb8151ac51b0c27ef3dbbe7 model (ready). +## rxode2 2.1.2.9000 model named rx_bbee652505e8bddf28c7e688baef12d9 model (ready). ## x$state: depot, center ## x$stateExtra: cp ## x$params: tka, tcl, tv, add.sd, eta.ka, eta.cl, eta.v, rxerr.cp @@ -373,12 +373,12 @@ mod$simulationIniModel ## from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2parse.h:33, ## from /home/matt/src/rxode2/inst/include/rxode2.h:9, ## from /home/matt/R/x86_64-pc-linux-gnu-library/4.3/rxode2parse/include/rxode2_model_shared.h:3, -## from rx_6f91c3aeabfdac8ced143f492e648867_.c:117: +## from rx_fe01e26500a95dd1ba4d8972bc71bc7a_.c:117: ## /usr/share/R/include/R_ext/Complex.h:80:6: warning: ISO C99 doesn’t support unnamed structs/unions [-Wpedantic] ## 80 | \}; ## | ^ -## rxode2 2.1.2 model named rx_6f91c3aeabfdac8ced143f492e648867 model (ready). +## rxode2 2.1.2.9000 model named rx_fe01e26500a95dd1ba4d8972bc71bc7a model (ready). ## x$state: depot, center ## x$stateExtra: cp ## x$params: tka, tcl, tv, add.sd, eta.ka, eta.cl, eta.v, rxerr.cp diff --git a/tests/testthat/test-rxFix.R b/tests/testthat/test-rxFix.R new file mode 100644 index 000000000..5afe3e6fb --- /dev/null +++ b/tests/testthat/test-rxFix.R @@ -0,0 +1,133 @@ +test_that("rxFix", { + + One.comp.transit.allo <- function() { + ini({ + # Where initial conditions/variables are specified + lktr <- log(1.15) #log k transit (/h) + lcl <- log(0.15) #log Cl (L/hr) + lv <- log(7) #log V (L) + ALLC <- fix(0.75) #allometric exponent cl + ALLV <- fix(1.00) #allometric exponent v + prop.err <- 0.15 #proportional error (SD/mean) + add.err <- 0.6 #additive error (mg/L) + eta.ktr ~ 0.5 + eta.cl ~ 0.1 + eta.v ~ 0.1 + }) + model({ + #Allometric scaling on weight + cl <- exp(lcl + eta.cl + ALLC * logWT70) + v <- exp(lv + eta.v + ALLV * logWT70) + ktr <- exp(lktr + eta.ktr) + # RxODE-style differential equations are supported + d/dt(depot) = -ktr * depot + d/dt(central) = ktr * trans - (cl/v) * central + d/dt(trans) = ktr * depot - ktr * trans + ## Concentration is calculated + cp = central/v + # And is assumed to follow proportional and additive error + cp ~ prop(prop.err) + add(add.err) + }) + } + + tmp <- rxFixPop(One.comp.transit.allo) + + expect_equal(tmp$theta, + c(lktr = 0.139761942375159, lcl = -1.89711998488588, lv = 1.94591014905531, + prop.err = 0.15, add.err = 0.6)) + + tmp2 <- rxFixPop(tmp) + + expect_equal(tmp2$theta, + c(lktr = 0.139761942375159, lcl = -1.89711998488588, lv = 1.94591014905531, + prop.err = 0.15, add.err = 0.6)) + + tmp3 <- rxFixPop(tmp2, returnNull = TRUE) + + expect_true(is.null(tmp3)) + + nlmixr_threecmt_mm_no_add_wtcl_pdtg_kout_delay2 <- function() { + ini({ + tf_sc <- log(999) + tf_infilt <- log(999) + tka_sc <- log(999) + tka_infilt <- log(999) + tcl_low <- log(999) + tcl_high <- log(999) + tcl_c50 <- log(3000) + e_wt_cl <- fixed(999) + tv <- log(999) + tq1 <- log(999) + tvp1 <- log(10) + tq2 <- log(999) + tvp2 <- log(20) + eta_cl~999 + eta_v~999 + prop_err <- 999 + tg_bl <- log(999) + eta_tg_bl~999 + tg_kel <- log(999) + tg_ec50 <- log(5000) + tg_emax_kel <- log(2) + ktr_tg <- log(999) + prop_err_tg <- 999 + }) + model({ + # PK setup + f_sc <- exp(tf_sc) + f_infilt <- exp(tf_infilt) + ka_sc <- exp(tka_sc) + ka_infilt <- exp(tka_infilt) + cl_low <- exp(tcl_low + eta_cl)*(WEIGHT_BL/85)^e_wt_cl + cl_high <- exp(tcl_high + eta_cl)*(WEIGHT_BL/85)^e_wt_cl + cl_c50 <- exp(tcl_c50) + v <- exp(tv + eta_v) + q1 <- exp(tq1) + vp1 <- exp(tvp1) + q2 <- exp(tq2) + vp2 <- exp(tvp2) + # PK micro-parameters + ke_low <- cl_low/v + ke_high <- cl_high/v + kc_p1 <- q1/v + kp1_c <- q1/vp1 + kc_p2 <- q2/v + kp2_c <- q2/vp2 + # TG setup + tgbl <- exp(tg_bl + eta_tg_bl) + kin_tg <- tgbl*exp(tg_kel) + ktr_TG <- exp(ktr_tg) + TG(0) <- tgbl + # differential equations + cp <- CENTRAL/v*1e3 # 1e3 is for unit conversion + ke <- ke_low + (ke_high - ke_low)*cp/(cp + cl_c50) + kout_tg <- exp(tg_kel) + exp(tg_emax_kel)*TG_TR/(TG_TR + exp(tg_ec50)) + d/dt(IVINFILT) = - ka_infilt * IVINFILT + d/dt(SC) = -ka_sc * SC + d/dt(CENTRAL) = ka_sc * SC + ka_infilt * IVINFILT - ke*CENTRAL - kc_p1*CENTRAL + kp1_c*P1 - kc_p2*CENTRAL + kp2_c*P2 + d/dt(P1) = kc_p1*CENTRAL - kp1_c*P1 + d/dt(P2) = kc_p2*CENTRAL - kp2_c*P2 + f(SC) <- f_sc + f(IVINFILT) <- f_infilt + # TG transit model + d/dt(TG_TR) = ktr_tg*cp - ktr_tg*TG_TR + d/dt(TG) = kin_tg - kout_tg*TG + # Residual error models + cp ~ prop(prop_err) + TG ~ prop(prop_err_tg) + }) + } + + tmp <- rxFixPop(nlmixr_threecmt_mm_no_add_wtcl_pdtg_kout_delay2) + + expect_equal(tmp$theta, + c(tf_sc = 6.90675477864855, tf_infilt = 6.90675477864855, tka_sc = 6.90675477864855, + tka_infilt = 6.90675477864855, tcl_low = 6.90675477864855, tcl_high = 6.90675477864855, + tcl_c50 = 8.00636756765025, tv = 6.90675477864855, tq1 = 6.90675477864855, + tvp1 = 2.30258509299405, tq2 = 6.90675477864855, tvp2 = 2.99573227355399, + prop_err = 999, tg_bl = 6.90675477864855, tg_kel = 6.90675477864855, + tg_ec50 = 8.51719319141624, tg_emax_kel = 0.693147180559945, + ktr_tg = 6.90675477864855, prop_err_tg = 999)) + + +})