Estimating nlmixr2 models with 'nlmixr2targets'

library(nlmixr2targets)

Introduction to `nlmixr2targets`

The nlmixr2targets package improves reproducibility by ensuring that your model is up-to-date with your data, and it speeds your workflow using the targets package to only run models when the model or data have changed.

There are two main functions that are used within the package:

tar_nlmixr() which runs a single model, and
tar_nlmixr_multimodel() which runs multiple models for a single dataset.

Using nlmixr2targets requires the use of the targets package. To learn about the targets package, see the targets website.

Initial conditions

The native nlmixr2 DSL form for a compartment initial value is cmt(0) <- value inside a model({...}) block. Inside tar_nlmixr() and tar_nlmixr_multimodel() you may also write the nlmixr2targets-only workaround cmt(initial) <- value, which is rewritten back to cmt(0) <- value before nlmixr2 ever sees the model. The cmt(initial) form is not understood by bare nlmixr2 and only fits when routed through nlmixr2targets. Internally, tar_nlmixr() also rewrites cmt(0) to cmt(initial) in env so that targets’ static analysis can walk the model body, then restores cmt(0) at runtime.

See vignette("initial-conditions", package = "nlmixr2targets") for the full cheatsheet, including the pipe forms pheno |> model({...}) and pheno |> ini(...), and the known codetools::findGlobals() edge case for functions in env that are never routed through tar_nlmixr().

Running one model with one dataset (`tar_nlmixr()`)

The tar_nlmixr() function allows you to estimate one model with one dataset. It will generate three targets: a simplified version of the model, a minimal version of the dataset, and the estimation step.

The simplified version of the model removes parts that are less reproducible but changes none of the model intent. (Advanced information: The parts that are removed are that the source references and the model name. Also, the model is modified at this step for setting initial values as described in the previous section of this vignette.)

library(targets)
library(tarchetypes)
library(nlmixr2targets)

pheno <- function() {
  ini({
    lcl <- log(0.008); label("Typical value of clearance")
    lvc <- log(0.6); label("Typical value of volume of distribution")
    etalcl + etalvc ~ c(1,
                        0.01, 1)
    cpaddSd <- 0.1; label("residual variability")
  })
  model({
    cl <- exp(lcl + etalcl)
    vc <- exp(lvc + etalvc)
    kel <- cl / vc
    d / dt(central) <- -kel * central
    cp <- central / vc
    cp ~ add(cpaddSd)
  })
}

plan_model <-
  tar_plan(
    myData = nlmixr2data::pheno_sd,
    tar_nlmixr(
      model_pheno,
      object = pheno,
      data = myData,
      est = "saem"
    )
  )

list(
  plan_model
)

Running multiple models with one dataset (`tar_nlmixr_multimodel()`)

A common use case is to generate multiple models using a single dataset and estimation method. tar_nlmixr_multimodel() allows the generation of a named list of models to allow subsequent analysis of all models.

Internally, tar_nlmixr_multimodel() passes all the models to tar_nlmixr() so that the data set simplification and equivalent steps run once per model, and no model is run more often than required for dataset or model changes.

library(targets)
library(tarchetypes)
library(nlmixr2targets)

pheno <- function() {
  ini({
    lcl <- log(0.008); label("Typical value of clearance")
    lvc <- log(0.6); label("Typical value of volume of distribution")
    etalcl + etalvc ~ c(1,
                        0.01, 1)
    cpaddSd <- 0.1; label("residual variability")
  })
  model({
    cl <- exp(lcl + etalcl)
    vc <- exp(lvc + etalvc)
    kel <- cl / vc
    d / dt(central) <- -kel * central
    cp <- central / vc
    cp ~ add(cpaddSd)
  })
}

pheno2 <- function() {
  ini({
    lcl <- log(0.008); label("Typical value of clearance")
    lvc <- log(0.6); label("Typical value of volume of distribution")
    etalcl + etalvc ~ c(2,
                        0.01, 2)
    cpaddSd <- 3.0; label("residual variability")
  })
  model({
    cl <- exp(lcl + etalcl)
    vc <- exp(lvc + etalvc)
    kel <- cl / vc
    d / dt(central) <- -kel * central
    cp <- central / vc
    cp ~ add(cpaddSd)
  })
}

