assemblerr is an R package helping you build pharmacometric models quickly and without hassle. The package allows you to assemble a model from a set of predefined building blocks and then convert it to NONMEM code. This document will show you the most important steps when working with assemblerr and is an excellent place to start if you are new to the package. We will use a few examples to discuss how to :
Let's start simple and assume you want to build a dose-response model using an Emax function. The following code would accomplish that:
m <- model() + prm_log_normal("emax", median = 10, var_log = 0.09) + prm_log_normal("ed50", median = 50, var_log = 0.09) + algebraic(effect~emax*dose/(ed50 + dose)) + obs_additive(~effect, var_add = 1)
The first line is the so-called constructor, something like the foundation of what you are trying to build. In this example, we are building a general pharmacometric model. Later, we will discuss when and why you would choose something else. All other lines specify building blocks and combine them using the
+ operator. What each of the building blocks defines is hopefully clear from its name, but let's go through them anyway to be sure:
Rformula. In these formulae, the
~should be understood like a
=. You can use formulae in assemblerr whenever you want to define a mathematical relationship.
Let's assume you want to create a model with a proportional error model instead and change the initial estimate for it, you would simply need to exchange the last line to something like
obs_proportional(~effect, var_prop = 0.05). Similarly, if you would want to switch one of the parameters shown above, you could also change it. The use of these exchangeable building blocks is the core principle in assemblerr.
When you are done defining your model you can check whether it contains any issues using the
check() function, e.g.,
In this case, the model definition used the variable 'dose' without it being declared. You can tell assemblerr that dose is part of the dataset by modifying the model as follows:
m <- m + input_variable("dose")
Now, the check function does not find any issues:
You can also replace an already existing building block by adding a new one with the same name. For example, the following code updates the parameter model for the Emax parameter to one without variability:
m <- m + prm_no_var("emax", value = 10)
The warning message indicates that an existing building block has been replace. Now, there won't be any random effect associated with Emax when you render the modified model:
After having defined your model, you probably want to turn it into something useful, in assemblerr this process is called rendering. Since we stored our assembled model above in the variable
m we can simply pass this variable to the
render function, which will write the resulting NONMEM code to the console, like so
While this is great to quickly check the resulting code, we generally want a model file. We can create it using the
filename= argument of the
render() function, i.e.,
setwd(tempdir()) render(m, filename = "run1.mod")
render function also has an options argument that allows you to modify some aspects of the rendering process. For example, the following code automatically adds the mu-referenced code for all parameters.
render( m, options = assemblerr_options( prm.use_mu_referencing = TRUE ) )
"All models are wrong" AND useless if we cannot do anything with them. So far, you only saw the description of a model. The
tasks= argument of the
render() function allows you to specify what the model should "do".
By default assemblerr adds an estimation task using the FOCE algorithm but you can change that, for example using:
render(m, tasks = tsk_estimation(algorithm = "imp", se = TRUE))
This example also requested the standard errors for the estimated parameters, which resulted in the addition of the covariance step in the NONMEM code.
You can combine multiple tasks using the
+ operator. For example, you might want to estimate parameters and output the parameter values for all parameters in the model. The following code accomplishes this:
render(m, tasks = tsk_estimation() + tsk_output(filename = "sdtab", variables = vars_prms()))
In the previous example, you saw the variable selection helper
vars_prms() in action. Variable selection helpers allow you to specify what variables should be used in a table. You can read more about them in the documentation for the
Many pharmacometric models are defined using ordinary differential equations. In assemblerr you can specify ordinary differential equations via compartments and flows (while this does not cover all types of ordinary differential equations, it will hopefully get you pretty far). To show compartment and flows in action, we will extend the dose-response model from the beginning to a concentration-response model. The following code chunk shows the necessary changes:
m2 <- model() + prm_log_normal("v", median = 100, var_log = 0.09) + prm_log_normal("cl", median = 10, var_log = 0.09) + prm_log_normal("emax", median = 10, var_log = 0.09) + prm_log_normal("ec50", median = 1, var_log = 0.09) + compartment("central", volume = "v") + flow(from = "central", definition = ~cl*C) + obs_additive(effect~emax*C["central"]/(ec50 + C["central"]), var_add = 1)
In this model definition, we used two additional parameters and, more importantly, a compartment as well as a flow from that compartment. You should note, how you can use the variable
C in the flow definition to refer to the concentration in the "from" compartment. Similarly, you can also use
C["central"] to refer to the concentration in the central compartment when defining the observational model. The
A variable provides a reference to the amount in the compartment.
Pay attention to how assemblerr replaces the
C variables with the corresponding compartment reference when we render the model like before:
Did you also notice what happened to the differential equation? It disappeared! assemblerr recognized that this simple one-compartment model could be written as a much quicker analytical solution using ADVAN1 and did the necessary replacements.
For PK models, assemblerr supports building blocks with a higher level of abstraction. The following snippet assembles a one-compartment model with linear elimination:
pkm <- pk_model() + pk_absorption_fo() + pk_distribution_1cmp() + pk_elimination_linear() + obs_additive(conc~C["central"])
You can see how
pk_model() instead of
model() was used as a constructor. The resulting model objects accepts building blocks such as:
as well as the additive observation model (line 5) that you saw before.
Despite its slightly different nature, you can turn the model into NONMEM code using the
render() function. Also, tasks and rendering options work like in the previous examples. You could therefore do the following to create NONMEM code with mu-referencing, importance sampling estimation, and a table containing all model parameters:
render( model = pkm, tasks = tsk_estimation("imp") + tsk_output("sdtab", variables = vars_prms()), options = assemblerr_options(prm.use_mu_referencing = TRUE) )
PK model building blocks come with a default parameter model for their parameters. You can change the default by providing the parameter model as an argument to the building blocks. For example, if you want to have no variability on the volume of distribution, you could use the following code (note also how VC was renamed to V1):
pkm <- pk_model() + pk_absorption_fo() + pk_distribution_1cmp(prm_vc = prm_no_var("v1", value = 100)) + pk_elimination_linear() + obs_additive(conc~C["central"]) render(pkm)
Alternatively, you can overwrite the default parameter model by adding a another one with the same name, i.e., :
pkm <- pk_model() + pk_absorption_fo() + pk_distribution_1cmp() + pk_elimination_linear() + prm_no_var("vc", value = 100) + obs_additive(conc~C["central"]) render(pkm)
There is more to learn about assemblerr and you have multiple sources available:
?. For example,
?prm_log_normalwill open the documentation for the log-normal parameter building block. You can find links to similar building blocks at the bottom of each help page.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.