Skip to contents

Low-level models and other input styles

Source: vignettes/articles/rxode2-classic-models.Rmd
rxode2-classic-models.Rmd
library(rxode2)
#> rxode2 5.1.8.9000 using 2 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

Most of the rxode2 documentation writes models as model functions, with an ini({}) block for the parameter values and a model({}) block for the equations, and builds dosing with et(). rxode2 also supports other styles, mostly kept from earlier versions and from other tools. This page collects them and shows how they behave differently.

Low-level models

Writing a low-level model

A low-level model gives rxode2() the model equations only, without ini({}). It can be written as a block of statements:

mod <- rxode2({
  ka <- exp(tka)
  cl <- exp(tcl)
  v <- exp(tv)
  d/dt(depot) <- -ka * depot
  d/dt(center) <- ka * depot - cl / v * center
  cp <- center / v
})
mod
#> rxode2 5.1.8.9000 model named rx_9b9374d1e585f14e028927a52e797350 model (✔ ready). 
#> x$state: depot, center
#> x$params: tka, tcl, tv
#> x$lhs: ka, cl, v, cp

as a character string:

modStr <- rxode2("
  d/dt(center) <- -cl / v * center
  cp <- center / v
")

or as a file:

modFile <- tempfile(fileext = ".rxode2")
writeLines(c("d/dt(center) <- -cl / v * center",
             "cp <- center / v"), modFile)
modFile <- rxode2(filename = modFile)
modFile
#> rxode2 5.1.8.9000 model named rx_4410256ee5a9f5195a432b1cf7b05480 model (✔ ready). 
#> x$state: center
#> x$params: cl, v
#> x$lhs: cp

Each call translates the model to C, compiles it and returns an rxode2 object.

Solving: params= and inits=

A low-level model has no parameter values, so every parameter must be supplied when solving, with params= (a named vector, or a data frame with one row per subject):

ev <- et(amt = 300) |> et(seq(0, 24, by = 2))
s <- rxSolve(mod, ev, params = c(tka = 0.45, tcl = 1, tv = 3.45))
head(s)
#>   time       ka       cl        v       cp        depot   center
#> 1    0 1.568312 2.718282 31.50039 0.000000 3.000000e+02   0.0000
#> 2    2 1.568312 2.718282 31.50039 8.042989 1.302872e+01 253.3573
#> 3    4 1.568312 2.718282 31.50039 7.117353 5.658247e-01 224.1994
#> 4    6 1.568312 2.718282 31.50039 6.004315 2.457319e-02 189.1383
#> 5    8 1.568312 2.718282 31.50039 5.053200 1.067188e-03 159.1778
#> 6   10 1.568312 2.718282 31.50039 4.252220 4.634651e-05 133.9466

A missing parameter is an error:

try(rxSolve(mod, ev))
#> Model:
#> 
#> ka=exp(tka);
#> cl=exp(tcl);
#> v=exp(tv);
#> d/dt(depot)=-ka*depot;
#> d/dt(center)=ka*depot-cl/v*center;
#> cp=center/v;
#> 
#> Error in rxSolveSEXP(object, .ctl, .nms, .xtra, params, .eventsForSolve,  : 
#>   The following parameter(s) are required for solving: tv, tcl, tka

The state values at time 0 can be given with inits= (or with state(0) <- in the model, as in a model function):

s <- rxSolve(mod, ev, params = c(tka = 0.45, tcl = 1, tv = 3.45),
             inits = c(depot = 10))
head(s)
#>   time       ka       cl        v       cp        depot   center
#> 1    0 1.568312 2.718282 31.50039 0.000000 3.100000e+02   0.0000
#> 2    2 1.568312 2.718282 31.50039 8.311088 1.346301e+01 261.8025
#> 3    4 1.568312 2.718282 31.50039 7.354598 5.846850e-01 231.6727
#> 4    6 1.568312 2.718282 31.50039 6.204459 2.539227e-02 195.4429
#> 5    8 1.568312 2.718282 31.50039 5.221640 1.102760e-03 164.4837
#> 6   10 1.568312 2.718282 31.50039 4.393961 4.789205e-05 138.4115

How they differ from model functions

The same model as a model function:

modFun <- function() {
  ini({
    tka <- 0.45
    tcl <- 1
    tv <- 3.45
    eta.cl ~ 0.1
    add.sd <- 0.7
  })
  model({
    ka <- exp(tka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv)
    d/dt(depot) <- -ka * depot
    d/dt(center) <- ka * depot - cl / v * center
    cp <- center / v
    cp ~ add(add.sd)
  })
}
ui <- modFun()

  • Parameter values. A model function takes its parameter values from ini({}), so it solves with rxSolve(ui, ev). params= still works and overrides the ini({}) values. A low-level model only has the values supplied with params=.

  • Between-subject variability. The eta.cl ~ 0.1 line in ini({}) is used to simulate the etas. A low-level model has no omega; it has to be given with omega= in rxSolve().

  • Residual error. cp ~ add(add.sd) defines an endpoint with residual error; solving simulates it in the sim column (ipredSim is the prediction). The low-level language does not have endpoints: there ~ is an assignment that is not included in the output, and residual error terms have to be written by hand with sigma= in rxSolve().

rxSolve(ui, ev, seed = 42)
#> ── Solved rxode2 object ──
#> ── Parameters (x$params): ──
#>        tka        tcl         tv     add.sd     eta.cl 
#> 0.45000000 1.00000000 3.45000000 0.70000000 0.07431574 
#> ── Initial Conditions (x$inits): ──
#>  depot center 
#>      0      0 
#> ── First part of data (object): ──
#> # A tibble: 13 × 9
#>    time    ka    cl     v    cp ipredSim   sim       depot center
#>   <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl> <dbl>       <dbl>  <dbl>
#> 1     0  1.57  2.93  31.5  0        0    0.815 300             0 
#> 2     2  1.57  2.93  31.5  7.97     7.97 7.46   13.0         251.
#> 3     4  1.57  2.93  31.5  6.96     6.96 6.62    0.566       219.
#> 4     6  1.57  2.93  31.5  5.80     5.80 5.47    0.0246      183.
#> 5     8  1.57  2.93  31.5  4.81     4.81 4.73    0.00107     152.
#> 6    10  1.57  2.93  31.5  4.00     4.00 3.21    0.0000463   126.
#> # ℹ 7 more rows
hidden <- rxode2({
  a ~ 3
  b <- a + 1
})
hidden$lhs # `a` is not output
#> [1] "b"
try(rxode2({
  cp <- center / v
  cp ~ add(add.sd)
}))
#> rxode2 model syntax error:
#> ================================================================================
#> :001:     cp <- center/v
#> :002: function 'add' is not supported; user function not found:
#>           cp ~ add(add.sd)
#>                    ^~~~~~
#> ================================================================================
#> Error : syntax errors (see above)
  • Turning the variability off. For a model function, omega = NA and sigma = NA give the typical-value prediction (in the pred column):
rxSolve(ui, ev, params = c(tcl = 2), omega = NA, sigma = NA)
#> ── Solved rxode2 object ──
#> ── Parameters (x$params): ──
#>      tka      tcl       tv   add.sd   eta.cl rxerr.cp 
#>     0.45     2.00     3.45     0.70     0.00     0.00 
#> ── Initial Conditions (x$inits): ──
#>  depot center 
#>      0      0 
#> ── First part of data (object): ──
#> # A tibble: 13 × 8
#>    time    ka    cl     v    cp  pred       depot center
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>       <dbl>  <dbl>
#> 1     0  1.57  7.39  31.5  0     0    300            0  
#> 2     2  1.57  7.39  31.5  6.52  6.52  13.0        205. 
#> 3     4  1.57  7.39  31.5  4.36  4.36   0.566      137. 
#> 4     6  1.57  7.39  31.5  2.74  2.74   0.0246      86.3
#> 5     8  1.57  7.39  31.5  1.71  1.71   0.00107     54.0
#> 6    10  1.57  7.39  31.5  1.07  1.07   0.0000463   33.8
#> # ℹ 7 more rows
  • Model piping. model() and ini() piping work on the model function form. Piping a low-level model first converts it to a model function:
mod |> ini(tka = 0.45, tcl = 1, tv = 3.45)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ promote `tka` to population parameter with initial estimate 0.45
#> ℹ change initial estimate of `tka` to `0.45`
#> ℹ promote `tcl` to population parameter with initial estimate 1
#> ℹ change initial estimate of `tcl` to `1`
#> ℹ promote `tv` to population parameter with initial estimate 3.45
#> ℹ change initial estimate of `tv` to `3.45`
#>  ── rxode2-based free-form 2-cmt ODE model ────────────────────────────────────── 
#>  ── Initalization: ──  
#> Fixed Effects ($theta): 
#>  tka  tcl   tv 
#> 0.45 1.00 3.45 
#> 
#> States ($state or $stateDf): 
#>   Compartment Number Compartment Name
#> 1                  1            depot
#> 2                  2           center
#>  ── Model (Normalized Syntax): ── 
#> function() {
#>     ini({
#>         tka <- 0.45
#>         tcl <- 1
#>         tv <- 3.45
#>     })
#>     model({
#>         ka = exp(tka)
#>         cl = exp(tcl)
#>         v = exp(tv)
#>         d/dt(depot) = -ka * depot
#>         d/dt(center) = ka * depot - cl/v * center
#>         cp = center/v
#>     })
#> }

Converting between forms

  • ui$simulationModel (or ui$simulationIniModel, which also includes the ini({}) values) gives the low-level model that rxode2 solves for a model function, including the generated code for the endpoints:
ui$simulationModel
#> rxode2 5.1.8.9000 model named rx_bcaf847b49c4f40092b7afc27fd1b419 model (✔ ready). 
#> x$state: depot, center
#> x$stateExtra: cp
#> x$params: tka, tcl, tv, add.sd, eta.cl, rxerr.cp
#> x$lhs: ka, cl, v, cp, ipredSim, sim
  • as.rxUi() converts a low-level model to a model function object (rxUi); it has no ini({}) values:
modUi <- as.rxUi(mod)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
modUi$iniDf
#>  [1] ntheta        neta1         neta2         name          lower        
#>  [6] est           upper         fix           label         backTransform
#> [11] condition     err          
#> <0 rows> (or 0-length row.names)
as.function(modUi)
#> function () 
#> {
#>     model({
#>         ka = exp(tka)
#>         cl = exp(tcl)
#>         v = exp(tv)
#>         d/dt(depot) = -ka * depot
#>         d/dt(center) = ka * depot - cl/v * center
#>         cp = center/v
#>     })
#> }
#> <environment: 0x5633d49f8bc8>

When a low-level model is still the right tool

  • Models in packages. rxPkg() and rxUse() store compiled low-level models in a package (see Pre-compiled models in packages).
  • Generated code. When a tool writes the model text, the string and file forms take it directly.
  • Options of the compiled model. Options such as calcJac= (a symbolic Jacobian, see stiff ODEs) are arguments of rxode2().
  • All parameters per subject. When every parameter comes from a data frame with params=, an ini({}) block adds nothing.
  • C-level functions. Packages that provide C functions to rxode2 models (custom C/C++ functions) document the low-level form.

Other event-table styles

The main way to build dosing and sampling is et() (see event tables). rxSolve() also accepts the styles below.

eventTable() methods

The classic RxODE interface creates an event table and then calls methods on it:

evt <- eventTable(amount.units = "mg", time.units = "hours")
evt$add.dosing(dose = 300, nbr.doses = 3, dosing.interval = 12)
evt$add.sampling(seq(0, 36, by = 4))
evt
#> ── EventTable with 11 records ──
#> 1 dosing records (see x$get.dosing(); add with add.dosing or et)
#> 10 observation times (see x$get.sampling(); add with add.sampling or et)
#> multiple doses in `addl` columns, expand with x$expand(); or etExpand(x)
#> ── First part of x: ──
#> # A tibble: 11 × 5
#>    time   amt  ii  addl evid         
#>     [h] <dbl> [h] <int> <evid>       
#>  1    0   300  12     2 1:Dose (Add) 
#>  2    0    NA  NA    NA 0:Observation
#>  3    4    NA  NA    NA 0:Observation
#>  4    8    NA  NA    NA 0:Observation
#>  5   12    NA  NA    NA 0:Observation
#>  6   16    NA  NA    NA 0:Observation
#>  7   20    NA  NA    NA 0:Observation
#>  8   24    NA  NA    NA 0:Observation
#>  9   28    NA  NA    NA 0:Observation
#> 10   32    NA  NA    NA 0:Observation
#> 11   36    NA  NA    NA 0:Observation

The same methods are functions that can be piped:

evt2 <- eventTable() |>
  add.dosing(dose = 300, nbr.doses = 3, dosing.interval = 12) |>
  add.sampling(seq(0, 36, by = 4))
rxSolve(ui, evt2, omega = NA, sigma = NA)
#> ── Solved rxode2 object ──
#> ── Parameters (x$params): ──
#>      tka      tcl       tv   add.sd   eta.cl rxerr.cp 
#>     0.45     1.00     3.45     0.70     0.00     0.00 
#> ── Initial Conditions (x$inits): ──
#>  depot center 
#>      0      0 
#> ── First part of data (object): ──
#> # A tibble: 10 × 8
#>    time    ka    cl     v    cp  pred     depot center
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>     <dbl>  <dbl>
#> 1     0  1.57  2.72  31.5  0     0    300           0 
#> 2     4  1.57  2.72  31.5  7.12  7.12   0.566     224.
#> 3     8  1.57  2.72  31.5  5.05  5.05   0.00107   159.
#> 4    12  1.57  2.72  31.5  3.58  3.58 300.        113.
#> 5    16  1.57  2.72  31.5  9.65  9.65   0.566     304.
#> 6    20  1.57  2.72  31.5  6.85  6.85   0.00107   216.
#> # ℹ 4 more rows

NONMEM-style data frames

A data frame with NONMEM columns (ID, TIME, AMT, EVID, CMT, RATE, II, ADDL, SS, DV, …) is used directly (see dataset columns):

nmData <- data.frame(ID = 1,
                     TIME = c(0, 1, 2, 4, 8, 12),
                     AMT = c(300, 0, 0, 0, 0, 0),
                     EVID = c(1, 0, 0, 0, 0, 0),
                     CMT = 1)
rxSolve(ui, nmData, omega = NA, sigma = NA)
#> ── Solved rxode2 object ──
#> ── Parameters (x$params): ──
#>      tka      tcl       tv   add.sd   eta.cl rxerr.cp 
#>     0.45     1.00     3.45     0.70     0.00     0.00 
#> ── Initial Conditions (x$inits): ──
#>  depot center 
#>      0      0 
#> ── First part of data (object): ──
#> # A tibble: 5 × 8
#>    time    ka    cl     v    cp  pred       depot center
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>       <dbl>  <dbl>
#> 1     1  1.57  2.72  31.5  7.14  7.14 62.5          225.
#> 2     2  1.57  2.72  31.5  8.04  8.04 13.0          253.
#> 3     4  1.57  2.72  31.5  7.12  7.12  0.566        224.
#> 4     8  1.57  2.72  31.5  5.05  5.05  0.00107      159.
#> 5    12  1.57  2.72  31.5  3.58  3.58  0.00000201   113.

deSolve-style event data frames

A deSolve event data frame (columns var, time, value and method, with method one of "add", "mult" or "rep") is converted by et() to dose (evid=1), multiply (evid=6) and replace (evid=5) events:

eventdat <- data.frame(var = c("depot", "center", "center"),
                       time = c(0, 6, 12),
                       value = c(300, 2, 0),
                       method = c("add", "mult", "rep"))
ev3 <- et(eventdat) |> et(seq(0, 24, by = 2))
rxSolve(ui, ev3, omega = NA, sigma = NA)
#> ── Solved rxode2 object ──
#> ── Parameters (x$params): ──
#>      tka      tcl       tv   add.sd   eta.cl rxerr.cp 
#>     0.45     1.00     3.45     0.70     0.00     0.00 
#> ── Initial Conditions (x$inits): ──
#>  depot center 
#>      0      0 
#> ── First part of data (object): ──
#> # A tibble: 13 × 8
#>    time    ka    cl     v    cp  pred       depot center
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>       <dbl>  <dbl>
#> 1     0  1.57  2.72  31.5  0     0    300             0 
#> 2     2  1.57  2.72  31.5  8.04  8.04  13.0         253.
#> 3     4  1.57  2.72  31.5  7.12  7.12   0.566       224.
#> 4     6  1.57  2.72  31.5 12.0  12.0    0.0246      378.
#> 5     8  1.57  2.72  31.5 10.1  10.1    0.00107     318.
#> 6    10  1.57  2.72  31.5  8.50  8.50   0.0000463   268.
#> # ℹ 7 more rows

Classic RxODE evid codes

rxode2 also accepts the classic RxODE event ids, which combine the compartment, the dose type and the infusion flags in one number. They are described in classic RxODE events.

Next steps