9 More complex models

9.1 Set up the tutorial environment

Tip

Make sure you have run the tutorial setup code in your R session before copying, pasting and running example code here.

9.2 Introduction

Here we provide examples of adding covariates and compartments to models.

9.3 Model with covariates

Let’s create another run object, this time using a model that include covariates.

First, recall the model from the NPAG example.

# model with covariates declared but not used
mod1 <- PM_model$new(
  pri = list(
    Ka = ab(0.1, 0.9),
    Ke = ab(0.001, 0.1),
    V = ab(30, 120),
    lag1 = ab(0, 4)
  ),
  cov = list(
    wt = interp(),
    africa = interp("none"),
    age = interp(),
    gender = interp("none"),
    height = interp()
  ),
  eqn = function(){
    two_comp_bolus
  },
  lag = function(){
    lag[1] = lag1
  },
  out = function(){
  Y[1] = X[2]/V
  },
  err = list(
    proportional(5, c(0.02, 0.05, -0.0002, 0))
  )
)

Next, we update it to change V to an new variable called V0 in the PRI list element. In the OUT list element, we include an equation that relates V0 to normalized weight. Note that we could also have included the V = V0 * (wt/70) in a SEC block to globally define this relationship for any other model list block. In this case, we only use V in the OUT block, so we define it there locally.

# declare covariates AND use them in the model
mod2 <- PM_model$new(
  pri = list(
    Ka = ab(0.1, 0.9),
    Ke = ab(0.001, 0.1),
    V0 = ab(30, 120),
    lag1 = ab(0, 4)
  ),
  cov = list(
    wt = interp(),
    africa = interp("none"),
    age = interp(),
    gender = interp("none"),
    height = interp()
  ),
  eqn = function(){
    two_comp_bolus
  },
  lag = function(){
    lag[1] = lag1
  },
  out = function(){
    V = V0 * (wt/70)
    Y[1] = X[2]/V
  },
  err = list(
    proportional(5, c(0.02, 0.05, -0.0002, 0))
  )
)

Ensure the data are loaded.

exData <- PM_data$new(data = "src/ex.csv")

Now, run that new model.

run2 <- mod2$fit(data = exData, run = 2, path = "Runs", overwrite = TRUE)

Right after the run, the results are in run2, but for future loading if you return to your script later, you can load it with PM_load().

run2 <- PM_load(2, path = "Runs")

You can examine the results as we’ve seen in previous chapters.

9.4 Add compartments to a model

Now, let’s create a model with an additional compartment. We’ll add a peripheral compartment to the previous two-compartment bolus model.

# model with covariates and an additional compartment
mod3 <- PM_model$new(
  pri = list(
    Ka = ab(0.1, 0.9),
    Ke = ab(0.001, 0.1),
    V0 = ab(30, 120),
    KCP = ab(0.001, 0.5),
    KPC = ab(0.001, 0.5),
    lag1 = ab(0, 4)
  ),
  cov = list(
    wt = interp(),
    africa = interp("none"),
    age = interp(),
    gender = interp("none"),
    height = interp()
  ),
  eqn = function(){
    three_comp_bolus
  },
  lag = function(){
    lag[1] = lag1
  },
  out = function(){
    V = V0 * (wt/70)
    Y[1] = X[2]/V
  },
  err = list(
    proportional(5, c(0.02, 0.05, -0.0002, 0))
  )
)

Notes:

We used model_lib to check the name of the algebraic 3-compartment model and the required parameters that should be declared (Ka, Ke, KCP, KPC and V).
Don’t forget that you can simply type the name of the library model into the console to get information about it.

# type the name to get information about the model
three_comp_bolus

Output is still in compartment 2.

Now, run that new model.

# run model with 3 compartments and covariates
run3 <- mod3$fit(data = exData, run = 3, path = "Runs", overwrite = TRUE)

Again, you can load it later with PM_load().

run3 <- PM_load(3, path = "Runs")

9.5 Differential equations

Here is mod3 again, but this time using differential equations instead of built-in compartmental models, just to show the structure of the ODEs.

# ODE model with covariates and an additional compartment
mod3b <- PM_model$new(
  pri = list(
    Ka = ab(0.1, 0.9),
    Ke = ab(0.001, 0.1),
    V0 = ab(30, 120),
    K12 = ab(0.001, 0.5),
    K21 = ab(0.001, 0.5),
    lag1 = ab(0, 4)
  ),
  cov = list(
    wt = interp(),
    africa = interp("none"),
    age = interp(),
    gender = interp("none"),
    height = interp()
  ),
  eqn = function(){
    dX[1] = B[1] - Ka * X[1]
    dX[2] = R[1] + Ka * X[1] - (K12 + Ke) * X[2] + K21 * X[3]
    dX[3] = K12 * X[2] - K21 * X[3]
  },
  lag = function(){
    lag[1] = lag1
  },
  out = function(){
    V = V0 * (wt/70)
    Y[1] = X[2]/V
  },
  err = list(
    proportional(5, c(0.02, 0.05, -0.0002, 0))
  )
)

Notes:

We have replaced KCP and KPC with K12 and K21, to emphasize that ODE models do not require specifically named parameters in the PRI block.
By including R[1] in the central compartment, X[2] in this model, we allow for infusions of drug 1 in addition to the bolus input for drug 1, B[1]. The equivalent algebraic model, three_comp_bolus, does the same.