Skip to contents

ETA and EPS resampling

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

Overview

rxode2 supports resampling of between-subject variability (ETAs) and residual variability (EPSs) within the model code itself using simeta() and simeps(). This allows you to implement rejection sampling, truncated distributions, or any custom sampling logic directly in the model definition – without post-processing the simulated output.

This is analogous to the SIMETA and SIMEPS functions available in NONMEM, though the implementation and the context in which they are used differs slightly between the tools.

Resampling ETAs with simeta()

A common use case is constraining a covariate derived from a lognormal distribution to a physiologically plausible range. For example, simulating body weight from a lognormal distribution but only accepting values between 60 and 80 kg:

mod <- function() {
  ini({
    eta.wt ~ 0.25 # SD ~ 0.5 on log scale: wide variability
  })
  model({
    wt  <- 70 * exp(eta.wt)   # individual weight (kg)
    i   <- 0                   # iteration counter (for safety)
    while ((wt < 60) || (wt > 80)) {
      i <- i + 1
      if (i > 100) break       # prevent infinite loops
      simeta()                 # resample all ETAs for this subject
      wt <- 70 * exp(eta.wt)
    }
  })
}

When simeta() is called, all ETA values for the current subject are re-drawn from the omega matrix (here from ini({})). The model then re-evaluates the downstream expressions. Because the resampling is done at the C level, parallel solving with OpenMP is fully supported.

e <- et(seq(0, 24, by = 4), id = 1:100)

f <- rxSolve(mod, e, seed = 1234)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments

summary(f$wt)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   60.49   65.30   69.51   69.78   74.90   79.97

All simulated weights fall within the 60-80 kg window:

stopifnot(all(f$wt >= 60), all(f$wt <= 80))
range(f$wt)
#> [1] 60.48968 79.97034

Interaction with nStud

simeta() works correctly when simulating multiple studies (nStud). Each study draws its own omega matrix (if dfSub is specified for a scaled chi-squared draw), and simeta() resamples within the study-specific omega:

f2 <- rxSolve(
  mod, e,
  nStud  = 5,
  dfSub  = 40,
  seed   = 42
)

range(f2$wt)
#> [1] 60.01110 79.97703

Resampling EPSs with simeps()

simeps() re-draws all EPS (residual variability) values for the current observation record. A practical use case is ensuring a simulated observation is non-negative without post-hoc truncation:

mod_conc <- function() {
  model({
    C <- 0 + err            # hypothetical: mean zero + residual noise
    i <- 0
    while (C < 0) {
      simeps()              # resample all EPSs for this record
      C <- 0 + err
      i <- i + 1
      if (i > 10) break
    }
  })
}
e2 <- et(0:9)

f3 <- rxSolve(mod_conc, e2, sigma = lotri(err ~ 1), seed = 10)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments

# Without resampling some records would be negative; with simeps() all are >= 0
all(f3$C >= 0)
#> [1] TRUE

Preserving original draws when the condition is already met

When the initial draw already satisfies the condition, simeps() is never called and the original value is preserved. This means the marginal distribution of accepted draws is the same as the original distribution restricted to the feasible region, without any hidden bias from extra calls:

# Solve without simeps as reference
mod_ref <- function() {
  model({
    C <- 0 + err
    i <- 0
  })
}

f_ref <- rxSolve(mod_ref, e2, sigma = lotri(err ~ 1), seed = 10)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments

f_resamp <- rxSolve(mod_conc, e2, sigma = lotri(err ~ 1), seed = 10)

# Records where no resampling was needed (i == 0) should match exactly
both <- merge(f_ref, f_resamp[f_resamp$i == 0, ], by = "time")
stopifnot(isTRUE(all.equal(both$C.x, both$C.y)))

Differences from a pure post-hoc truncation

Using simeta() / simeps() inside the model has two advantages over filtering or replacing values after rxSolve():

  1. Self-consistency: downstream model quantities (e.g., PK parameters computed from body weight) always reflect the accepted draw – no risk of a filtered row having inconsistent derived variables.
  2. Parallelism: resampling happens inside the C solver loop, so OpenMP parallelism across subjects is fully preserved.

Next steps