Getting started

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 :

Assemble a simple model

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:

Line 2 and 3 define the parameters “emax” and “ed50” as being log-normally distributed with a specific median and variance on the log-scale.
Line 4 defines the algebraic relationship for the variable effect using a classical R formula. In these formulae, the ~ should be understood like a =. You can use formulae in assemblerr whenever you want to define a mathematical relationship.
Line 5 defines the observational model, i.e., what is the observation in the model (the variable “effect” in this case) as well as the model for residual error (additive) and the variance of the error term.

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.

Check and modify

When you are done defining your model you can check whether it contains any issues using the check() function, e.g.,

check(m)

## ! 1 critical issue

##   1. Undefined variable 'dose' in algebraics

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:

check(m)

## ✓ No 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)

## Warning: Building block replaced
## x The building block [emax: log-normal] has been replaced with [emax: no variability]

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:

render(m)

## $PROBLEM
## $INPUT DOSE ID DV
## $DATA data.csv IGNORE=@
## $PRED
## EMAX = THETA(1)
## ED50 = THETA(2) * EXP(ETA(1))
## EFFECT = EMAX * DOSE/(ED50 + DOSE)
## Y = EFFECT + EPS(1)
## $ESTIMATION METHOD=COND MAXEVAL=999999 
## $THETA (0, 10, Inf) ; POP_EMAX
## $THETA (0, 50, Inf) ; POP_ED50
## $OMEGA 0.09 ; IIV_ED50
## $SIGMA 1; RUV_ADD

Rendering

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

render(m)

## $PROBLEM
## $INPUT DOSE ID DV
## $DATA data.csv IGNORE=@
## $PRED
## EMAX = THETA(1)
## ED50 = THETA(2) * EXP(ETA(1))
## EFFECT = EMAX * DOSE/(ED50 + DOSE)
## Y = EFFECT + EPS(1)
## $ESTIMATION METHOD=COND MAXEVAL=999999 
## $THETA (0, 10, Inf) ; POP_EMAX
## $THETA (0, 50, Inf) ; POP_ED50
## $OMEGA 0.09 ; IIV_ED50
## $SIGMA 1; RUV_ADD

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")

The 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
  )
)

## $PROBLEM
## $INPUT DOSE ID DV
## $DATA data.csv IGNORE=@
## $PRED
## EMAX = THETA(1)
## MU_1 = LOG(THETA(2))
## ED50 = EXP(MU_1 + ETA(1))
## EFFECT = EMAX * DOSE/(ED50 + DOSE)
## Y = EFFECT + EPS(1)
## $ESTIMATION METHOD=COND MAXEVAL=999999 
## $THETA (0, 10, Inf) ; POP_EMAX
## $THETA (0, 50, Inf) ; POP_ED50
## $OMEGA 0.09 ; IIV_ED50
## $SIGMA 1; RUV_ADD

Tasks

“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))

## $PROBLEM
## $INPUT DOSE ID DV
## $DATA data.csv IGNORE=@
## $PRED
## EMAX = THETA(1)
## ED50 = THETA(2) * EXP(ETA(1))
## EFFECT = EMAX * DOSE/(ED50 + DOSE)
## Y = EFFECT + EPS(1)
## $ESTIMATION METHOD=IMP AUTO=1 
## $COVARIANCE PRINT=E
## $THETA (0, 10, Inf) ; POP_EMAX
## $THETA (0, 50, Inf) ; POP_ED50
## $OMEGA 0.09 ; IIV_ED50
## $SIGMA 1; RUV_ADD

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()))

## $PROBLEM
## $INPUT DOSE ID DV
## $DATA data.csv IGNORE=@
## $PRED
## EMAX = THETA(1)
## ED50 = THETA(2) * EXP(ETA(1))
## EFFECT = EMAX * DOSE/(ED50 + DOSE)
## Y = EFFECT + EPS(1)
## $ESTIMATION METHOD=COND MAXEVAL=999999 
## $TABLE EMAX ED50 FILE=sdtab NOAPPEND NOPRINT
## $THETA (0, 10, Inf) ; POP_EMAX
## $THETA (0, 50, Inf) ; POP_ED50
## $OMEGA 0.09 ; IIV_ED50
## $SIGMA 1; RUV_ADD

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 tsk_output() function.

Ordinary differential equations

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:

