nlmixr2/rxode2 user functions to modify code

By Matthew Fidler in rxode2 nlmixr2

January 24, 2025

This month I will talk about a new type of user function. Previously, I spoke of user functions that you can use in your code to extend the functionality of rxode2 and nlmixr2.

Recently we released the ability for certain functions to generate code.

I will go over examples that could be helpful to extend rxode2 ui models:

  • An example that allows arguments to be named inside of the rxode2 ui models (though classic models still do not allow this).

  • An example that allows linear models to be expanded and possibly even automatically initialized.

If you are interested, this is what was added to allow neural network ODE to be integrated directly in rxode2/nlmixr2 using pmxNODE.

As a note, the user interface/function models are also flexible in generating specific code and is also used to allow posologyr (here) to use nlmixr2/rxode2 models directly for all of its clinical optimization routines. The model flexibility is also what is used to generate code for NONMEM, Monolix and PopED in babelmixr2.

In addition to this flexibility there is still a undocumented method to register low level C functions from another package to extend rxode2/nlmixr2 even further (which I will add in further documentation or a blog post).

The collection of these features and programming interface is unique to nlmixr2 and rxode2 (compared to say Monolix or NONMEM). As in posologyr and pmxNODE, this allows the language to be expanded in so many interesting and unique ways.

Also, this blog post is a bit more focused on how to add code expansion to a package, though it can work with your own code too; for many users this is simply a way to add Neural Network ODEs (after loading pmxNODE) to rxode2/nlmixr2, linear model expansions or other similar sorts of functions.

If you are not interested in the details for using this directly in your code, but are interested in ramifications of generating code from the UI, checkout the following package websites (in order of UI integration):

  • babelmixr2 generates PopED, NONMEM, and Monolix code from a nlmixr2/rxode2 model (as well as NONMEM/Monolix model imports). It also uses the interface to automatically start population models using PK estimates generated from PKNCA.

  • posologyr generates lower level rxode2 code and uses internal functions to personalize dosing regimens

  • pmxNODE which allows nlmixr2 to use neural network ODEs inside of rxode2 and nlmixr2 functions.

rxode2/nlmixr2 Extensions don’t need to use the ui features to be useful generally. Some good projects that extend nlmixr2 are (also in order of appearance):

  • shinyMixR for a shiny run manager for nlmixr2 models.

  • PharmTeX for generating pdf reports from a nlmixr2 model.

  • PharmPy for generating (some) models with nlmixr2.

  • nlmixr2rpt for generating Word/PowerPoint reports from a nlmixr2 model.

  • ruminate which includes a shiny simulator for rxode2.

  • nlmixr2autoinit automated initial estimates for Pop PK analyses using nlmixr2

  • nlmixr2auto automated model selection for Pop PK analyses using nlmixr2.

Let me know if I am missing any excellent tools!

There are other tools that also translate to nlmixr2 (which I am not including at this time, but perhaps later).

Now, back to adding new user functions that generate code in rxode2/nlmixr2.

Functions to insert rxode2 code into the current model

The first example is allowing named arguments in a function that translates itself into a non-named form:

f <- function() {
  model({
    a <- rxpois(lambda=lam)
  })
}

# Which will evaluate into a standard rxode2 function that does not
# support named arguments (since it is translated to C)
f()

# Which is still true in the standard rxode2:

try(rxode2({
  a <- rxpois(lambda=lam)
}))

This is accomplished by a combination of two functions, which are highly commented:

rxUdfUi.rxpois <- function(fun) {
  # Fun is the language object (ie quoted R object) to be evaluated or
  # changed in the code
  .fun <- fun
  # Since the `rxpois` function is built into the rxode2 we need to
  # have a function with a different conflicts.  In this case, I take
  # the function name (fun[[1]]), and prepend a ".", which follows
  # `rxode2`'s naming convention of un-exported functions.
  #
  # This next evaluation changes the expression function to .rxpois()
  .fun[[1]] <- str2lang(paste0(".", deparse1(fun[[1]])))
  # Since this is still a R expression, you can then evaluate the
  # function .rxpois to produce the proper code:
  eval(.fun)
}

