From 50cd5824ca8c43c5704c6ec841e472b0eac8cf1c Mon Sep 17 00:00:00 2001 From: "Matthew L. Fidler" Date: Mon, 18 Sep 2023 17:20:31 -0500 Subject: [PATCH 1/2] Add rbind method --- NAMESPACE | 1 + R/rbind.R | 107 ++++++++++++++++++++++++++++++++++++ R/rxsolve.R | 63 +++++++++++---------- man/meanProbs.Rd | 35 ++++++++++-- man/reexports.Rd | 1 - tests/testthat/test-rbind.R | 53 ++++++++++++++++++ 6 files changed, 224 insertions(+), 36 deletions(-) create mode 100644 R/rbind.R create mode 100644 tests/testthat/test-rbind.R diff --git a/NAMESPACE b/NAMESPACE index a56557c2d..d7b410224 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -86,6 +86,7 @@ S3method(print,rxSolveSimType) S3method(print,rxSymInvCholEnv) S3method(print,rxUi) S3method(print,rxode2) +S3method(rbind,rxSolve) S3method(rxChain2,EventTable) S3method(rxChain2,default) S3method(rxCompile,character) diff --git a/R/rbind.R b/R/rbind.R new file mode 100644 index 000000000..daea1c538 --- /dev/null +++ b/R/rbind.R @@ -0,0 +1,107 @@ +.rbind2rxSove <- function(v, v2) { + .tmp <- v$add.dosing # make sure info is calculated + .tmp <- v2$add.dosing + .cls <- class(v) + .env1 <- attr(class(v), ".rxode2.env") + .env2 <- attr(class(v2), ".rxode2.env") + if (length(.env1$.check.names) != length(.env2$.check.names) && + !all(.env1$.check.names == .env2$.check.names)) { + stop("cannot rbind these 2 rxSolve objects", call.=FALSE) + } + ## if (is.null(.env1$.params.single) || is.null(.env2$.params.single)) { + ## stop("cannot rbind single solve environments", call.=FALSE) + ## } + .cloneEnv <- new.env(parent=emptyenv()) + for (.v in ls(.env1, all=TRUE)) { + assign(.v, get(.v, envir=.env1), envir=.cloneEnv) + } + .nStud1 <- .cloneEnv$.args$nStud + .nStud2 <- .env2$.args$nStud + if (.nStud1 > 1 && .nStud2 > 1) { + .args <- .cloneEnv$.args + .args$nStud <- .nStud1 + .nStud2 + .cloneEnv$.args <- .args + .v <- v + class(.v) <- "data.frame" + .v2 <- v2 + .v2$sim.id <- .nStud1 + .v2$sim.id + class(.v2) <- "data.frame" + .v <- rbind(.v, .v2) + .cloneEnv$.args <- NULL + .cloneEnv$.nsub <- .cloneEnv$.nsub + .env2$.nsub + .cloneEnv$.et <- NULL + .cloneEnv$.args.params <- NULL + .cloneEnv$.args.inits <- NULL + .cloneEnv$.args.object <- rxModelVars(.cloneEnv$.args.object) + .cloneEnv$.args.par0 <- NULL + .cloneEnv$.args.params <- NULL + .cloneEnv$.check.nrow <- .cloneEnv$.check.nrow + .env2$.check.nrow + .cloneEnv$nobs <- .cloneEnv$nobs + .env2$nobs + .cloneEnv$.nsim <- .cloneEnv$.nsim + .env2$.nsim + .cloneEnv$.nsub <- NULL + .cloneEnv$.dadt.counter <- 0L + .cloneEnv$.init.dat <- setNames(rep(NA_real_, length(.cloneEnv$.init.dat)), + names(.cloneEnv$.init.dat)) + .cloneEnv$.jac.counter <- 0L + .cloneEnv$.nsub <- NA_integer_ + .cloneEnv$.par.pos <- NULL + .cloneEnv$.par.pos.ini <- NULL + .pd <- .env2$.params.dat + .pd$sim.id <- .nStud1 + .pd$sim.id + + .cloneEnv$.params.dat <- rbind(.cloneEnv$.params.dat, .pd) + .cloneEnv$.params.single <- NULL + .cloneEnv$.par.pos <- NULL + .cloneEnv$.par.pos.ini <- NULL + .cloneEnv$.et <- NULL + if (inherits(.v$id, "factor")) { + .cloneEnv$.idLevels <- levels(.v$id) + } + .cloneEnv$.jac.counter <- .cloneEnv$.slvr.counter <- 0L + .cloneEnv$.nsim <- .cloneEnv$.nsim + .env2$.nsim + .cloneEnv$.real.update <- FALSE + .cloneEnv$.sigma <- NULL + .fun <- function(...) { + stop("functions don't work on rbound rxSolve", call.=FALSE) + } + .cloneEnv$.replace.sampling <- .fun + .cloneEnv$add.dosing <- .fun + .cloneEnv$add.sampling <- .fun + .cloneEnv$clear.dosing <- .fun + .cloneEnv$clear.sampling <- .fun + .cloneEnv$get.dosing <- .fun + .cloneEnv$get.EventTable <- .fun + .cloneEnv$get.nobs <- .fun + .cloneEnv$get.obs.rec <- .fun + .cloneEnv$get.sampling <- .fun + .cloneEnv$get.units <- .fun + .cloneEnv$import <- .fun + .cloneEnv$counts.EventTable <- NULL + .cloneEnv$get.units <- NULL + .cloneEnv$units <- c(dosing="NA", time="NA") + .cloneEnv$dll <- NULL + attr(.cls, ".rxode2.env") <- .cloneEnv + class(.v) <- .cls + return(.v) + } +} + +#' @export +rbind.rxSolve <- function(..., deparse.level = 1) { + .lst <- list(...) + if (length(.lst) >= 2) { + .ret <- .rbind2rxSove(.lst[[1]], .lst[[2]]) + if (length(.lst) == 2) { + return(.ret) + } + return(do.call(rbind, + c(list(.ret), + lapply(seq_along(.lst)[-(1:2)], + function(i){ + .lst[[.]] + })))) + } + if (length(.lst) == 1) return(.lst[[1]]) + stop("called rbind.rxSolve() with no arguments", + call.=FALSE) +} diff --git a/R/rxsolve.R b/R/rxsolve.R index 34a52723e..8a48b2779 100644 --- a/R/rxsolve.R +++ b/R/rxsolve.R @@ -1,4 +1,4 @@ - #' Options, Solving & Simulation of an ODE/solved system +#' Options, Solving & Simulation of an ODE/solved system #' #' This uses rxode2 family of objects, file, or model specification to #' solve a ODE system. There are many options for a solved rxode2 @@ -1212,8 +1212,8 @@ rxSolve.function <- function(object, params = NULL, events = NULL, inits = NULL, .rx <- object$simulationModel } list(list(object=.rx, params = params, events = events, inits = inits), - .rxControl, - list(theta = theta, eta = eta)) + .rxControl, + list(theta = theta, eta = eta)) } #' @rdname rxSolve @@ -1238,7 +1238,7 @@ rxSolve.rxUi <- function(object, params = NULL, events = NULL, inits = NULL, ... if (is.null(.lst$omega) && is.null(.lst$sigma)) { .pred <- TRUE if (!.hasIpred && any(rxModelVars(.lst[[1]])$lhs == "ipredSim")) { - .lst$drop <- c(.lst$drop, "ipredSim") + .lst$drop <- c(.lst$drop, "ipredSim") } } .ret <- do.call("rxSolve.default", .lst) @@ -1301,6 +1301,9 @@ rxSolve.nlmixr2FitCore <- rxSolve.nlmixr2FitData #' @export rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL, ..., theta = NULL, eta = NULL) { + if (inherits(object, "rxSolve") && is.null(.v$env$.ev)) { + stop("cannot solve from more than one rxSolve object bound together", call.=FALSE) + } on.exit({ .clearPipe() .asFunctionEnv$rx <- NULL @@ -1310,13 +1313,13 @@ rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL, if (rxIs(object, "rxEt")) { if (!is.null(events)) { stop("events can be pipeline or solving arguments not both", - call. = FALSE - ) + call. = FALSE + ) } if (is.null(rxode2et::.pipeRx(NA))) { stop("need an rxode2 compiled model as the start of the pipeline", - call. = FALSE - ) + call. = FALSE + ) } else { events <- object object <- rxode2et::.pipeRx(NA) @@ -1328,16 +1331,16 @@ rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL, } if (is.null(rxode2et::.pipeRx(NA))) { stop("need an rxode2 compiled model as the start of the pipeline", - call. = FALSE - ) + call. = FALSE + ) } else { .rxParams <- object object <- rxode2et::.pipeRx(NA) } if (is.null(rxode2et::.pipeEvents(NA))) { stop("need an rxode2 events as a part of the pipeline", - call. = FALSE - ) + call. = FALSE + ) } else { events <- rxode2et::.pipeEvents(NA) rxode2et::.pipeEvents(NULL) @@ -1347,29 +1350,29 @@ rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL, events <- rxode2et::.pipeEvents(NA) } else if (!is.null(rxode2et::.pipeEvents(NA)) && !is.null(events)) { stop("'events' in pipeline AND in solving arguments, please provide just one", - call. = FALSE - ) + call. = FALSE + ) } else if (!is.null(rxode2et::.pipeEvents(NA)) && !is.null(params) && - rxIs(params, "event.data.frame")) { + rxIs(params, "event.data.frame")) { stop("'events' in pipeline AND in solving arguments, please provide just one", - call. = FALSE - ) + call. = FALSE + ) } if (!is.null(rxode2et::.pipeParams(NA)) && is.null(params)) { params <- rxode2et::.pipeParams(NA) } else if (!is.null(rxode2et::.pipeParams(NA)) && !is.null(params)) { stop("'params' in pipeline AND in solving arguments, please provide just one", - call. = FALSE - ) + call. = FALSE + ) } if (!is.null(rxode2et::.pipeInits(NA)) && is.null(inits)) { inits <- rxode2et::.pipeInits(NA) } else if (!is.null(rxode2et::.pipeInits(NA)) && !is.null(inits)) { stop("'inits' in pipeline AND in solving arguments, please provide just one", - call. = FALSE - ) + call. = FALSE + ) } if (.applyParams) { @@ -1380,13 +1383,13 @@ rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL, .xtra <- list(...) if (any(duplicated(names(.xtra)))) { stop("duplicate arguments do not make sense", - call. = FALSE - ) + call. = FALSE + ) } if (any(names(.xtra) == "covs")) { stop("covariates can no longer be specified by 'covs'\n include them in the event dataset\n\nindividual covariates: Can be specified by a 'iCov' dataset\n each each individual covariate has a value\n\ntime varying covariates: modify input event data-frame or\n 'eventTable' to include covariates(https://tinyurl.com/y52wfc2y)\n\nEach approach needs the covariates named to match the variable in the model", - call. = FALSE - ) + call. = FALSE + ) } .nms <- names(as.list(match.call())[-1]) .lst <- list(...) @@ -1400,9 +1403,9 @@ rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL, .both <- intersect(.mv$params, .ctl$keep) if (length(.both) > 0) { .keep <- .ctl$keep[!(.ctl$keep %in% .both)] - if (length(.keep) == 0L) { - .keep <- NULL - } + if (length(.keep) == 0L) { + .keep <- NULL + } .w <- which(names(.ctl) == "keep") .ctl[[.w]] <- .keep .ctl <- do.call(rxControl, @@ -1536,8 +1539,8 @@ rxSolve.default <- function(object, params = NULL, events = NULL, inits = NULL, events <- c(.theta, .eta) } else { stop("cannot specify 'params' and 'theta'/'eta' at the same time", - call. = FALSE - ) + call. = FALSE + ) } } if (!is.null(.ctl$iCov)) { diff --git a/man/meanProbs.Rd b/man/meanProbs.Rd index 0dcec476a..0a72bbb16 100644 --- a/man/meanProbs.Rd +++ b/man/meanProbs.Rd @@ -13,7 +13,7 @@ meanProbs(x, ...) na.rm = FALSE, names = TRUE, useT = TRUE, - useProbs = FALSE + onlyProbs = TRUE ) } \arguments{ @@ -32,8 +32,15 @@ vectors unless ‘na.rm’ is ‘TRUE’.} the confidence-based estimates. If false use the normal distribution to calculate the confidence based estimates.} -\item{useProbs}{logical; if true, only return the probability -based confidence interval estimates, otherwise return} +\item{onlyProbs}{logical; if true, only return the probability based +confidence interval estimates, otherwise return} +} +\value{ +By default the return has the probabilities as names (if +named) with the points where the expected distribution are +located given the sampling mean and standard deviation. If +\code{onlyProbs=FALSE} then it would prepend mean, variance, standard +deviation, minimum, maximum and number of non-NA observations. } \description{ The generic function \code{meanProbs} produces expected quantiles under @@ -53,11 +60,29 @@ The smallest observation corresponds to a probability of 0 and the largest to a probability of 1 and the mean corresponds to 0.5. The mean and standard deviation of the sample is calculated based -on Welford's method for a single pass +on Welford's method for a single pass. + +This is meant to perform in the same way as \code{quantile()} so it can +be a drop in replacement for code using \code{quantile()} but using +distributional assumptions. } \examples{ -meanProbs(x <- rnorm(1001)) +quantile(x<- rnorm(1001)) +meanProbs(x) + +# Can get some extra statistics if you request onlyProbs=FALSE +meanProbs(x, onlyProbs=FALSE) + +x[2] <- NA_real_ + +meanProbs(x, onlyProbs=FALSE) + +quantile(x<- rnorm(42)) + +meanProbs(x) + +meanProbs(x, useT=FALSE) } \author{ diff --git a/man/reexports.Rd b/man/reexports.Rd index 365c14ab7..f1310e9e9 100644 --- a/man/reexports.Rd +++ b/man/reexports.Rd @@ -82,4 +82,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/tests/testthat/test-rbind.R b/tests/testthat/test-rbind.R new file mode 100644 index 000000000..292587c32 --- /dev/null +++ b/tests/testthat/test-rbind.R @@ -0,0 +1,53 @@ +test_that("rbind", { + + mod2 <- function() { + ini({ + KA <- 2.94E-01 + TCL <- 1.86E+01 + V2 <- 4.02E+01 + Q <- 1.05E+01 + V3 <- 2.97E+02 + Kin <- 1 + Kout <- 1 + EC50 <- 200 + eta.Cl ~ 0.2 + eff.sd <- sqrt(0.05) + c2.sd <- sqrt(0.05) + }) + model({ + CL <- TCL * exp(eta.Cl) + d/dt(depot) <- -KA * depot + d/dt(centr) <- KA * depot - CL * (centr / V2) - Q * (centr / V2) + Q * (peri / V3) + d/dt(peri) <- Q * (centr / V2) - Q * (peri / V3) + C2 <- centr/V2 + d/dt(eff) <- Kin - Kout * (1 - C2 / (EC50 + C2)) * eff + eff(0) <- 1000 + eff ~ add(eff.sd) + C2 ~ prop(c2.sd) + }) + } + + ev <- eventTable() %>% + add.dosing(dose = 10000, nbr.doses = 10, dosing.interval = 12, dosing.to = 2) %>% + add.dosing(dose = 20000, nbr.doses = 5, start.time = 120, dosing.interval = 24, dosing.to = 2) %>% + add.sampling(0:240) %>% + et(id=1:2) %>% + as.data.frame() + + ev$dvid <- ifelse(ev$evid == 0, 1, NA_real_) + + ev2 <- ev[which(ev$dvid == 1), ] + ev2$dvid <- 2 + + ev <- rbind(ev, ev2) %>% + dplyr::arrange(id, time) + + rownames(ev) <- NULL + + pk1 <- + suppressWarnings(rxSolve(mod2, nStud = 4, ev, cores = 2)) + + pk2 <- + suppressWarnings(rxSolve(mod2, nStud = 4, ev, cores = 2)) + +}) From c4b9ac39e6011ef9bbbaa2f0c88326895761d461 Mon Sep 17 00:00:00 2001 From: "Matthew L. Fidler" Date: Mon, 18 Sep 2023 17:31:53 -0500 Subject: [PATCH 2/2] rbind fixes and tests --- R/rbind.R | 5 ++--- man/reexports.Rd | 1 + tests/testthat/test-rbind.R | 7 +++++++ 3 files changed, 10 insertions(+), 3 deletions(-) diff --git a/R/rbind.R b/R/rbind.R index daea1c538..493aaefa8 100644 --- a/R/rbind.R +++ b/R/rbind.R @@ -27,7 +27,6 @@ .v2$sim.id <- .nStud1 + .v2$sim.id class(.v2) <- "data.frame" .v <- rbind(.v, .v2) - .cloneEnv$.args <- NULL .cloneEnv$.nsub <- .cloneEnv$.nsub + .env2$.nsub .cloneEnv$.et <- NULL .cloneEnv$.args.params <- NULL @@ -94,11 +93,11 @@ rbind.rxSolve <- function(..., deparse.level = 1) { if (length(.lst) == 2) { return(.ret) } - return(do.call(rbind, + return(do.call(rbind.rxSolve, c(list(.ret), lapply(seq_along(.lst)[-(1:2)], function(i){ - .lst[[.]] + .lst[[i]] })))) } if (length(.lst) == 1) return(.lst[[1]]) diff --git a/man/reexports.Rd b/man/reexports.Rd index f1310e9e9..365c14ab7 100644 --- a/man/reexports.Rd +++ b/man/reexports.Rd @@ -82,3 +82,4 @@ 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/tests/testthat/test-rbind.R b/tests/testthat/test-rbind.R index 292587c32..8dcc8b615 100644 --- a/tests/testthat/test-rbind.R +++ b/tests/testthat/test-rbind.R @@ -50,4 +50,11 @@ test_that("rbind", { pk2 <- suppressWarnings(rxSolve(mod2, nStud = 4, ev, cores = 2)) + b1 <- rbind.rxSolve(pk1, pk2) + + b2 <- rbind.rxSolve(pk1, pk2, pk1) + + expect_true(inherits(b1, "rxSolve")) + expect_true(inherits(b2, "rxSolve")) + expect_error(rbind.rxSolve()) })