knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
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.
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
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:
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:
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:
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:
If this is not functioning correctly you can change it to a covariate which you can add a type of initial estimate to later:
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):
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
:
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.
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()
:
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
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:
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:
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:
rxSetPipingAuto() mod10 <- mod |> model({ emax <- exp(temax) e0 <- exp(te0 + eta.e0) ec50 <- exp(tec50) kin <- exp(tkin) kout <- exp(tkout) }, append=FALSE)
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
:
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:
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
:
mod6 <- mod5 |> model({ emax <- exp(temax) e0 <- exp(te0 + eta.e0) ec50 <- exp(tec50) kin <- exp(tkin) kout <- exp(tkout) }, append=FALSE) print(mod6)
The lines in a model can be removed in one of 2 ways either use
-param
or param <- NULL
in model piping:
mod7 <- mod6 |> model(-emax) print(mod7) # Equivalently mod8 <- mod6 |> model(emax <- NULL) print(mod8)
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:
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):
library(dplyr) mod12 <- mod10 |> rename(drug1kout=kout, tv.drug1kout=tkout) print(mod12)
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:
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:
try(rxAppendModel(ocmt, idr))
If you want to combine the models without respecting the having the
variables in common, you can use common=FALSE
:
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)
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:
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()
:
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:
clLine <- modelExtract(f, cl, lines=TRUE) line <- attr(clLine, "lines")
Now I wish to change the line to "cl \<- exp(tcl+eta.cl)"
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
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
Remove covariances between all parameters or a group of parameters
You may wish to create your own functions; we will discuss this too.
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:
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:
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)
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:
print(f2) f3 <- f2 |> ini(tka=unfix) print(f3)
You can also unfix and change the initial estimate with
ini(parameter=unfix(newEst))
:
print(f2) f3 <- f2 |> ini(tka=unfix(10)) print(f3)
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:
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:
# 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:
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)
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:
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:
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:
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)
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:
# 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)
If you want to change/add a parameter label you assign the parameter
to label("label to add")
. For example:
f4 <- f |> ini(tka=label("Typical Ka (1/hr)")) print(f4)
You can also change the order while performing operations:
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:
f6 <- f |> ini(tka=label(NULL)) print(f6)
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:
f7 <- f |> ini(tka=backTransform(exp)) print(f7)
If you wish to remove them you can also do that with |>
ini(tka=backTransform(NULL))
:
f8 <- f |> ini(tka=backTransform(NULL)) print(f8)
There are two approaches to removing covarinaces for between subject
variabilities: diag()
and -cov(var1, var2)
The diag()
removes either all covariance elements (with no
arguments) or any covariance elements included in the argument list:
fd <- function() { ini({ tka <- 0.45 tcl <- 1 tv <- 3.45 add.sd <- c(0, 0.7) eta.ka ~ 0.6 eta.v ~ c(0.01, 0.1) eta.cl ~ c(0.01, 0.01, 1) }) 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 you want to remove all covariances you can use diag() with no # arguments fd %>% ini(diag()) # if you want to remove only covariances with eta.ka you can use: fd %>% ini(diag(eta.ka)) # if you want to remove only the covariances with eta.ka and eta.v you can use: fd %>% ini(-cov(eta.ka, eta.v))
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
):
f <- rxode2(f) ini <- f$iniDf print(ini) ini$est[1] <- 7 ini(f) <- ini print(f)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.