render(m2)

## $PROBLEM
## $INPUT ID TIME DV AMT
## $DATA data.csv IGNORE=@
## $SUBROUTINES ADVAN1  TRANS2 
## $PK
## V = THETA(1) * EXP(ETA(1))
## CL = THETA(2) * EXP(ETA(2))
## EMAX = THETA(3) * EXP(ETA(3))
## EC50 = THETA(4) * EXP(ETA(4))
## $ERROR
## EFFECT = EMAX * (A(1)/V)/(EC50 + A(1)/V)
## Y = EFFECT + EPS(1)
## $ESTIMATION METHOD=COND MAXEVAL=999999 
## $THETA (0, 100, Inf) ; POP_V
## $THETA (0, 10, Inf) ; POP_CL
## $THETA (0, 10, Inf) ; POP_EMAX
## $THETA (0, 1, Inf) ; POP_EC50
## $OMEGA 0.09 ; IIV_V
## $OMEGA 0.09 ; IIV_CL
## $OMEGA 0.09 ; IIV_EMAX
## $OMEGA 0.09 ; IIV_EC50
## $SIGMA 1; RUV_ADD

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.

Pharmacokinetic models

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:

A first-order absorption model (line 2)
A compartment distribution model (line 3)
A linear elimination model (line 4)

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)
)

## $PROBLEM
## $INPUT ID TIME DV AMT
## $DATA data.csv IGNORE=@
## $SUBROUTINES ADVAN2  TRANS2 
## $PK
## MU_1 = LOG(THETA(1))
## MAT = EXP(MU_1 + ETA(1))
## MU_2 = LOG(THETA(2))
## VC = EXP(MU_2 + ETA(2))
## MU_3 = LOG(THETA(3))
## CL = EXP(MU_3 + ETA(3))
## KA = 1/MAT
## V = VC
## $ERROR
## CONC = A(2)/VC
## Y = CONC + EPS(1)
## $ESTIMATION METHOD=IMP AUTO=1 
## $TABLE MAT VC CL FILE=sdtab NOAPPEND NOPRINT
## $THETA (0, 0.5, Inf) ; POP_MAT
## $THETA (0, 100, Inf) ; POP_VC
## $THETA (0, 50, Inf) ; POP_CL
## $OMEGA 0.1 ; IIV_MAT
## $OMEGA 0.1 ; IIV_VC
## $OMEGA 0.1 ; IIV_CL
## $SIGMA 1; RUV_ADD

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)

## $PROBLEM
## $INPUT ID TIME DV AMT
## $DATA data.csv IGNORE=@
## $SUBROUTINES ADVAN2  TRANS2 
## $PK
## MAT = THETA(1) * EXP(ETA(1))
## V1 = THETA(2)
## CL = THETA(3) * EXP(ETA(2))
## KA = 1/MAT
## V = V1
## $ERROR
## CONC = A(2)/V1
## Y = CONC + EPS(1)
## $ESTIMATION METHOD=COND MAXEVAL=999999 
## $THETA (0, 0.5, Inf) ; POP_MAT
## $THETA (0, 100, Inf) ; POP_V1
## $THETA (0, 50, Inf) ; POP_CL
## $OMEGA 0.1 ; IIV_MAT
## $OMEGA 0.1 ; IIV_CL
## $SIGMA 1; RUV_ADD

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"])

## Warning: Building block replaced
## x The building block [vc: log-normal] has been replaced with [vc: no variability]

render(pkm)

## $PROBLEM
## $INPUT ID TIME DV AMT
## $DATA data.csv IGNORE=@
## $SUBROUTINES ADVAN2  TRANS2 
## $PK
## MAT = THETA(1) * EXP(ETA(1))
## VC = THETA(2)
## CL = THETA(3) * EXP(ETA(2))
## KA = 1/MAT
## V = VC
## $ERROR
## CONC = A(2)/VC
## Y = CONC + EPS(1)
## $ESTIMATION METHOD=COND MAXEVAL=999999 
## $THETA (0, 0.5, Inf) ; POP_MAT
## $THETA (0, 100, Inf) ; POP_VC
## $THETA (0, 50, Inf) ; POP_CL
## $OMEGA 0.1 ; IIV_MAT
## $OMEGA 0.1 ; IIV_CL
## $SIGMA 1; RUV_ADD