diff --git a/NAMESPACE b/NAMESPACE index 562719b8..a751a813 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -144,6 +144,7 @@ export(as.nlmixr) export(as.nlmixr2) export(as.nonmem2rx) export(babel.poped.database) +export(babelBpopIdx) export(bblDatToMonolix) export(bblDatToMrgsolve) export(bblDatToNonmem) diff --git a/NEWS.md b/NEWS.md index 9d0d5460..c554b364 100644 --- a/NEWS.md +++ b/NEWS.md @@ -13,6 +13,10 @@ `PopED` code between two different `PopED` databases, like in issue #131. +* Added a new function `babelBpopIdx(poped.db, "par")` which will get + the poped index for a model generated from `babelmixr2`, which is + useful when calculating the power (as in example 11). + # babelmixr2 0.1.4 * Added experimental `PopED` integration diff --git a/R/poped.R b/R/poped.R index 81d36423..a8c08c15 100644 --- a/R/poped.R +++ b/R/poped.R @@ -3042,7 +3042,13 @@ babel.poped.database <- function(popedInput, ..., optTime=NA) { if (is.environment(popedInput$babelmixr2)) { .babelmixr2 <- popedInput$babelmixr2 .db <- PopED::create.poped.database(popedInput=popedInput, ...) - .db$babelmixr2 <- .babelmixr2 + .b2 <- new.env(parent=emptyenv()) + for (v in ls(envir=.babelmixr2, all=TRUE)) { + assign(v, get(v, envir=.babelmixr2), envir=.b2) + } + .b2$modelNumber <- .poped$modelNumber + .poped$modelNumber <- .poped$modelNumber + 1 + .db$babelmixr2 <- .b2 if (is.logical(optTime) && length(optTime) == 1L && !is.na(optTime)) { popedMultipleEndpointResetTimeIndex() @@ -3054,3 +3060,68 @@ babel.poped.database <- function(popedInput, ..., optTime=NA) { call.=FALSE) } } + +#' Get the bpop_idx by variable name for a poped database created by `babelmixr2` +#' +#' This may work for other poped databases if the population parameters are named. +#' +#' @param popedInput The babaelmixr2 created database +#' @param var variable to query +#' @return index of the variable +#' @export +#' @author Matthew L. Fidler +#' +#' @examples +#' +#' if (requireNamespace("PopED", quietly=TRUE)) { +#' +#' f <- function() { +#' ini({ +#' tV <- 72.8 +#' tKa <- 0.25 +#' tCl <- 3.75 +#' tF <- fix(0.9) +#' pedCL <- 0.8 +#' +#' eta.v ~ 0.09 +#' eta.ka ~ 0.09 +#' eta.cl ~0.25^2 +#' +#' prop.sd <- fix(sqrt(0.04)) +#' add.sd <- fix(sqrt(5e-6)) +#' +#' }) +#' model({ +#' V<-tV*exp(eta.v) +#' KA<-tKa*exp(eta.ka) * (pedCL**isPediatric) # add covariate for pediatrics +#' CL<-tCl*exp(eta.cl) +#' Favail <- tF +#' +#' N <- floor(t/TAU)+1 +#' y <- (DOSE*Favail/V)*(KA/(KA - CL/V)) * +#' (exp(-CL/V * (t - (N - 1) * TAU)) * +#' (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - +#' exp(-KA * (t - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU))) +#' +#' y ~ prop(prop.sd) + add(add.sd) +#' }) +#' } +#' +#' e <- et(c( 1,8,10,240,245)) +#' +#' babel.db <- nlmixr2(f, e, "poped", +#' popedControl(m = 2, +#' groupsize=20, +#' bUseGrouped_xt=TRUE, +#' a=list(c(DOSE=20,TAU=24,isPediatric = 0), +#' c(DOSE=40, TAU=24,isPediatric = 0)))) +#' +#' babelBpopIdx(babel.db, "pedCL") +#' +#' } +babelBpopIdx <- function(popedInput, var) { + .w <- which(rownames(popedInput$parameters$bpop) == var) + if (length(.w) == 1L) return(.w) + stop("cannot find '", var, "' in the baelmixr2 poped model", + call.=FALSE) +} diff --git a/man/babelBpopIdx.Rd b/man/babelBpopIdx.Rd new file mode 100644 index 00000000..1210a249 --- /dev/null +++ b/man/babelBpopIdx.Rd @@ -0,0 +1,71 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/poped.R +\name{babelBpopIdx} +\alias{babelBpopIdx} +\title{Get the bpop_idx by variable name for a poped database created by \code{babelmixr2}} +\usage{ +babelBpopIdx(popedInput, var) +} +\arguments{ +\item{popedInput}{The babaelmixr2 created database} + +\item{var}{variable to query} +} +\value{ +index of the variable +} +\description{ +This may work for other poped databases if the population parameters are named. +} +\examples{ + +if (requireNamespace("PopED", quietly=TRUE)) { + +f <- function() { + ini({ + tV <- 72.8 + tKa <- 0.25 + tCl <- 3.75 + tF <- fix(0.9) + pedCL <- 0.8 + + eta.v ~ 0.09 + eta.ka ~ 0.09 + eta.cl ~0.25^2 + + prop.sd <- fix(sqrt(0.04)) + add.sd <- fix(sqrt(5e-6)) + + }) + model({ + V<-tV*exp(eta.v) + KA<-tKa*exp(eta.ka) * (pedCL**isPediatric) # add covariate for pediatrics + CL<-tCl*exp(eta.cl) + Favail <- tF + + N <- floor(t/TAU)+1 + y <- (DOSE*Favail/V)*(KA/(KA - CL/V)) * + (exp(-CL/V * (t - (N - 1) * TAU)) * + (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - + exp(-KA * (t - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU))) + + y ~ prop(prop.sd) + add(add.sd) + }) +} + +e <- et(c( 1,8,10,240,245)) + +babel.db <- nlmixr2(f, e, "poped", + popedControl(m = 2, + groupsize=20, + bUseGrouped_xt=TRUE, + a=list(c(DOSE=20,TAU=24,isPediatric = 0), + c(DOSE=40, TAU=24,isPediatric = 0)))) + +babelBpopIdx(babel.db, "pedCL") + +} +} +\author{ +Matthew L. Fidler +} diff --git a/src/poped.cpp b/src/poped.cpp index 9ca203fe..7bf1a7ca 100644 --- a/src/poped.cpp +++ b/src/poped.cpp @@ -218,6 +218,9 @@ Rcpp::NumericVector popedMultipleEndpointParam(Rcpp::NumericVector p, Rcpp::IntegerVector modelSwitch, int maxMT, bool optTime=true) { + if (optTime && globalTimeIndexer.isInitialized()) { + globalTimeIndexer.reset(); + } globalTimeIndexer.initialize(modelSwitch, times, optTime); Rcpp::NumericVector ret(p.size()-1+maxMT); std::fill(ret.begin(), ret.end(), globalTimeIndexer.getMaxTime()); @@ -467,7 +470,9 @@ void popedSolveFidMat(arma::mat &matMT, NumericVector &theta, int id, int nrow, setIndIdx(ind, j); kk = getIndIx(ind, j); curT = getTime(kk, ind); - isMT = getIndEvid(ind, kk) >= 10 && getIndEvid(ind, kk) <= 99; + int evid = getIndEvid(ind, kk); + isMT = evid >= 10 && evid <= 99; + // REprintf("curT: %f; evid: %d; isMT: %d; kk: %d\n", curT, evid, isMT, kk); if (isDose(getIndEvid(ind, kk))) { rxInner.calc_lhs(id, curT, getOpIndSolve(op, ind, j), lhs); continue; diff --git a/tests/testthat/test-poped.R b/tests/testthat/test-poped.R index 07d88e01..76de2a50 100644 --- a/tests/testthat/test-poped.R +++ b/tests/testthat/test-poped.R @@ -207,16 +207,16 @@ if (requireNamespace("PopED", quietly=TRUE)) { library(PopED) - expect_equal(list(ofv = 20.947993269923, - fim = lotri({ - tcl ~ 39.1196937156543 - tv ~ c(-0.386317381405573, 40.0037481563051) - sig_var_add.err ~ 800097.977314259 - }), - rse = c(tcl = 3.31152092974365, - tv = 30.9526397537534, - sig_var_add.err = 11.1796553130823)), - evaluate_design(db), tolerance = 1e-4) + expect_equal( + list(ofv = 20.9455595871912, + fim = lotri({ + tcl ~ 39.0235064308582 + tv ~ c(-0.385137678397696, 40.0037346438455) + sig_var_add.err ~ 800120.483866746 + }), + rse = c(tcl = 3.31559905061048, tv = 30.9526395967922, + sig_var_add.err = 11.1794980759702)), + evaluate_design(db), tolerance = 1e-4) ## v <- poped_optim(db, opt_xt=TRUE) @@ -379,27 +379,22 @@ if (requireNamespace("PopED", quietly=TRUE)) { tCl <- 3.75 tF <- fix(0.9) pedCL <- 0.8 - eta.v ~ 0.09 eta.ka ~ 0.09 eta.cl ~0.25^2 - prop.sd <- fix(sqrt(0.04)) add.sd <- fix(sqrt(5e-6)) - }) model({ V<-tV*exp(eta.v) - KA<-tKa*exp(eta.ka) * (pedCL**isPediatric) # add covariate for pediatrics - CL<-tCl*exp(eta.cl) + KA<-tKa*exp(eta.ka) + CL<-tCl*exp(eta.cl) * (pedCL^isPediatric) # add covariate for pediatrics Favail <- tF - N <- floor(t/TAU)+1 y <- (DOSE*Favail/V)*(KA/(KA - CL/V)) * (exp(-CL/V * (t - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - exp(-KA * (t - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU))) - y ~ prop(prop.sd) + add(add.sd) }) } @@ -436,6 +431,11 @@ if (requireNamespace("PopED", quietly=TRUE)) { vdiffr::expect_doppelganger("poped-ped-next", plot_model_prediction(babel.db.ped, model_num_points = 300)) + + expect_equal(babelBpopIdx(babel.db.ped, "pedCL"), 4L) + + expect_error(babelBpopIdx(babel.db, "matt")) + }) ## The tests run interactively runs OK