plan_model <-
  tar_nlmixr_multimodel(
    all_models,
    data = nlmixr2data::pheno_sd,
    est = "saem",
    "Base model; additive residual error = 1" = pheno,
    "Base model; additive residual error = 3" = pheno2
  )

plan_report <-
  tar_plan(
    # Determine the AIC for all tested models
    aic_list = sapply(X = all_models, FUN = AIC)
  )

list(
  plan_model,
  plan_report
)

Model piping for multiple models estimated with one dataset

Model piping for nlmixr2 models (see vignette("modelPiping", package = "nlmixr2")) is possible within the multiple models being estimated with tar_nlmixr_multimodel(). It simplifies examples like the one above so that you can focus on the model content and avoid rewriting models, as with all nlmixr2 model piping.

To use model piping, simply refer to the model by its name like a named list. Behind the scenes, nlmixr2targets will work out the dependencies between the models and only rerun the dependent model if it or the dependent model changes.

library(targets)
library(tarchetypes)
library(nlmixr2targets)
library(nlmixr2)

pheno <- function() {
  ini({
    lcl <- log(0.008); label("Typical value of clearance")
    lvc <- log(0.6); label("Typical value of volume of distribution")
    etalcl + etalvc ~ c(1,
                        0.01, 1)
    cpaddSd <- 0.1; label("residual variability")
  })
  model({
    cl <- exp(lcl + etalcl)
    vc <- exp(lvc + etalvc)
    kel <- cl / vc
    d / dt(central) <- -kel * central
    cp <- central / vc
    cp ~ add(cpaddSd)
  })
}

plan_model <-
  tar_nlmixr_multimodel(
    all_models,
    data = nlmixr2data::pheno_sd,
    est = "saem",
    "Base model; additive residual error = 1" = pheno,
    "Base model; additive residual error = 3" =
    all_models[["Base model; additive residual error = 1"]] |> ini(cpaddSd = 3)
  )

list(
  plan_model
)

Continuing the pipeline when a model fails

By default, if a model fails during estimation the error propagates and targets::tar_make() stops, just as any other target error would. This is usually what you want for a single model, but when you are fitting many models at once (for example with tar_nlmixr_multimodel()) one failing model would otherwise halt the whole pipeline and prevent you from seeing the models that did succeed.

Both tar_nlmixr() and tar_nlmixr_multimodel() accept an error argument to control this:

error = "stop" (the default) lets the estimation error propagate and halt the pipeline.
error = "continue" catches the estimation error and stores a failure sentinel on the target instead, so the rest of the pipeline still runs.

The sentinel is an object of class nlmixr2targetsError that also inherits from "try-error", and it carries the original error message. Because it is clearly not an nlmixr2 fit object, you can detect a failed model with a simple inherits() check.

library(targets)
library(tarchetypes)
library(nlmixr2targets)

plan_model <-
  tar_nlmixr_multimodel(
    all_models,
    data = nlmixr2data::pheno_sd,
    est = "saem",
    error = "continue",
    "Base model; additive residual error = 1" = pheno,
    "Base model; additive residual error = 3" = pheno2
  )

plan_report <-
  tar_plan(
    # Keep only the models that estimated successfully
    successful_models = Filter(
      f = function(fit) !inherits(fit, "try-error"),
      x = all_models
    ),
    # Compute AIC for the successful models only
    aic_list = sapply(X = successful_models, FUN = AIC)
  )

list(
  plan_model,
  plan_report
)

After tar_make(), a target that failed to estimate can be inspected directly:

fit <- tar_read(all_models)[["Base model; additive residual error = 3"]]
inherits(fit, "try-error")  # TRUE if this model failed
print(fit)                  # shows the captured error message

Introduction to nlmixr2targets