diff --git a/NAMESPACE b/NAMESPACE index a3797a01b..d3f312655 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -474,6 +474,7 @@ export(scale_y_discrete) export(setRxThreads) export(stat_amt) export(stat_cens) +export(toTrialDuration) export(uppergamma) export(waiver) export(write.template.server) @@ -534,6 +535,7 @@ importFrom(rxode2et,rxEvid) importFrom(rxode2et,rxRateDur) importFrom(rxode2et,rxReq) importFrom(rxode2et,rxStack) +importFrom(rxode2et,toTrialDuration) importFrom(rxode2parse,.getLastIdLvl) importFrom(rxode2parse,forderForceBase) importFrom(rxode2parse,rxDerived) diff --git a/NEWS.md b/NEWS.md index 8a70b08dc..607dab26f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -169,6 +169,8 @@ mu-referencing style to run the optimization. `ini(param=NULL)` changes the parameter to a covariate to align with this idiom of dropping parameters +- `rxRename` has been refactored to run faster + ## Internal new features - Add `as.model()` for list expressions, which implies `model(ui) <- @@ -188,6 +190,9 @@ mu-referencing style to run the optimization. ## Bug fixes +- Simulating/solving from functions/ui now prefers params over `omega` + and `sigma` in the model (#632) + - Piping does not add constants to the initial estimates - When constants are specified in the `model({})` block (like `k <- 1`), they will not diff --git a/R/piping-ini.R b/R/piping-ini.R index 7bdbc8627..eb9802ead 100644 --- a/R/piping-ini.R +++ b/R/piping-ini.R @@ -402,18 +402,18 @@ # Do nothing return() } else if (is.logical(append)) { - checkmate::assert_logical(append, any.missing = FALSE, len = 1) + checkmate::assertLogical(append, any.missing = FALSE, len = 1) if (isTRUE(append)) { appendClean <- Inf } else if (isFALSE(append)) { appendClean <- 0 } } else if (is.numeric(append)) { - checkmate::assert_number(append, null.ok = FALSE, na.ok = FALSE) + checkmate::assertNumber(append, null.ok = FALSE, na.ok = FALSE) appendClean <- append } else if (is.character(append)) { - checkmate::assert_character(append, any.missing = FALSE, len = 1, null.ok = FALSE) - checkmate::assert_choice(append, choices = ini$name) + checkmate::assertCharacter(append, any.missing = FALSE, len = 1, null.ok = FALSE) + checkmate::assertChoice(append, choices = ini$name) appendClean <- which(ini$name == append) } else { stop("'append' must be NULL, logical, numeric, or character/expression of variable in model", diff --git a/R/plot.R b/R/plot.R index 2bb1f6175..f6cc66a2d 100644 --- a/R/plot.R +++ b/R/plot.R @@ -176,8 +176,8 @@ rxTheme <- function(base_size = 11, base_family = "", stopifnot(length(log) == 1) stopifnot(is.character(log)) stopifnot(log %in% c("", "x", "y", "xy", "yx")) - useLogX <- nchar(log) == 2 | log == "x" - useLogY <- nchar(log) == 2 | log == "y" + useLogX <- nchar(log) == 2L || log == "x" + useLogY <- nchar(log) == 2L || log == "y" useXgxr <- getOption("rxode2.xgxr", TRUE) && requireNamespace("xgxr", quietly = TRUE) diff --git a/R/reexport.R b/R/reexport.R index 3dd96e1a6..3b65d4e47 100644 --- a/R/reexport.R +++ b/R/reexport.R @@ -272,3 +272,7 @@ rxode2parse::.getLastIdLvl #' @importFrom rxode2random .expandPars #' @export rxode2random::.expandPars + +#' @importFrom rxode2et toTrialDuration +#' @export +rxode2et::toTrialDuration diff --git a/R/rxsolve.R b/R/rxsolve.R index 6af91c0e6..f6aacc4a1 100644 --- a/R/rxsolve.R +++ b/R/rxsolve.R @@ -1157,16 +1157,17 @@ rxSolve.function <- function(object, params = NULL, events = NULL, inits = NULL, .lst <- list(...) .nlst <- names(.lst) .w <- which(vapply(names(.ctl), function(x) { - !(x %in% names(.nlst)) && exists(x, envir=.meta) + !(x %in% .nlst) && exists(x, envir=.meta) }, logical(1), USE.NAMES=FALSE)) .extra <- NULL if (length(.w) > 0) { .v <- names(.ctl)[.w] .minfo(paste0("rxControl items read from fun: '", - paste(.v, collapse="', '"), "'")) + paste(.v, collapse="', '"), "'")) .extra <- setNames(lapply(.v, function(x) { get(x, envir=.meta) }), .v) + } do.call(rxSolve, c(list(NULL, params = NULL, events = NULL, inits = NULL), .lst, .extra, @@ -1220,6 +1221,24 @@ rxSolve.function <- function(object, params = NULL, events = NULL, inits = NULL, .rxControl$omega <- NULL } } + if (inherits(.rxControl$omega, "matrix")) { + .omega <- .rxControl$omega + .v <- vapply(dimnames(.omega)[[1]], + function(v) { + !(v %in% names(params)) + }, logical(1), USE.NAMES = FALSE) + if (length(.v) == 1L) { + if (!.v) .rxControl$omega <- NULL + } else { + .omega <- .omega[.v, .v] + if (all(dim(.omega) == c(0L, 0L))) { + .rxControl$omega <- NULL + } else { + .rxControl$omega <- .omega + } + } + + } if (inherits(.rxControl$omega, "matrix") && all(dim(.rxControl$omega) == c(0,0))) { .rxControl$omega <- NULL @@ -1235,6 +1254,24 @@ rxSolve.function <- function(object, params = NULL, events = NULL, inits = NULL, .rxControl$sigma <- NULL } } + if (inherits(.rxControl$sigma, "matrix")) { + .sigma <- .rxControl$sigma + .v <- vapply(dimnames(.sigma)[[1]], + function(v) { + !(v %in% names(params)) + }, logical(1), USE.NAMES = FALSE) + if (length(.v) == 1L) { + if (!.v) .rxControl$sigma <- NULL + } else { + .sigma <- .sigma[.v, .v, drop = FALSE] + if (all(dim(.sigma) == c(0L, 0L))) { + .rxControl$sigma <- NULL + } else { + .rxControl$sigma <- .sigma + } + } + + } if (inherits(.rxControl$sigma, "matrix") && all(dim(.rxControl$sigma) == c(0,0))) { .rxControl$sigma <- NULL diff --git a/R/ui-rename.R b/R/ui-rename.R index 8dca1ff1a..b376d144e 100644 --- a/R/ui-rename.R +++ b/R/ui-rename.R @@ -41,68 +41,149 @@ } list(line[[2]], line[[3]], .var.name, .var.name2) } - -#' Rename variables in the expression +#' Renames everything in one function #' -#' @param item Expression to recursively rename -#' @param new New name -#' @param old Old name -#' @return new expression with variable renamed -#' @author Matthew L. Fidler +#' @param item language item to process +#' @param lst list of renaming from .assertRenameErrorModelLine, ie list(new, old, newChar, oldChar) +#' @param isLhs is the expression the left handed side of the equation +#' @return expression renamed #' @noRd -.rxRenameRecursive <- function(item, new, old, isLhs=FALSE) { +#' @author Matthew L. Fidler +.rxRenameRecursiveAll <- function(item, lst, isLhs=FALSE) { if (is.atomic(item)) { return(item) } if (is.name(item)) { - if (identical(item, old)) { - return(new) - } else { - return(item) + .env <- new.env(parent=emptyenv()) + .env$new <- NULL + lapply(seq_along(lst), function(i) { + if (!is.null(.env$new)) return(NULL) + .curLst <- lst[[i]] + .old <- .curLst[[2]] + if (identical(item, .old)) { + .env$new <- .curLst[[1]] + } + return(NULL) + }) + if (!is.null(.env$new)) { + return(.env$new) } + 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, .rxRenameRecursive, new=new, old=old, isLhs=TRUE)) - if (is.call(.denom)) .denom <- as.call(lapply(.denom, .rxRenameRecursive, new=new, old=old, isLhs=TRUE)) + if (is.call(.num)) .num <- as.call(lapply(.num, .rxRenameRecursiveAll, lst=lst, isLhs=TRUE)) + if (is.call(.denom)) .denom <- as.call(lapply(.denom, .rxRenameRecursiveAll, lst=lst, isLhs=TRUE)) return(as.call(c(list(item[[1]]), .num, .denom))) - } else if (isLhs && identical(item[[1]], old) && length(item) == 2L && + } else if (isLhs && length(item) == 2L && is.numeric(item[[2]])) { - # handle x(0) = items - return(as.call(c(new, lapply(item[-1], .rxRenameRecursive, new=new, old=old, isLhs=isLhs)))) + .env <- new.env(parent=emptyenv()) + .env$new <- NULL + lapply(seq_along(lst), + function(i) { + if (!is.null(.env$new)) return(NULL) + .curLst <- lst[[i]] + .old <- .curLst[[2]] + if (identical(item[[1]], .old)) { + .env$new <- .curLst[[1]] + } + return(NULL) + }) + if (!is.null(.env$new)) { + # handle x(0) = items + return(as.call(c(.env$new, lapply(item[-1], .rxRenameRecursiveAll, lst=lst, isLhs=isLhs)))) + } } if (identical(item[[1]], quote(`=`)) || identical(item[[1]], quote(`<-`)) || identical(item[[1]], quote(`~`))) { - .elhs <- lapply(item[c(-1, -3)], .rxRenameRecursive, new=new, old=old, isLhs=TRUE) - .erhs <- lapply(item[c(-1, -2)], .rxRenameRecursive, new=new, old=old, isLhs=FALSE) + .elhs <- lapply(item[c(-1, -3)], .rxRenameRecursiveAll, lst=lst, isLhs=TRUE) + .erhs <- lapply(item[c(-1, -2)], .rxRenameRecursiveAll, lst=lst, isLhs=FALSE) return(as.call(c(item[[1]], .elhs, .erhs))) } else { - return(as.call(c(list(item[[1]]), lapply(item[-1], .rxRenameRecursive, new=new, old=old, isLhs=isLhs)))) + return(as.call(c(list(item[[1]]), lapply(item[-1], .rxRenameRecursiveAll, lst=lst, isLhs=isLhs)))) } } else { stop("unknown expression", call.=FALSE) } } -#' Rename one item in the rxui +#' Rename all items in matrix dimnames #' -#' @param rxui rxui for renaming -#' @param lst list with (new, old, newChr, oldChr) -#' @return Nothing, called for side effects +#' @param mat matrix +#' @param lst list for renaming +#' @return renamed matrix +#' @noRd #' @author Matthew L. Fidler +.rxRenameAllMat <- function(mat, lst) { + .d <- dimnames(mat)[[1]] + .d <- vapply(seq_along(.d), function(i) { + .env <- new.env(parent=emptyenv()) + .env$new <- NULL + .cur <- .d[i] + lapply(seq_along(lst), + function(j) { + if (!is.null(.env$new)) return(NULL) + .curLst <- lst[[j]] + .old <- .curLst[[4]] + if (.cur == .old) { + .env$new <- .curLst[[3]] + } + }) + if (!is.null(.env$new)) return(.env$new) + return(.cur) + }, character(1), USE.NAMES=FALSE) + dimnames(mat) <- list(.d, .d) + mat +} + +#' Rename all the items in the initialization data frame and model +#' +#' @param rxui the ui to process +#' @param lst the list of old and new expressions (like above) +#' @return Called for side effects #' @noRd -.rxRename1 <- function(rxui, lst) { +#' @author Matthew L. Fidler +.rxRenameAll <- function(rxui, lst) { + rxui <- rxUiDecompress(rxui) .iniDf <- rxui$iniDf - .w <- which(.iniDf$name == lst[[4]]) - if (length(.w) == 1) { - .iniDf$name[.w] <- lst[[3]] - rxui$iniDf <- .iniDf + .iniDf$name <- vapply(seq_along(.iniDf$name), + function(i) { + .env <- new.env(parent=emptyenv()) + .env$new <- NULL + .cur <- .iniDf$name[i] + lapply(seq_along(lst), + function(j) { + if (!is.null(.env$new)) return(NULL) + .curLst <- lst[[j]] + .old <- .curLst[[4]] + if (.cur == .old) { + .env$new <- .curLst[[3]] + } + }) + if (!is.null(.env$new)) return(.env$new) + return(.cur) + }, character(1), USE.NAMES=FALSE) + rxui$iniDf <- .iniDf + if (exists("sigma", rxui)) { + assign("sigma", .rxRenameAllMat(get("sigma", envir=rxui), lst), envir=rxui) + } + if (exists("thetaMat", rxui)) { + assign("thetaMat", .rxRenameAllMat(get("thetaMat", envir=rxui), lst), envir=rxui) + } + if (exists("meta", rxui)) { + .meta <- get("meta", rxui) + if (exists("sigma", .meta)) { + assign("sigma", .rxRenameAllMat(get("sigma", envir=.meta), lst), envir=.meta) + } + if (exists("thetaMat", .meta)) { + assign("thetaMat", .rxRenameAllMat(get("thetaMat", envir=.meta), lst), envir=.meta) + } } rxui$lstExpr <- lapply(seq_along(rxui$lstExpr), function(i) { - .rxRenameRecursive(rxui$lstExpr[[i]], new=lst[[1]], old=lst[[2]]) + .rxRenameRecursiveAll(rxui$lstExpr[[i]], lst=lst) }) } @@ -164,21 +245,29 @@ rxRename <- function(.data, ..., envir=parent.frame()) { #' @rdname rxRename #' @export .rxRename <- function(.data, ..., envir=parent.frame()) { + .inCompress <- FALSE + if (inherits(.data, "rxUi") && + inherits(.data, "raw")) { + .inCompress <- TRUE + } rxui <- assertRxUi(.data) + if (inherits(rxui, "raw")) { + rxui <- rxUiDecompress(rxui) + } .vars <- unique(c(rxui$mv0$state, rxui$mv0$params, rxui$mv0$lhs, rxui$predDf$var, rxui$predDf$cond, rxui$iniDf$name)) .modelLines <- .quoteCallInfoLines(match.call(expand.dots = TRUE)[-(1:2)], envir=envir) .lst <- lapply(seq_along(.modelLines), function(i) { .assertRenameErrorModelLine(.modelLines[[i]], .vars) }) rxui <- .copyUi(rxui) # copy ui so effects do not affect original - lapply(seq_along(.lst), function(i) { - .rxRename1(rxui, .lst[[i]]) - }) + .rxRenameAll(rxui, .lst) .ret <- rxui$fun() if (inherits(.data, "rxUi")) { - .x <- rxUiDecompress(.data) - .ret <- .newModelAdjust(.ret, .x, rename=TRUE) - .ret <- rxUiCompress(.ret) + ## .x <- rxUiDecompress(.data) + .ret <- .newModelAdjust(.ret, rxui, rename=TRUE) + if (.inCompress) { + .ret <- rxUiCompress(.ret) + } .cls <- setdiff(class(.data), class(.ret)) if (length(.cls) > 0) { class(.ret) <- c(.cls, class(.ret)) @@ -218,4 +307,4 @@ rxRename.default <- function(.data, ...) { .lst <- as.list(match.call()[-1]) .lst$.data <- .data do.call(.rxRename, c(.lst, list(envir=parent.frame(2)))) -} +} diff --git a/_pkgdown.yml b/_pkgdown.yml index 532ceaf32..da8be4eda 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -207,7 +207,9 @@ navbar: - text: "Getting started" - text: "rxode2 mini language syntax" href: articles/rxode2-syntax.html - - text: "rxode2 with piping ie %>%" + - text: "rxode2 model modification %>%" + href: articles/Modifying-Models.html + - text: "rxode2 pipeline ie %>%" href: articles/rxode2-pipeline.html - text: "Speeding up rxode2 ODE solving" href: articles/rxode2-speed.html diff --git a/man/reexports.Rd b/man/reexports.Rd index 365c14ab7..d6ab11a0c 100644 --- a/man/reexports.Rd +++ b/man/reexports.Rd @@ -62,6 +62,7 @@ \alias{rxSetIni0} \alias{.getLastIdLvl} \alias{.expandPars} +\alias{toTrialDuration} \title{Objects exported from other packages} \keyword{internal} \description{ @@ -75,7 +76,7 @@ below to see their documentation. \item{magrittr}{\code{\link[magrittr:pipe]{\%>\%}}} - \item{rxode2et}{\code{\link[rxode2et:dot-clearPipe]{.clearPipe}}, \code{\link[rxode2et:dot-collectWarnings]{.collectWarnings}}, \code{\link[rxode2et:dot-s3register]{.s3register}}, \code{\link[rxode2et]{add.dosing}}, \code{\link[rxode2et]{add.sampling}}, \code{\link[rxode2et]{as.et}}, \code{\link[rxode2et:rxEvid]{as.rxEvid}}, \code{\link[rxode2et:rxRateDur]{as.rxRateDur}}, \code{\link[rxode2et]{et}}, \code{\link[rxode2et]{etExpand}}, \code{\link[rxode2et]{etRbind}}, \code{\link[rxode2et]{etRep}}, \code{\link[rxode2et]{etSeq}}, \code{\link[rxode2et]{eventTable}}, \code{\link[rxode2et]{rxCbindStudyIndividual}}, \code{\link[rxode2et]{rxEtDispatchSolve}}, \code{\link[rxode2et]{rxEvid}}, \code{\link[rxode2et]{rxRateDur}}, \code{\link[rxode2et]{rxReq}}, \code{\link[rxode2et]{rxStack}}} + \item{rxode2et}{\code{\link[rxode2et:dot-clearPipe]{.clearPipe}}, \code{\link[rxode2et:dot-collectWarnings]{.collectWarnings}}, \code{\link[rxode2et:dot-s3register]{.s3register}}, \code{\link[rxode2et]{add.dosing}}, \code{\link[rxode2et]{add.sampling}}, \code{\link[rxode2et]{as.et}}, \code{\link[rxode2et:rxEvid]{as.rxEvid}}, \code{\link[rxode2et:rxRateDur]{as.rxRateDur}}, \code{\link[rxode2et]{et}}, \code{\link[rxode2et]{etExpand}}, \code{\link[rxode2et]{etRbind}}, \code{\link[rxode2et]{etRep}}, \code{\link[rxode2et]{etSeq}}, \code{\link[rxode2et]{eventTable}}, \code{\link[rxode2et]{rxCbindStudyIndividual}}, \code{\link[rxode2et]{rxEtDispatchSolve}}, \code{\link[rxode2et]{rxEvid}}, \code{\link[rxode2et]{rxRateDur}}, \code{\link[rxode2et]{rxReq}}, \code{\link[rxode2et]{rxStack}}, \code{\link[rxode2et]{toTrialDuration}}} \item{rxode2parse}{\code{\link[rxode2parse:dot-getLastIdLvl]{.getLastIdLvl}}, \code{\link[rxode2parse]{forderForceBase}}, \code{\link[rxode2parse]{rxDerived}}, \code{\link[rxode2parse]{rxSetIni0}}} diff --git a/man/rxode2.Rd b/man/rxode2.Rd index 43b14e371..b2d6cb892 100644 --- a/man/rxode2.Rd +++ b/man/rxode2.Rd @@ -325,7 +325,7 @@ compilation model. \if{html}{\out{
}}\preformatted{mod$simulationModel }\if{html}{\out{
}} -\if{html}{\out{
}}\preformatted{## rxode2 2.0.14.9000 model named rx_85848c9248e14e8cbf9a9b4a606b2010 model (ready). +\if{html}{\out{
}}\preformatted{## rxode2 2.0.14.9000 model named rx_22699ef487579a63baec999e5bee4f82 model (ready). ## x$state: depot, center ## x$stateExtra: cp ## x$params: tka, tcl, tv, add.sd, eta.ka, eta.cl, eta.v, rxerr.cp @@ -336,7 +336,7 @@ compilation model. mod$simulationIniModel }\if{html}{\out{
}} -\if{html}{\out{
}}\preformatted{## rxode2 2.0.14.9000 model named rx_40b4a6f44b577c0e14f674ccc10f618e model (ready). +\if{html}{\out{
}}\preformatted{## rxode2 2.0.14.9000 model named rx_85bfff7e67355fe5de08b65b7a424b1e 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-plot.R b/tests/testthat/test-plot.R index c716381b1..e13c15aa9 100644 --- a/tests/testthat/test-plot.R +++ b/tests/testthat/test-plot.R @@ -1,5 +1,5 @@ rxTest({ - + expect_plotlog <- function(o, timex, logx, logy, dat) { # Checking for the correct type for logx and logy is nontrivial, so j expect_named(o, c("timex", "logx", "logy", "dat")) @@ -77,7 +77,7 @@ rxTest({ test_that("plot() with invalid component throws an error", { skip_on_os("mac") - + pheno2 <- function() { ini({ tcl <- log(0.008) @@ -104,5 +104,5 @@ rxTest({ expect_warning(plot(sim, sim, "foo")) expect_error(plot(sim), NA) }) - + }) diff --git a/tests/testthat/test-ui-rename.R b/tests/testthat/test-ui-rename.R index 8f27e6040..d405f161a 100644 --- a/tests/testthat/test-ui-rename.R +++ b/tests/testthat/test-ui-rename.R @@ -1,5 +1,6 @@ rxTest({ test_that("rename for ui makes sense", { + ocmt <- function() { ini({ tka <- exp(0.45) @@ -147,8 +148,113 @@ rxTest({ expect_equal(tmp$lstExpr[[6]], quote(rate(dcmt) <- 1)) expect_equal(tmp$lstExpr[[7]], quote(dur(dcmt) <- 1)) expect_equal(tmp$lstExpr[[8]], quote(alag(dcmt) <- 1)) - - }) + + test_that("rename with sigma and thetaMat", { + + f <- function() { + description <- "BOLUS_2CPT_CLV1QV2 SINGLE DOSE FOCEI (120 Ind/2280 Obs) runODE032" + dfObs <- 2280 + dfSub <- 120 + sigma <- lotri({ + eps1 ~ 1 + }) + thetaMat <- lotri({ + theta1 + theta2 + theta3 + theta4 + RSV + eps1 + eta1 + + omega.2.1 + eta2 + omega.3.1 + omega.3.2 + eta3 + + omega.4.1 + omega.4.2 + omega.4.3 + eta4 ~ + c(0.000887681, + -0.00010551, 0.000871409, 0.000184416, -0.000106195, + 0.00299336, -0.000120234, -5.06663e-05, 0.000165252, + 0.00121347, 5.2783e-08, -1.56562e-05, 5.99331e-06, + -2.53991e-05, 9.94218e-06, 0, 0, 0, 0, 0, 0, -4.71273e-05, + 4.69667e-05, -3.64271e-05, 2.54796e-05, -8.16885e-06, + 0, 0.000169296, 0, 0, 0, 0, 0, 0, 0, 0, -7.37156e-05, + 2.56634e-05, -8.08349e-05, 1.37e-05, -4.36564e-06, + 0, 8.75181e-06, 0, 0.00015125, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6.63383e-05, + -8.19002e-05, 0.000548985, 0.000168356, 1.59122e-06, + 0, 3.48714e-05, 0, 4.31593e-07, 0, 0, 0.000959029, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, -9.49661e-06, 0.000110108, + -0.000306537, -9.12897e-05, 3.1877e-06, 0, 1.36628e-05, + 0, -1.95096e-05, 0, 0, -0.00012977, 0, 0, 0, 0.00051019) + }) + ini({ + theta1 <- 1.37034036528946 + label("log Cl") + theta2 <- 4.19814911033061 + label("log Vc") + theta3 <- 1.38003493562413 + label("log Q") + theta4 <- 3.87657341967489 + label("log Vp") + RSV <- c(0, 0.196446108190896, 1) + label("RSV") + eta1 ~ 0.101251418415006 + eta2 ~ 0.0993872449483344 + eta3 ~ 0.101302674763154 + eta4 ~ 0.0730497519364148 + }) + model({ + cmt(CENTRAL) + cmt(PERI) + cl <- exp(theta1 + eta1) + v <- exp(theta2 + eta2) + q <- exp(theta3 + eta3) + v2 <- exp(theta4 + eta4) + v1 <- v + scale1 <- v + k21 <- q/v2 + k12 <- q/v + d/dt(CENTRAL) <- k21 * PERI - k12 * CENTRAL - cl * CENTRAL/v1 + d/dt(PERI) <- -k21 * PERI + k12 * CENTRAL + f <- CENTRAL/scale1 + ipred <- f + rescv <- RSV + w <- ipred * rescv + ires <- DV - ipred + iwres <- ires/w + y <- ipred + w * eps1 + }) + } + + f <- f() + + f2 <- rxRename(f, eta.cl=eta1, eta.v=eta2, eta.q=eta3, eta.v2=eta4, eps=eps1) + + expect_equal(sort(intersect(dimnames(f2$omega)[[1]], + c("eta.cl", "eta.v", "eta.q", "eta.v2"))), + c("eta.cl", "eta.q", "eta.v", "eta.v2")) + + expect_equal(sort(intersect(dimnames(f2$thetaMat)[[1]], + c("eta.cl", "eta.v", "eta.q", "eta.v2", "eps"))), + c("eps", "eta.cl", "eta.q", "eta.v", "eta.v2")) + + expect_equal(dimnames(f2$sigma)[[1]], "eps") + + f <- rxUiDecompress(f) + f$sigma <- f$meta$sigma + f$thetaMat <- f$meta$thetaMat + rm("sigma", envir=f$meta) + rm("thetaMat", envir=f$meta) + f <- rxUiCompress(f) + + f3 <- rxRename(f, eta.cl=eta1, eta.v=eta2, eta.q=eta3, eta.v2=eta4, eps=eps1) + + expect_equal(sort(intersect(dimnames(f3$omega)[[1]], + c("eta.cl", "eta.v", "eta.q", "eta.v2"))), + c("eta.cl", "eta.q", "eta.v", "eta.v2")) + + expect_equal(sort(intersect(dimnames(f3$thetaMat)[[1]], + c("eta.cl", "eta.v", "eta.q", "eta.v2", "eps"))), + c("eps", "eta.cl", "eta.q", "eta.v", "eta.v2")) + + expect_equal(dimnames(f3$sigma)[[1]], "eps") + + }) + + }) diff --git a/vignettes/articles/Modifying-Models.Rmd b/vignettes/articles/Modifying-Models.Rmd new file mode 100644 index 000000000..ff8fd7f00 --- /dev/null +++ b/vignettes/articles/Modifying-Models.Rmd @@ -0,0 +1,851 @@ +--- +title: "Modifying Models" +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +```{r setup} +library(rxode2) +``` + +There are two fundamental operations that you may wish to do in +`rxode2`/`nlmixr2`. First you might want to modify your model (ie add +covariate effects, add between subject variability, etc). The second +thing you may wish to do is change initial estimates, change the +boundaries of the problem, fix/unfix the initial estimates, etc. + +## Modifying model + +There are a few tasks you might want to do with the overall model: + +- Change a line in the model + +- Add a line to the model + +- Rename parameters in the model + +- Combine different models + +- Create functions to add certain model features to the model + +We will go over the model piping and other functions that you can use +to modify models and even add your own functions that modify models. + +We will not cover any of the model modification functions in `nlmixr2lib` + +### Modifying a model line + +In my opinion, modifying lines in a model is likely the most common +task in modifying a model. We may wish to modify the model to have a +between subject variability or add a covariate effects. + +To begin of course you need a base model to modify. Let's start with a +very simple PK example, using the single-dose theophylline dataset +generously provided by Dr. Robert A. Upton of the University of +California, San Francisco: + +```{r iniModel} +one.compartment <- function() { + ini({ + tka <- 0.45; label("Ka") + tcl <- 1; label("Cl") + tv <- 3.45; label("V") + 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) + 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) + }) +} +``` + +If we believed we did not have enough absorption to support between +subject variability you can change the line to drop the between +subject by modifying a single line. To do this simply type the line +you want in the model piping expression: + +```{r removeEta} +mod <- one.compartment |> + model(ka <- exp(tka)) + +print(mod) +``` + +As expected, the line is modified. Also you can notice that the +initial estimate for the between subject variability is dropped since +it is no longer part of the model. + +If for some reason you wanted to add it back to the model you can +modify the model and add it back: + +```{r addEta} +mod2 <- mod |> + model(ka <- tka * exp(eta.ka)) + +print(mod2) +``` + +In this modification, the `eta.ka` is automatically assumed to be a +between subject variability parameter. Also since `eta.ka` is not +mu-referenced `rxode2` points this out. + +The automatic detection of `eta.ka` is because the name follows a +convention. Parameters starting or ending with the following names are +assumed to be between subject variability parameters: + +- eta (from NONMEM convention) +- ppv (per patient variability) +- psv (per subject variability) +- iiv (inter-individual variability) +- bsv (between subject variability) +- bpv (between patient variability) + +If this is not functioning correctly you can change it to a covariate +which you can add a type of initial estimate to later: + +```{r addCov} +mod2 <- mod |> + model(ka <- tka * exp(eta.ka) + WT * covWt, cov="eta.ka") + +print(mod2) +``` + +As seen above, the `eta.ka` in the above model is assumed to be a +data-input parameter or covariate instead of an estimated parameter. + +You can also note that `WT` is automatically recognized as a covariate +and `covWt` is automatically recognized as a covariate parameter. + +In general covariates and typical/population parameters are +automatically converted to estimated parameters based on the parameter +name starting with (or ending with): + +- tv (for typical value) +- t (also for typical value) +- pop (for population parameter) +- err (for error parameter) +- eff (for effect parameter) +- cov (for covariate parameters) + +This has a few notable exceptions for parameters like (`wt`, `sex` and +`crcl`) which are assumed to be covariates. + +If you don't want any automatic variable conversion, you can also use +`auto=FALSE`: + +```{r addWithoutAuto} +mod3 <- mod |> + model(ka <- tka * exp(eta.ka) + WT * covWt, auto=FALSE) + +print(mod3) +``` + +In this case all additional parameters (`eta.ka`, `WT`, and `covWt`) +are assumed to be parameters in the dataset. + +### Note on automatic detection of variables + +The automatic detection of variables is convenient for many models but +may not suit your style; If you do not like it you can always change +it by using `options()`: + +```{r option1} +options(rxode2.autoVarPiping=FALSE) +``` + +With this option disabled, all variables will be assumed to be +covariates and you will have to promote them to population parameters +with the `ini` block + +In the last example with this option enabled none of the variables +starting with `t` will be added to the model + +```{r option1ignore} +mod7 <- mod3 |> + model({ + emax <- exp(temax) + e0 <- exp(te0 + eta.e0) + ec50 <- exp(tec50) + kin <- exp(tkin) + kout <- exp(tkout) + }, append=FALSE) + +print(mod7) +``` + +Of course you could use it and then turn it back on: + +```{r option1resume} +options(rxode2.autoVarPiping=TRUE) +mod8 <- mod |> + model({ + emax <- exp(temax) + e0 <- exp(te0 + eta.e0) + ec50 <- exp(tec50) + kin <- exp(tkin) + kout <- exp(tkout) + }, append=FALSE) + +print(mod8) +``` + +Or you can use the +`withr::with_options(list(rxode2.autoVarPiping=FALSE), ...)` to turn +the option on temporarily. + +If you don't like the defaults for changing variables you could change +them as well with `rxSetPipingAuto()` + +For example if you only wanted variables starting or ending with `te` +you can change this with: + +```{r autoVarsChange} +rxSetPipingAuto(thetamodelVars = rex::rex("te")) + +mod9 <- mod |> + model({ + emax <- exp(temax) + e0 <- exp(te0 + eta.e0) + ec50 <- exp(tec50) + kin <- exp(tkin) + kout <- exp(tkout) + }, append=FALSE) + +print(mod9) +``` + +And as requested only the population parameters starting with `te` are +added to the `ini` block. + +If you want to reset the defaults you simply call `rxSetPipingAuto()` +without any arguments: + +```{r autoVarsReset} +rxSetPipingAuto() +mod10 <- mod |> + model({ + emax <- exp(temax) + e0 <- exp(te0 + eta.e0) + ec50 <- exp(tec50) + kin <- exp(tkin) + kout <- exp(tkout) + }, append=FALSE) +``` + +### Adding model lines + +There are three ways to insert lines in a `rxode2`/`nlmixr2` +model. You can add lines to the end of the model, after an expression +or to the beginning of the model all controlled by the `append` +option. + +Let's assume that there are two different assays that were run with +the same compound and you have noticed that they both have different +variability. + +You can modify the model above by adding some lines to the end of the +model by using `append=TRUE`: + +```{r appendLines} +mod4 <- mod |> + model({ + cp2 <- cp + cp2 ~ lnorm(lnorm.sd) + }, append=TRUE) + +print(mod4) +``` + +Perhaps instead you may want to add an indirect response model in +addition to the concentrations, you can choose where to add this: with +`append=lhsVar` where `lhsVar` is the left handed variable above where +you want to insert the new lines: + +```{r insertLines} +mod5 <- mod |> + model({ + PD <- 1-emax*cp/(ec50+cp) + ## + effect(0) <- e0 + kin <- e0*kout + d/dt(effect) <- kin*PD -kout*effect + }, append=d/dt(center)) +``` + +The last type of insertion you may wish to do is to add lines to the +beginning of the model by using `append=FALSE`: + +```{r prependLines} +mod6 <- mod5 |> + model({ + emax <- exp(temax) + e0 <- exp(te0 + eta.e0) + ec50 <- exp(tec50) + kin <- exp(tkin) + kout <- exp(tkout) + }, append=FALSE) + +print(mod6) +``` + +### Remove lines in the model + +The lines in a model can be removed in one of 2 ways either use +`-param` or `param <- NULL` in model piping: + +```{r} +mod7 <- mod6 |> + model(-emax) + +print(mod7) + +# Equivalently + +mod8 <- mod6 |> + model(emax <- NULL) + +print(mod8) +``` + + +### Rename parameters in a model + +You may want to rename parameters in a model, which is easy to do with +`rxRename()`. When `dplyr` is loaded you can even replace it with +`rename()`. The semantics are similar between the two functions, that +is you assigning `newVar=oldVar`. For example: + +```{r rename1} +mod11 <- mod10 |> + rxRename(drug1kout=kout, tv.drug1kout=tkout) + +print(mod11) +``` + +You can see every instance of the variable is named in the model is +renamed inside the `model` and `ini` block. + +For completeness you can see this with the `dplyr` verb (since it is a +S3 method): + +```{r rename2} +library(dplyr) +mod12 <- mod10 |> + rename(drug1kout=kout, tv.drug1kout=tkout) + +print(mod12) +``` + +### Combine different models + +You can also combine different models with `rxAppendModel()`. In +general they need variables in common to combine. This is because you +generally want the models to link between each other. In the below +example a pk and pd model this is done by renaming `cp` in the first +model to `ceff` in the second model: + +```{r append1} +ocmt <- function() { + ini({ + tka <- exp(0.45) # Ka + tcl <- exp(1) # Cl + tv <- exp(3.45); # log V + ## the label("Label name") works with all models + add.sd <- 0.7 + }) + model({ + ka <- tka + cl <- tcl + v <- tv + d/dt(depot) <- -ka * depot + d/dt(center) <- ka * depot - cl / v * center + cp <- center / v + cp ~ add(add.sd) + }) +} + +idr <- function() { + ini({ + tkin <- log(1) + tkout <- log(1) + tic50 <- log(10) + gamma <- fix(1) + idr.sd <- 1 + }) + model({ + kin <- exp(tkin) + kout <- exp(tkout) + ic50 <- exp(tic50) + d/dt(eff) <- kin - kout*(1-ceff^gamma/(ic50^gamma+ceff^gamma)) + eff ~ add(idr.sd) + }) +} + +rxAppendModel(ocmt %>% rxRename(ceff=cp), idr) +``` + +You will get an error if you try to combine models without variables +in common: + +```{r append2} +try(rxAppendModel(ocmt, idr)) +``` + +If you want to combine the models without respecting the having the +variables in common, you can use `common=FALSE`: + +```{r append3} +mod2 <- rxAppendModel(ocmt, idr, common=FALSE) |> + model(ceff=cp, append=ic50) # here we add the translation after the + # ic50 line to make it reasonable + +print(mod2) +``` + +### Creating more complex model modification functions + +These are pretty flexible, but you may want to do even more, so there +are some helper functions to help you create functions to do more. We +will discuss how to extract the model from the function and how to +update it. + +Lets start with a model: + +```{r modelExtractFn} +f <- function() { + ini({ + tka <- 0.45 + tcl <- 1 + tv <- 3.45 + eta.ka ~ 0.6 + eta.v ~ 0.1 + }) + model({ + ka <- exp(tka + eta.ka) + cl <- exp(tcl) + v <- exp(tv + eta.v) + d/dt(depot) <- -ka * depot + d/dt(center) <- ka * depot - cl/v * center + cp <- center/v + }) +} +``` + +Lets assume for a moment you want to remove an eta to `cl`. First you probably want to get all the model lines. You can do that with `modelExtract()`: + +```{r modelExtract} +totLines <- modelExtract(f, endpoint=NA) # endpoints should be included + +print(totLines) +``` + +Now you want to only worry about the `cl` line, you can subset here: + +```{r modelExtract2} +clLine <- modelExtract(f, cl, lines=TRUE) +line <- attr(clLine, "lines") +``` + +Now I wish to change the line to "cl \<- exp(tcl+eta.cl)" + +```{r modelExtract3} +totLines[line] <- "cl <- exp(tcl+eta.cl)" + +# For now lets remove the entire `ini` block (so you don't have to +# worry about syncing parameters). + +# + +ini(f) <- NULL + +model(f) <- totLines + +print(f) +``` + +Note that these functions do not modify the `ini({})` block. You may +have to modify the ini block first to make it a valid +`rxode2`/`nlmixr2` model. + +In this particular case, using model piping would be easier, but it +simply demonstrates two different way to extract model information and +a way to add information to the final model. + +These methods can be tricky because when using them you have to have +model that is parsed correctly. This means you have to make sure the +parameters and endpoints follow the correct rules + +## Modifying initial estimates + +The common items you want to do with initial estimates are: + +- Fix/Unfix a parameter + +- Change the initial condition values and bounds + +- Change the initial condition type + +- Change labels and transformations + +- Reorder parameters + +You may wish to create your own functions; we will discuss this too. + + +### Fixing or unfixing a parameter + +You can fix model estimates in two ways. The first is to fix the value +to whatever is in the model function, this is done by piping the model +parameter name (like `tka`) and setting it equal to `fix` (`%>% +ini(tka=fix)`). Below is a full example: + +```{r fixEstimate} +f <- function() { + ini({ + tka <- 0.45 + tcl <- 1 + tv <- 3.45 + add.sd <- c(0, 0.7) + eta.ka ~ 0.6 + eta.v ~ 0.1 + }) + model({ + ka <- exp(tka + eta.ka) + cl <- exp(tcl) + 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) + }) +} + +f2 <- f |> + ini(tka=fix) + +print(f2) +``` + +You can also fix the parameter to a different value if you wish; This +is very similar you can specify the value to fix inside of a `fix` +pseudo-function as follows: `%>% ini(tka=fix(0.1))`. A fully worked +example is below: + +```{r fixEstimate2} +f <- function() { + ini({ + tka <- 0.45 + tcl <- 1 + tv <- 3.45 + add.sd <- c(0, 0.7) + eta.ka ~ 0.6 + eta.v ~ 0.1 + }) + model({ + ka <- exp(tka + eta.ka) + cl <- exp(tcl) + 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) + }) +} + +f2 <- f |> + ini(tka=fix(0.1)) + +print(f2) +``` + +### Unfixing parameters + +You an unfix parameters very similarly to fixing. Instead of using the +`fix` keyword, you use the `unfix` keyword. So to unfix a parameter +(keeping its value) you would pipe the model using (`|> +ini(tka=unfix)`). Starting with the fixed model above a fully worked +example is: + +```{r unfix0} +print(f2) + +f3 <- f2 |> ini(tka=unfix) + +print(f3) +``` + +You can also unfix and change the initial estimate with +`ini(parameter=unfix(newEst))`: + +```{r unfix1} +print(f2) + +f3 <- f2 |> + ini(tka=unfix(10)) + +print(f3) +``` + + +### Changing the parameter values and possibly bounds + +#### Multiple parameter assignment + +You can also assign multiple parameters by providing them: + +- As a vector/list + +- As multiple lines in a piped `ini()` block + +- Using a covariance matrix + +In the case of a vector you can specify them and then pipe the model. + +For example: + +```{r pipeModel1} +ini1 <- c(tka=0.1, tcl=1, tv=3) + +f4 <- f |> ini(ini1) + +print(f4) + +# or equivalently + +ini1 <- list(tka=0.1, tcl=1, tv=3) + +f4a <- f |> ini(ini1) + +print(f4) +``` + +This can also be added with multiple lines or commas separating estimates: + +```{r pipeIni2} +# commas separating values: +f4 <- f |> ini(tka=0.1, tcl=1, tv=3) +print(f4) + +# multiple lines in {} +f4 <- f |> + ini({ + tka <- 0.2 + tcl <- 2 + tv <- 6 + }) + +print(f4) +``` + +You could also use a matrix to specify the covariance: + +```{r iniCov} +ome <- lotri(eta.ka + eta.v ~ c(0.6, + 0.01, 10.1)) + +f4 <- f |> ini(ome) + +print(f4) + +# or equavialtly use the lotri-type syntax for the omega: + +f4 <- f |> ini(eta.ka + eta.v ~ c(0.6, + 0.01, 0.2)) +print(f4) +``` + +#### Single parameter assignment + +The simplest way to change the initial parameter estimates is to +simply use `ini(parameter=newValue)`. You can also use `<-` or `~` to +change the value: + +A fully worked example showing all three types of initial value +modification is: + +```{r iniAssign1} +f3 <- f |> + ini(tka <- 0.1) + +f4 <- f |> + ini(tka=0.1) + +f5 <- f |> + ini(tka ~ 0.1) + +print(f5) +``` + +You can change the bounds like you do in the model specification by +using a numeric vector of `c(low, estimate)` or `c(low, estimate, +hi)`. Here is a worked example: + +```{r iniAssign2} +f3 <- f |> + ini(tka <- c(0, 0.1, 0.2)) + +print(f3) + + +f3 <- f |> + ini(tka <- c(0, 0.1)) + +print(f3) +``` + +Note by changing the parameters to their default values they might not +show up in the parameter printout: + +```{r iniAssign3} +f3 <- f |> + ini(tka <- c(0, 0.1, 0.2)) + +print(f3) + +# Now reassign +f4 <- f3 |> + ini(tka <- c(-Inf, 0.1, Inf)) + +print(f4) +``` + +#### Changing parameter types + +You can change the parameter type by two operators either by using +`-par` to convert the parameter to a covariate or `~par` to toggle +between population and individual parameters. + +Here is an example that does all 3: + +```{r parType} +# Switch population parameter to between subject variability parameter: +f4 <- f |> + ini( ~ tcl) + +print(f4) + +# Switch back to population parameter +f5 <- f4 |> + ini( ~ tcl) + +print(f5) + +# Change the variable to a covariate parameter (ie it doesn't have an +# initial estimate so remove it with the `-` operator): + +f6 <- f4 |> + ini(-tcl) + +print(f6) + +# You can change the covariate or remove the parameter estimate by +# `tcl <- NULL`: + +f6 <- f4 |> + ini(tcl <- NULL) + +print(f6) + +# to add it back as a between subject variability or population +# parameter you can pipe it as follows: + +f7 <- f6 |> + ini(tcl=4) + +print(f7) + + +f8 <- f6 |> + ini(tcl ~ 0.1) + +print(f8) +``` + +#### Changing parameter labels + +If you want to change/add a parameter label you assign the parameter +to `label("label to add")`. For example: + +```{r label0} +f4 <- f |> + ini(tka=label("Typical Ka (1/hr)")) + +print(f4) +``` + +You can also change the order while performing operations: + +```{r label1} +f5 <- f |> + ini(tka=label("Typical Ka (1/hr)"), append=tcl) + +print(f5) +``` + +If you want to remove the labels you can remove them with +`ini(par=label(NULL))`; For example: + +```{r label2} +f6 <- f |> + ini(tka=label(NULL)) + +print(f6) +``` + +#### Changing parameter transformations + +Back-transformations over-ride the back transformations in `nlmixr2` +models. They are very similar to the modification of the labels. + +Here you use `|> ini(tka=backTransform(exp))` to add an exponential +back-transformation for data: + +```{r parTrans0} +f7 <- f |> + ini(tka=backTransform(exp)) + +print(f7) +``` + +If you wish to remove them you can also do that with `|> +ini(tka=backTransform(NULL))`: + +```{r parTrans1} +f8 <- f |> + ini(tka=backTransform(NULL)) + +print(f8) +``` + +### More granular access of initial conditions + +Just like with `model()` you can modify the underlying data frame that +represents the `ini()` block. In this case I will simply change the +initial estimate of the first parameter (`tka`): + +```{r getIni} +f <- rxode2(f) + +ini <- f$iniDf + +print(ini) + +ini$est[1] <- 7 + +ini(f) <- ini + +print(f) +```