knitr::opts_chunk$set( collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE, out.width = "100%" ) ## These options cache the models and the model simulations in R ## To run the actual models on your system, take the save options off. ## nlmixrVersion <- as.character(packageVersion("nlmixr")); ## options(nlmixr.save=TRUE, ## nlmixr.save.dir=file.path(system.file(package="nlmixr"), nlmixrVersion)); ## if (!dir.exists(getOption("nlmixr.save.dir"))) ## dir.create(getOption("nlmixr.save.dir"))
Joint PK/PD models, or PK/PD models where you fix certain components are common in pharmacometrics. A classic example, (provided by Tomoo Funaki and Nick Holford) is Warfarin.
library(nlmixr) library(ggplot2)
In this example, we have a transit-compartment (from depot to gut to central volume) PK model and an effect compartment for the PCA measurement.
Below is an illustrated example of a model that can be applied to the data:
pk.turnover.emax <- function() { ini({ tktr <- log(1) tka <- log(1) tcl <- log(0.1) tv <- log(10) ## eta.ktr ~ 1 eta.ka ~ 1 eta.cl ~ 2 eta.v ~ 1 prop.err <- 0.1 pkadd.err <- 0.1 ## temax <- logit(0.8) #temax <- 7.5 tec50 <- log(0.5) tkout <- log(0.05) te0 <- log(100) ## eta.emax ~ .5 eta.ec50 ~ .5 eta.kout ~ .5 eta.e0 ~ .5 ## pdadd.err <- 10 }) model({ ktr <- exp(tktr + eta.ktr) ka <- exp(tka + eta.ka) cl <- exp(tcl + eta.cl) v <- exp(tv + eta.v) ## #poplogit = log(temax/(1-temax)) emax=expit(temax+eta.emax) #logit=temax+eta.emax ec50 = exp(tec50 + eta.ec50) kout = exp(tkout + eta.kout) e0 = exp(te0 + eta.e0) ## DCP = center/v PD=1-emax*DCP/(ec50+DCP) ## effect(0) = e0 kin = e0*kout ## d/dt(depot) = -ktr * depot d/dt(gut) = ktr * depot -ka * gut d/dt(center) = ka * gut - cl / v * center d/dt(effect) = kin*PD -kout*effect ## cp = center / v cp ~ prop(prop.err) + add(pkadd.err) effect ~ add(pdadd.err) }) }
Notice there are two endpoints in the model cp
and effect
. Both
are modeled in nlmixr using the ~
"modeled by" specification.
To see more about how nlmixr will handle the multiple compartment model, it is quite informative to parse the model and print the information about that model. In this case an initial parsing would give:
ui <- nlmixr(pk.turnover.emax) ui
In the middle of the printout, it shows how the data must be formatted
(using the cmt
and dvid
data items) to allow nlmixr to model the
multiple endpoint appropriately.
Of course if you are interested you can directly access the
information in ui$multipleEndpoint
.
ui$multipleEndpoint
Notice that the cmt
and dvid
items can use the named variables
directly as either the cmt
or dvid
specification. This flexible
notation makes it so you do not have to rename your compartments to
run nlmixr model functions.
The other thing to note is that the cp
is specified by an ODE
compartment above the number of compartments defined in the RxODE
part of the nlmixr
model. This is because cp
is not a defined
compartment, but a related variable cp
.
The last thing to notice that the cmt
items are numbered cmt=5
for
cp
or cmt=4
for effect
even though they were specified in the
model first by cp
and cmt
. This ordering is because effect
is a
compartment in the RxODE
system. Of course cp
is related to the
compartment central
, and it may make more sense to pair cp
with
the central
compartment.
If this is something you want to have you can specify the compartment
to relate the effect to by the |
operator. In this case you would
change
cp ~ prop(prop.err) + add(pkadd.err)
to
cp ~ prop(prop.err) + add(pkadd.err) | central
With this change, the model could be updated to:
pk.turnover.emax2 <- function() { ini({ tktr <- log(1) tka <- log(1) tcl <- log(0.1) tv <- log(10) ## eta.ktr ~ 1 eta.ka ~ 1 eta.cl ~ 2 eta.v ~ 1 prop.err <- 0.1 pkadd.err <- 0.1 ## temax <- logit(0.8) tec50 <- log(0.5) tkout <- log(0.05) te0 <- log(100) ## eta.emax ~ .5 eta.ec50 ~ .5 eta.kout ~ .5 eta.e0 ~ .5 ## pdadd.err <- 10 }) model({ ktr <- exp(tktr + eta.ktr) ka <- exp(tka + eta.ka) cl <- exp(tcl + eta.cl) v <- exp(tv + eta.v) ## emax=expit(temax+eta.emax) ec50 = exp(tec50 + eta.ec50) kout = exp(tkout + eta.kout) e0 = exp(te0 + eta.e0) ## DCP = center/v PD=1-emax*DCP/(ec50+DCP) ## effect(0) = e0 kin = e0*kout ## d/dt(depot) = -ktr * depot d/dt(gut) = ktr * depot -ka * gut d/dt(center) = ka * gut - cl / v * center d/dt(effect) = kin*PD -kout*effect ## cp = center / v cp ~ prop(prop.err) + add(pkadd.err) | center effect ~ add(pdadd.err) }) } ui2 <- nlmixr(pk.turnover.emax2) ui2$multipleEndpoint
Notice in this case the cmt
variables are numbered sequentially and
the cp
variable matches the center
compartment.
When dvid
and cmt
are combined in the same dataset, the cmt
data
item is always used on the event information and the dvid
is used on
the observations. nlmixr
expects the cmt
data item to match the
dvid
item for observations OR to be either zero or one for the
dvid
to replace the cmt
information.
If you do not wish to use dvid
items to define multiple endpoints in
nlmixr, you can set the following option:
options(RxODE.combine.dvid=FALSE) ui2$multipleEndpoint
Then only cmt
items are used for the multiple endpoint models. Of
course you can turn it on or off for different models if you wish:
options(RxODE.combine.dvid=TRUE) ui2$multipleEndpoint
With this information, we can use the built-in warfarin dataset in nlmixr:
summary(warfarin)
Since dvid specifies pca
as the effect endpoint, you can update the
model to be more explicit making one last change:
cp ~ prop(prop.err) + add(pkadd.err) effect ~ add(pdadd.err)
to
cp ~ prop(prop.err) + add(pkadd.err) effect ~ add(pdadd.err) | pca
pk.turnover.emax3 <- function() { ini({ tktr <- log(1) tka <- log(1) tcl <- log(0.1) tv <- log(10) ## eta.ktr ~ 1 eta.ka ~ 1 eta.cl ~ 2 eta.v ~ 1 prop.err <- 0.1 pkadd.err <- 0.1 ## temax <- logit(0.8) tec50 <- log(0.5) tkout <- log(0.05) te0 <- log(100) ## eta.emax ~ .5 eta.ec50 ~ .5 eta.kout ~ .5 eta.e0 ~ .5 ## pdadd.err <- 10 }) model({ ktr <- exp(tktr + eta.ktr) ka <- exp(tka + eta.ka) cl <- exp(tcl + eta.cl) v <- exp(tv + eta.v) emax = expit(temax+eta.emax) ec50 = exp(tec50 + eta.ec50) kout = exp(tkout + eta.kout) e0 = exp(te0 + eta.e0) ## DCP = center/v PD=1-emax*DCP/(ec50+DCP) ## effect(0) = e0 kin = e0*kout ## d/dt(depot) = -ktr * depot d/dt(gut) = ktr * depot -ka * gut d/dt(center) = ka * gut - cl / v * center d/dt(effect) = kin*PD -kout*effect ## cp = center / v cp ~ prop(prop.err) + add(pkadd.err) effect ~ add(pdadd.err) | pca }) }
fit.TOS <- nlmixr(pk.turnover.emax3, warfarin, "saem", control=list(print=0), table=list(cwres=TRUE, npde=TRUE)) print(fit.TOS)
plot(fit.TOS); v1s <- nlmixr::vpc(fit.TOS, show=list(obs_dv=T), scales="free_y") + ylab("Warfarin Cp [mg/L] or PCA") + xlab("Time [h]") v2s <- nlmixr::vpc(fit.TOS, show=list(obs_dv=T), pred_corr = TRUE) + ylab("Prediction Corrected Warfarin Cp [mg/L] or PCA") + xlab("Time [h]") library(patchwork) v1s / v2s
## FOCEi fit/vpcs fit.TOF <- nlmixr(pk.turnover.emax3, warfarin, "focei", control=list(print=0), table=list(cwres=TRUE, npde=TRUE));
print(fit.TOF) plot(fit.TOF) v1f <- nlmixr::vpc(fit.TOF, show=list(obs_dv=T), scales="free_y") + ylab("Warfarin Cp [mg/L] or PCA") + xlab("Time [h]") v2f <- nlmixr::vpc(fit.TOF, show=list(obs_dv=T), pred_corr = TRUE) + ylab("Prediction Corrected Warfarin Cp [mg/L] or PCA") + xlab("Time [h]") v1f / v2f
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.