knitr::opts_chunk$set( collapse = TRUE, comment = "#>") # Setup the moodels from PopED's # https://andrewhooker.github.io/PopED/articles/model_def_other_pkgs.html library(babelmixr2) library(PopED) library(mrgsolve) library(PKPDsim) library(rxode2) ff_analytic <- function(model_switch,xt,parameters,poped.db){ with(as.list(parameters),{ y=xt N = floor(xt/TAU)+1 f=(DOSE/V)*(KA/(KA - CL/V)) * (exp(-CL/V * (xt - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) - exp(-KA * (xt - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU))) return(list( f=f,poped.db=poped.db)) }) } code <- ' $PARAM CL=3.75, V=72.8, KA=0.25 $CMT DEPOT CENT $ODE dxdt_DEPOT = -KA*DEPOT; dxdt_CENT = KA*DEPOT - (CL/V)*CENT; $TABLE double CP = CENT/V; $CAPTURE CP ' moda <- mrgsolve::mcode("optim", code, atol=1e-8, rtol=1e-8,maxsteps=5000) ff_ode_mrg <- function(model_switch, xt, p, poped.db){ times_xt <- drop(xt) dose_times <- seq(from=0,to=max(times_xt),by=p[["TAU"]]) time <- sort(unique(c(0,times_xt,dose_times))) is.dose <- time %in% dose_times data <- tibble::tibble(ID = 1, time = time, amt = ifelse(is.dose,p[["DOSE"]], 0), cmt = ifelse(is.dose, 1, 0), evid = cmt, CL = p[["CL"]], V = p[["V"]], KA = p[["KA"]]) out <- mrgsolve::mrgsim_q(moda, data=data) f <- out$CP f <- f[match(times_xt,out$time)] return(list(f=matrix(f,ncol=1),poped.db=poped.db)) } sfg <- function(x,a,bpop,b,bocc){ parameters=c( KA=bpop[1]*exp(b[1]), CL=bpop[2]*exp(b[2]), V=bpop[3]*exp(b[3]), DOSE=a[1], TAU=a[2]) return( parameters ) } feps <- function(model_switch,xt,parameters,epsi,poped.db){ f <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db))[[1]] y = f*(1+epsi[,1])+epsi[,2] return(list(y=y,poped.db=poped.db)) } poped_db_analytic <- create.poped.database( ff_fun =ff_analytic, fg_fun =sfg, fError_fun=feps, bpop=c(KA=0.25,CL=3.75,V=72.8), d=c(KA=0.09,CL=0.25^2,V=0.09), sigma=c(prop=0.04,add=0.0025), m=2, groupsize=20, xt=c( 1,2,8,240,245), minxt=c(0,0,0,240,240), maxxt=c(10,10,10,248,248), bUseGrouped_xt=1, a=cbind(DOSE=c(20,40),TAU=c(24,24)), maxa=c(DOSE=200,TAU=24), mina=c(DOSE=0,TAU=24)) poped_db_ode_mrg <- create.poped.database(poped_db_analytic,ff_fun = ff_ode_mrg) pk1cmtoral <- PKPDsim::new_ode_model("pk_1cmt_oral") # take from library ff_ode_pkpdsim <- function(model_switch, xt, p, poped.db){ #Set up time points for the ODE times_xt <- drop(xt) dose_times <- seq(from=0,to=max(times_xt),by=p[["TAU"]]) times <- sort(unique(c(0,times_xt,dose_times))) # N = length(dose_times) regimen = PKPDsim::new_regimen(amt=p[["DOSE"]],n=N,interval=p[["TAU"]]) design <- PKPDsim::sim( ode = pk1cmtoral, parameters = c(CL=p[["CL"]],V=p[["V"]],KA=p[["KA"]]), regimen = regimen, only_obs = TRUE, t_obs = times, checks = FALSE, return_design = TRUE) tmp <- PKPDsim::sim_core(sim_object = design, ode = pk1cmtoral) f <- tmp$y m_tmp <- match(round(times_xt,digits = 6),tmp[,"t"]) if(any(is.na(m_tmp))){ stop("can't find time points in solution\n", "try changing the digits argument in the match function") } # f <- f[m_tmp] return(list(f = f, poped.db = poped.db)) } modrx <- rxode2::rxode2({ d/dt(DEPOT) = -KA*DEPOT; d/dt(CENT) = KA*DEPOT - (CL/V)*CENT; CP=CENT/V; }) ff_ode_rx <- function(model_switch, xt, p, poped.db){ times_xt <- drop(xt) et(0,amt=p[["DOSE"]], ii=p[["TAU"]], until=max(times_xt)) %>% et(times_xt) -> data out <- rxode2::rxSolve(modrx, p, data, atol=1e-8, rtol=1e-8,maxsteps=5000, returnType="data.frame") f <- out$CP[match(times_xt,out$time)] return(list(f=matrix(f,ncol=1),poped.db=poped.db)) } poped_db_ode_rx <- create.poped.database(poped_db_analytic,ff_fun = ff_ode_rx) poped_db_ode_pkpdsim <- create.poped.database(poped_db_analytic,ff_fun = ff_ode_pkpdsim)
babelmixr2
now introduces a new method that takes rxode2
/nlmixr2
models converts them to a PopED
database to help with optimal design.
As in the PopED vignette comparing ODE solvers (and their speeds), this section will:
take the model described and adapt it in two different rxode2
model functions, the solved and ode cases (this is done by the
nlmixr()
call which creates a PopED
database)
compare these examples to the pharmacometric solvers in the PopED
vignette (mrgsolve
, original rxode2
and PKPDsim
)
The first step of a design using babelmixr2 is to tell babelmixr2
about the design being optimized. This is a bit different than what
is done in PopED
directly. Below I am using the et()
function to
create the event table like a typical rxode2
simulation, but it is
used to specify the study design:
library(babelmixr2) library(PopED) e <- et(amt=1, ii=24, until=250) %>% et(time=c(1,2,8,240,245)) %>% as.data.frame() %>% dplyr::mutate(low=c(NA_real_, 0, 0, 0, 240, 240), high=c(NA_real_, 10, 10, 10, 248, 248)) print(e)
Here note that time
is the design times for the PopED
designs,
they can include dosing; only observations are considered the
time-points. They become the xt
parameter in the PopED
database
(excluding the doses).
We also build on the structure of the rxode2
event table with
simulations. In simulations the sampling windows cause random times
to be generated inside the sampling windows. For this reason, the
last line of code fixes the times to where we want to have the
multiple endpoint design.
Therefore, in this dataset the low
becomes minxt
and high
becomes maxxt
.
We chose these because they build on what is already know from
nlmixr2
and used in and do not require any extra coding.
Other things you may have to include in your PopED
model data frame are:
dvid
which gives the integer of the model endpoint measured (like
rxode2
but has to be an integer). This becomes model_switch
in
the PopED
dataset.
G_xt
which is the PopED
grouping variable; This will be put into
the PopED
database as G_xt
id
becomes an ID for a design (which you can use as a covariate to
pool different designs or different regimens for optimal design).
nlmixr2
/rxode2
ui functionOnce the design is setup, we need to specify a model. It is easy to
specify the model using the nlmixr2
/rxode2
function/ui below:
# model f <- function() { ini({ tKA <- 0.25 tCL <- 3.75 tV <- 72.8 Favail <- fix(0.9) eta.ka ~ 0.09 eta.cl ~ 0.25 ^ 2 eta.v ~ 0.09 prop.sd <- sqrt(0.04) add.sd <- sqrt(0.0025) }) model({ ka <- tKA * exp(eta.ka) v <- tV * exp(eta.v) cl <- tCL * exp(eta.cl) d/dt(depot) <- -ka * depot d/dt(central) <- ka * depot - cl / v * central cp <- central / v f(depot) <- DOSE * Favail cp ~ add(add.sd) + prop(prop.sd) }) } f <- f() # compile/check nlmixr2/rxode2 model # Create a PopED database for `nlmixr2`: poped_db_ode_babelmixr2 <- nlmixr(f, e, "poped", popedControl(a=list(c(DOSE=20), c(DOSE=40)), maxa=c(DOSE=200), mina=c(DOSE=0)))
Note when creating a PopED
database with a model and a design event
table, many of the PopED
database components are generated for you.
These are not things that are hidden, but things you can access
directly from the model or even from the compiled ui
. Much of the
other options for optimal design can be specified with the
popedControl()
function.
This can help you understand what babelmixr2
is doing, we will show
what is being added:
ff_fun
from babelmixr2This is the function that is run to generate the predictions:
# The ff_fun can be retrieved from the ui with f$popedFfFun f$popedFfFun
Some things to note in this function:
The model changes based on the number of time-points requested. In
this case it is 5
since there were 5
design points in the design above.
There are some babelmixr2
specific functions here:
babelmixr2::popedMultipleEndpointParam
, which indexes the time
input and model_input to make sure the input matches as requested
(in many PopED functions they use match.time
and this is a bit
similar).
.popedRxRunSetup
/.popedRxRunFullSetupMe
which runs the rxode2
setup including loading the data and model into memory (and is a
bit different depending on the number of time points you are
using)
.popedSolveIdME
/.popedSolveIdME2
which solves the rxode2 model and uses the
indexes to give the solve used in the model.
babelmixr2
By describing this, you can also see that there are 2 rxode2
models
generated for the PopED
database. You can see these inside of the
PopED database as well.
The first model uses model times to solve for arbitrary times based on design:
summary(poped_db_ode_babelmixr2$babelmixr2$modelMT)
You can also see the model used for solving scenarios with a number of time points greater than the design specification:
summary(poped_db_ode_babelmixr2$babelmixr2$modelF)
You can also see that the models are identical with the exception of
requesting modeled times. You can see the base/core rxode2
model
form the UI here:
f$popedRxmodelBase
fg_fun
babelmixr2
also generates PopED
s fg_fun
, which translates
covariates and parameters into the parameters required in the ff_fun
and used in solving the rxode2
model.
# You can see the PopED fg_fun from the model UI with # f$popedFgFun: f$popedFgFun
fError_fun
You can see the babelmixr2
generated error function as well with:
f$popedFErrorFun
One really important note to keep in mind is that PopED
works with
variances instead of standard deviations (which is a key difference
between nlmixr2
and PopED
).
This means that the exported model from babelmixr2
operates on
variances instead of standard deviations and care must be taken in
using these values to not misinterpret the two.
The export also tries to flag this to make it easier to remember.
PopED
f$popedBpop # PopED bpop f$popedNotfixedBpop # PopED notfixed_bpop f$popedD # PopED d f$popedNotfixedD # PopED notfixed_d f$popedCovd # PopED covd f$popedNotfixedCovd # PopED notfixed_covd f$popedSigma # PopED sigma (variance is exported, not SD) f$popedNotfixedSigma # PopED notfixed_sigma
The rest of the parameters are generated in conjunction with the
popedControl()
.
babelmixr2
You can also specify the models using the linCmt()
solutions as
below:
f2 <- function() { ini({ tV <- 72.8 tKA <- 0.25 tCL <- 3.75 Favail <- fix(0.9) eta.ka ~ 0.09 eta.cl ~ 0.25 ^ 2 eta.v ~ 0.09 prop.sd <- sqrt(0.04) add.sd <- fix(sqrt(5e-6)) }) model({ ka <- tKA * exp(eta.ka) v <- tV * exp(eta.v) cl <- tCL * exp(eta.cl) cp <- linCmt() f(depot) <- DOSE cp ~ add(add.sd) + prop(prop.sd) }) } poped_db_analytic_babelmixr2 <- nlmixr(f, e, popedControl(a=list(c(DOSE=20), c(DOSE=40)), maxa=c(DOSE=200), mina=c(DOSE=0)))
library(ggplot2) library(microbenchmark) compare <- microbenchmark( evaluate_design(poped_db_analytic), evaluate_design(poped_db_analytic_babelmixr2), evaluate_design(poped_db_ode_babelmixr2), evaluate_design(poped_db_ode_mrg), evaluate_design(poped_db_ode_pkpdsim), evaluate_design(poped_db_ode_rx), times = 100L) autoplot(compare) + theme_bw()
Note that the babelmixr2
ode solution is the fastest ode solver in this
comparison. Among other things, this is because the model is loaded into
memory and does not need to be setup each time. (As benchmarks, the
mrgsolve
, and PKPDsim
implementations on the PopED
's
website are included).
Also to me, the speed of all the tools are reasonable. In my opinion,
the benefit of the babelmixr2
interface to PopED
is the simplicity
of using nlmixr2
/ rxode2
functional models or fits directly in
PopED
without relying on conversions.
The interface is a bit different than the traditional PopED
interface,
and requires a design data-set as well as a popedControl()
to setup a
PopED
database to run all of the PopED
tasks. This is because
traditionally nlmixr2
takes a dataset, "estimation" method and
controls to change estimation method options.
babelmixr2
adopts the same paradigm of model, data, control to be
applied to PopED
. This should allow easy translation between the
systems. With easier translation, hopefully optimal design in clinical
trials will be easier to achieve.
babelmixr2
loads models into memory and needs to keep track of
which model is loaded. To help this you need to use
babel.poped.database
in place of create.poped.database
when
modifying babelmixr2 generated PopED
databases. If this isn't
done, there is a chance that the model loaded will not be the
expected loaded model and may either crash R or possibly give incorrect
results.
babelmixr2
translates all error components to variances instead of
the standard deviations in the nlmixr2
/rxode2
model
When there are covariances in the omega
specification, they will
be identified as D[#,#]
in the PopED
output. To see what these
numbers refer to it is helpful to see the name translations with
$popedD
.
Depending on your options, babelmixr2
may literally fix the model
components, which means indexes may be different than you
expect. The best way to get the correct index is use the
babelmixr2
function babelBpopIdx()
which is useful for using
PopED
babelmixr2
uses model times in creating PopED
databases;
therefore models with modeling times in them cannot be used in this translation
babelmixr2
does not yet support inter-occasion variability models.
The examples from PopED
have been converted to work with babelmixr2
and are available in the package and on GitHub
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.