# The above s3 method can be registered in a package or you can use
# the following code to register it in your session:
rxode2::.s3register("rxode2::rxUdfUi", "rxpois")

# This is the function that changes the code as needed
.rxpois <- function(lambda) {
  # The first part of this code tries to change the value into a
  # character.  This handles cases like rxpois(lambda=lam),
  # rxpois(lam), rxpois("lam").  It also tries to evaluate the
  # argument supplied to lambda in case it comes from a different
  # location.
  .lam <- as.character(substitute(lambda))
  .tmp <- try(force(lambda), silent=TRUE)
  if (!inherits(.tmp, "try-error")) {
    if (is.character(.tmp)) {
      .lam <- lambda
    }
  }
  # This part creates a list with the replacement text, in this case
  # it woulb be rxpois(lam) where there is no equals included, as
  # required by `rxode2`:
  list(replace=paste0("rxpois(", .lam, ")"))
}

In general the list that the function needs to return can have:

  • $replace – The text that will be replaced

  • $before – lines that will be placed in the model before the current function is found

  • $after – lines that are added in the model after the current function is found

  • $iniDf – the initial estimates data.frame for this problem. While calling this function, you can retrieve the initial conditions currently used parsing you can get the prior value with rxUdfUiIniDf() and then you can modify it inside the function and return the new data.frame in this list element. This allows you add/delete initial estimates from the model as well as modify the model lines themselves.

  • $uiUseData – when TRUE, this instructs rxode2 and nlmixr2est to re-parse this function in the presence of data, this means a bit more function setup will need to be done.

  • $uiUseMv – when TRUE this instructs rxode2 to re-parse the function after the initial model variables are calculated.

In addition to the rxUdfUiIniDf() you can get information about the parser:

  • rxUdfUiParsing() returns if the rxode2 ui function is being parsed currently (this allows a function to be overloaded as a udf for calling from rxode2 as well as a function for modifying the model).

  • rxUdfUiNum() during parsing the function you are calling (in the example above rxpois() can be called multiple times. This gives the number of the function in the model in order (the first would give 1 the second, 2, etc). This can be used to create unique variables with functions like rxIntToLetter() or rxIntToBase().

  • rxUdfUiIniLhs() which gives the left-handed side of the equation where the function is found. This is also a R language object.

  • rxUdfUiIniMv() gives the model variables for parsing (can be used in functions like linCmt())

  • rxUdfUiData() which specifies the data that are being used to simulate, estimate, etc.

  • rxUdfUiEst() which gives the estimation/simulation method that is being used with the model. For example, with simulation it would be rxSolve.

Using model variables in rxode2 ui models

You can also take and change the model and take into consideration the rxode2 model variables before the full ui has completed its parsing. These rxode2 model variables has information that might change what variables you make or names of variables. For example it has what is on the left hand side of the equations ($lhs), what are the input parameters ($params) and what is the ODE states ($state)).

If you are using this approach, you will likely need to do the following steps:

  • When data are not being processed, you need to put the function in an rxode2 acceptable form, no named arguments, no strings, and only numbers or variables in the output.

  • The number of arguments of this output needs to be declared in the S3 method by adding the attribute "nargs" to method. For example, the built in testMod1() ui modification function uses only one argument when parsed

Below is a commented example of the model variables example:

testMod1 <- function(val=1) {
  # This converts the val to a character if it is somthing like testMod1(b)
  .val <- as.character(substitute(val))
  .tmp <- suppressWarnings(try(force(val), silent = TRUE))
  if (!inherits(.tmp, "try-error")) {
    if (is.character(.tmp)) {
      .val <- val
    }
  }
  # This does the UI parsing
  if (rxUdfUiParsing()) {
    # See if the model variables are available
    .mv <- rxUdfUiMv()
    if (is.null(.mv)) {
      # Put this in a rxode2 low level acceptible form, no complex
      # expressions, no named arguments, something that is suitable
      # for C.
      #
      # The `uiUsMv` tells the parser this needs to be reparsed when
      # the model variables become avaialble during parsing.
      return(list(replace=paste0("testMod1(", .val, ")"),
                  uiUseMv=TRUE))
    } else {
      # Now that we have the model variables, we can then do something
      # about this
      .vars <- .mv$params
      if (length(.vars) > 0) {
        # If there is parameters available, this dummy function times
        # the first input function by the value specified
        return(list(replace=paste0(.vars[1], "*", .val)))
      } else {
        # If the value isn't availble, simply replace the function
        # with the value.
        return(list(replace=.val))
      }
    }
  }
  stop("This function is only for use in rxode2 ui models",
       call.=FALSE)

}

rxUdfUi.testMod1 <- function(fun) {
  eval(fun)
}

# To allow this to go to the next step, you need to declare how many
# arguments this argument has, in this case 1.  Bu adding the
# attribute "nargs", rxode2 lower level parser knows how to handle
# this new function.  This allows rxode2 to generate the model
# variables and send it to the next step.

attr(rxUdfUi.testMod1, "nargs") <- 1L

# If you are in a package, you can use the rxoygen tag @export to
# register this as a rxode2 model definition.
#
# If you are using this in your own script, you need to register the s3 function
# One way to do this is:
rxode2::.s3register("rxode2::rxUdfUi", "testMod1")

## These are some examples of this function in use:

f <- function() {
  model({
    a <- b + testMod1(3)
  })
}

f <- f()

print(f)
##  ── rxode2-based Pred model ───────────────────────────────────────────────────── 
##  ── Model (Normalized Syntax): ── 
## function() {
##     model({
##         a <- b + (b * 3)
##     })
## }
f <- function() {
  model({
    a <- testMod1(c)
  })
}

f <- f()

print(f)
##  ── rxode2-based Pred model ───────────────────────────────────────────────────── 
##  ── Model (Normalized Syntax): ── 
## function() {
##     model({
##         a <- (c * c)
##     })
## }
f <- function() {
  model({
    a <- testMod1(1)
  })
}

f <- f()

print(f)
##  ── rxode2-based Pred model ───────────────────────────────────────────────────── 
##  ── Model (Normalized Syntax): ── 
## function() {
##     model({
##         a <- 1
##     })
## }

Using data for rxode2 ui modification models

The same steps are needed to use the data in the model replacement; You can then use the data and the model to replace the values inside the model. A worked example linMod() is included that has the ability to use:

  • model variables,
  • put lines before or after the model,
  • add initial conditions
  • And use data in the initial estimates
# You can print the code:
linMod
## function (variable, power, dv = "dv", intercept = TRUE, type = c("replace", 
##     "before", "after"), num = NULL, iniDf = NULL, data = FALSE, 
##     mv = FALSE) 
## {
##     .dv <- as.character(substitute(dv))
##     .tmp <- suppressWarnings(try(force(dv), silent = TRUE))
##     if (!inherits(.tmp, "try-error")) {
##         if (is.character(.tmp)) {
##             .dv <- dv
##         }
##     }
##     .var <- as.character(substitute(variable))
##     .tmp <- try(force(variable), silent = TRUE)
##     .doExp3 <- FALSE
##     if (!inherits(.tmp, "try-error")) {
##         if (is.character(.tmp)) {
##             .var <- variable
##         }
##         else if (!inherits(.tmp, "formula")) {
##             .dv <- as.character(substitute(dv))
##             .tmp <- suppressWarnings(try(force(dv), silent = TRUE))
##             if (!inherits(.tmp, "try-error")) {
##                 if (is.character(.tmp)) {
##                   .dv <- dv
##                 }
##             }
##         }
##         else if (length(variable) == 2) {
##             if (!identical(variable[[1]], quote(`~`))) {
##                 stop("unexpected formula, needs to be the form ~x^3", 
##                   call. = FALSE)
##             }
##             .doExp3 <- TRUE
##             .exp3 <- variable[[2]]
##         }
##         else {
##             if (length(variable) != 3) {
##                 stop("unexpected formula, needs to be the form dv~x^3", 
##                   call. = FALSE)
##             }
##             if (!identical(variable[[1]], quote(`~`))) {
##                 stop("unexpected formula, needs to be the form dv~x^3", 
##                   call. = FALSE)
##             }
##             .dv <- as.character(variable[[2]])
##             data <- TRUE
##             .exp3 <- variable[[3]]
##             .doExp3 <- TRUE
##         }
##         if (.doExp3) {
##             if (length(.exp3) == 1) {
##                 .var <- variable <- as.character(.exp3)
##                 power <- 1
##             }
##             else if (length(.exp3) == 3) {
##                 if (!identical(.exp3[[1]], quote(`^`))) {
##                   stop("unexpected formula, needs to be the form dv~x^3", 
##                     call. = FALSE)
##                 }
##                 if (!is.numeric(.exp3[[3]])) {
##                   stop("unexpected formula, needs to be the form dv~x^3", 
##                     call. = FALSE)
##                 }
##                 .var <- variable <- as.character(.exp3[[2]])
##                 power <- .exp3[[3]]
##             }
##             else {
##                 stop("unexpected formula, needs to be the form dv~x^3", 
##                   call. = FALSE)
##             }
##         }
##     }
##     checkmate::assertCharacter(.var, len = 1L, any.missing = FALSE, 
##         pattern = "^[.]*[a-zA-Z]+[a-zA-Z0-9._]*$", min.chars = 1L, 
##         .var.name = "variable")
##     checkmate::assertCharacter(.dv, len = 1L, any.missing = FALSE, 
##         pattern = "^[.]*[a-zA-Z]+[a-zA-Z0-9._]*$", min.chars = 1L, 
##         .var.name = "dv")
##     checkmate::assertLogical(intercept, len = 1L, any.missing = FALSE)
##     checkmate::assertIntegerish(power, lower = ifelse(intercept, 
##         0L, 1L), len = 1L)
##     if (is.null(num)) {
##         num <- rxUdfUiNum()
##     }
##     checkmate::assertIntegerish(num, lower = 1, any.missing = FALSE, 
##         len = 1)
##     if (mv && is.null(rxUdfUiMv())) {
##         if (intercept) {
##             return(list(replace = paste0("linModM(", .var, ", ", 
##                 power, ")"), uiUseMv = TRUE))
##         }
##         else {
##             return(list(replace = paste0("linModM0(", .var, ", ", 
##                 power, ")"), uiUseMv = TRUE))
##         }
##     }
##     if (data && is.null(rxUdfUiData())) {
##         if (intercept) {
##             return(list(replace = paste0("linModD(", .var, ", ", 
##                 power, ", ", .dv, ")"), uiUseData = TRUE))
##         }
##         else {
##             return(list(replace = paste0("linModD0(", .var, ", ", 
##                 power, ",", .dv, ")"), uiUseData = TRUE))
##         }
##     }
##     if (is.null(iniDf)) {
##         iniDf <- rxUdfUiIniDf()
##     }
##     assertIniDf(iniDf, null.ok = TRUE)
##     type <- match.arg(type)
##     .mv <- rxUdfUiMv()
##     if (!is.null(.mv)) {
##         .varsMv <- c(.mv$lhs, .mv$params, .mv$state)
##         .pre <- paste0(.var, num, rxIntToLetter(seq_len(power + 
##             ifelse(intercept, 1L, 0L)) - 1L))
##         .pre <- vapply(.pre, function(v) {
##             if (v %in% .varsMv) {
##                 paste0("rx.linMod.", v)
##             }
##             else {
##                 v
##             }
##         }, character(1), USE.NAMES = FALSE)
##     }
##     else {
##         .pre <- paste0("rx.linMod.", .var, num, rxIntToLetter(seq_len(power + 
##             ifelse(intercept, 1L, 0L)) - 1L))
##     }
##     if (!is.null(iniDf)) {
##         .theta <- iniDf[!is.na(iniDf$ntheta), , drop = FALSE]
##         if (length(.theta$ntheta) > 0L) {
##             .maxTheta <- max(.theta$ntheta)
##             .theta1 <- .theta[1, ]
##         }
##         else {
##             .maxTheta <- 0L
##             .theta1 <- .rxBlankIni("theta")
##         }
##         .theta1$lower <- -Inf
##         .theta1$upper <- Inf
##         .theta1$fix <- FALSE
##         .theta1$label <- NA_character_
##         .theta1$backTransform <- NA_character_
##         .theta1$condition <- NA_character_
##         .theta1$err <- NA_character_
##         .est <- rep(0, length(.pre))
##         if (data) {
##             .dat <- rxUdfUiData()
##             .wdv <- which(tolower(names(.dat)) == tolower(.dv))
##             if (length(.wdv) == 0L) {
##                 warning(.dv, "not found in data, so no initial estimates will be set to zero")
##             }
##             else {
##                 names(.dat)[.wdv] <- .dv
##                 .model <- stats::lm(stats::as.formula(paste0(.dv, 
##                   " ~ stats::poly(", .var, ",", power, ")", ifelse(intercept, 
##                     "", "+0"))), data = rxUdfUiData())
##                 .est <- coef(.model)
##             }
##         }
##         .cur <- c(list(.theta), lapply(seq_along(.pre), function(i) {
##             .cur <- .theta1
##             .cur$name <- .pre[i]
##             .cur$est <- .est[i]
##             .cur$ntheta <- .maxTheta + i
##             .cur
##         }))
##         .theta <- do.call(rbind, .cur)
##         .eta <- iniDf[is.na(iniDf$neta), , drop = FALSE]
##         .iniDf <- rbind(.theta, .eta)
##     }
##     else {
##         .iniDf <- NULL
##     }
##     .linMod <- paste(vapply(seq_along(.pre), function(i) {
##         if (intercept) {
##             if (i == 1) 
##                 return(.pre[i])
##             if (i == 2) 
##                 return(paste0(.pre[i], "*", .var))
##             paste0(.pre[i], "*", paste0(.var, "^", i - 1L))
##         }
##         else {
##             if (i == 1) 
##                 return(paste0(.pre[i], "*", .var))
##             paste0(.pre[i], "*", paste0(.var, "^", i))
##         }
##     }, character(1)), collapse = "+")
##     if (type == "replace") {
##         list(replace = .linMod, iniDf = .iniDf)
##     }
##     else if (type == "before") {
##         .replace <- paste0("rx.linMod.", .var, ".f", num)
##         list(before = paste0(.replace, " <- ", .linMod), replace = .replace, 
##             iniDf = .iniDf)
##     }
##     else if (type == "after") {
##         .replace <- paste0("rx.linMod.", .var, ".f", num)
##         list(after = paste0(.replace, " <- ", .linMod), replace = "0", 
##             iniDf = .iniDf)
##     }
## }
## <bytecode: 0x650dd9f133f8>
## <environment: namespace:rxode2>
# You can also print the s3 method that is used for this method

rxode2:::rxUdfUi.linMod
## function (fun) 
## {
##     eval(fun)
## }
## <bytecode: 0x650dd9f96dc8>
## <environment: namespace:rxode2>
## attr(,"nargs")
## [1] 2
Posted on:
January 24, 2025
Length:
14 minute read, 2860 words
Categories:
rxode2 nlmixr2
See Also: