Control for vae (variational autoencoder) estimation method in nlmixr2

Variational-autoencoder NLME estimation (Rohleff et al., CPT:PSP 2025): an LSTM encoder learns the individual posterior q(eta|y) and an rxode2 decoder reconstructs the observations, trained on an ELBO / BICc-ELBO objective for simultaneous population-parameter estimation and covariate selection.

Usage

vaeControl(
  seed = 42L,
  itersBurnIn = 100L,
  klWarmup = 50L,
  gammaIter = 250L,
  iters = 300L,
  nGradStep = 5L,
  hiddenDim = 25L,
  learningRate = 0.005,
  burnInLearningRate = 0.008,
  sigma0 = NULL,
  covariateSelection = TRUE,
  pinCovariates = TRUE,
  muRefCovAlg = TRUE,
  shapes = c("power", "lin", "log", "identity", "center", "hockey"),
  covCenterType = c("median", "mean"),
  covCenter = NULL,
  catCutoff = 0.05,
  covSelectAlpha = 2,
  covSelectSmooth = TRUE,
  gammaSeries = c("reference", "saem"),
  sigma0Interp = c("sd", "reference"),
  residOptimize = c("twoStage", "moment", "optimize"),
  residRhoend = NULL,
  omegaUpdate = c("suffStat", "blend"),
  perNoCor = 0.75,
  inputScale = c("reference", "observed"),
  covSelectMethod = c("auto", "bnb", "l0learn"),
  covSelectMaxExact = 17L,
  bnbStrategy = c("lifo", "fifo", "lc"),
  parEncoderBackward = !isTRUE(getOption("nlmixr2.identical", FALSE)),
  nonMuTheta = c("regress", "grad", "eta", "fix", "none"),
  nonMuEtaOmega = 0.01,
  mStepObjective = c("outer", "elbo"),
  likelihood = c("focei", "foce", "focep", "laplace"),
  objf = c("importanceSampling", "linear"),
  nIsSample = 3000L,
  returnVae = FALSE,
  print = 1L,
  useColor = NULL,
  printNcol = NULL,
  covMethod = c("r,s", "analytic", "r", "s", "sa", "imp", ""),
  optExpression = TRUE,
  sumProd = FALSE,
  literalFix = TRUE,
  literalFixRes = TRUE,
  addProp = c("combined2", "combined1"),
  calcTables = TRUE,
  compress = FALSE,
  adjObf = TRUE,
  ci = 0.95,
  sigdig = 3,
  sigdigTable = NULL,
  rhoend = NULL,
  stickyRecalcN = 4,
  maxOdeRecalc = 5,
  odeRecalcFactor = 10^(0.5),
  outerStickyRecalcN = 4,
  outerMaxOdeRecalc = 5,
  outerOdeRecalcFactor = 10^(0.5),
  indTolRelax = TRUE,
  eventSens = c("jump", "fd"),
  rxControl = NULL,
  ...
)

Arguments

seed

Random seed for the VAE training (encoder init, Adam, reparameterization sampling); default 42. Training is stochastic, so a fixed seed makes every fit reproducible.

itersBurnIn

Number of burn-in iterations (encoder-only, tiny KL weight) before the main EM phase.

klWarmup

Number of KL-annealing iterations over which the KL weight is ramped from a small value to 1 (prevents posterior collapse).

gammaIter

Number of main iterations before the EMA-smoothing phase of the population-parameter update begins.

iters

Total number of main-loop iterations (after burn-in).

nGradStep

Number of Adam gradient steps per EM outer iteration (the reference `L_iter`).

hiddenDim

LSTM hidden dimension (the reference `h_dim`).

learningRate

Adam learning rate used in the main training phase.

burnInLearningRate

Adam learning rate used during burn-in.

sigma0

Encoder prior standard deviation(s) at initialization (a small value giving a sharp initial posterior). `NULL` uses a small default per individual parameter. This is distinct from the `ini()` omega.

covariateSelection

When `TRUE` (default) perform automated BICc-ELBO covariate selection during training; when `FALSE` fit only the covariate structure written in the model. In the `FALSE` case the model-declared covariate coefficients (both linear `beta*WT` effects and transformed ones such as `beta*log(WT/70)`) are estimated in place by the regress M-step regardless of `nonMuTheta`; a `ini(... ~ fix())` coefficient stays fixed.

pinCovariates

When `TRUE` (default) and the model already declares covariate effects, restrict the automatic covariate selection to only the covariate/parameter pairs written in the model – the branch-and-bound search may still drop a declared covariate, but can never add one the model did not specify. A declared covariate that is not a valid search candidate (time-varying, or a raw-linear form that does not match the `log`/centered encoding) is estimated in place by the regress M-step instead, with a note in `$runInfo`. When the model declares no covariates there is nothing to pin and the full search runs. Has no effect when `covariateSelection` is `FALSE`.

muRefCovAlg

When `TRUE` (default) an algebraic/centered covariate effect written in the model (e.g. `wt.cl*(WT/70)` or `wt.cl*log(WT/70)`) is handled as a mu2/mu3 reference: the covariate expression – including its centering – is evaluated into an internal `nlmixrMuDerCov#` data column and the model uses the linear `wt.cl*nlmixrMuDerCov#` form during fitting, so the VAE covariate search never re-centers it. The original expression is restored in the reported model.

shapes

Which parameterizations ("shapes") of a continuous covariate the automatic search may consider, using the same vocabulary as `nlmixr2scm::runSCM()`: `"power"` (`beta*log(COV/ctr)`), `"lin"` (`beta*(COV - ctr)`), `"log"` (`beta*log(COV)`), `"identity"` (`beta*COV`) and `"center"` (`beta*(COV/ctr)`). `"hockey"` is a two-armed piecewise linear relationship knotted at the centering value, written as `beta.low*(COV < ctr)*(COV - ctr) + beta.hi*(COV >= ctr)*(COV - ctr)`; both arms enter or neither does, and it is skipped for a covariate with fewer than `catCutoff` of the subjects on one side of the knot. At most one shape of a covariate may enter a given parameter. Because the selection objective is an ordinary least squares fit with a free intercept, `"power"`/`"log"` span the same model, as do `"lin"`/`"identity"`/`"center"`; the shape therefore decides how an accepted relationship is written back, and when several eligible shapes span the same model the one listed first wins. `"hockey"` spans a strictly larger model than the linear shapes and costs two coefficients rather than one.

May also be a **list**, whose elements are dispatched individually so the two forms mix freely: an element named by covariate (`WT = "power"`) is shorthand for the covariate-wide rule `list(covar = "WT", shapes = "power")`, and a `list(var=, covar=, shapes=)` element restricts one parameter/covariate pair. The most specific rule wins – `var`+`covar` beats `covar`, which beats `var`, which beats a rule naming neither – and ties go to the rule listed last.

In the list form, **naming a covariate also puts it in the search**. `fixCov = TRUE` (the default, given as an element of the list) fixes the searched set to exactly the covariates named, so `shapes = list(WT = "power")` searches `WT` and nothing else. Add `fixCov = FALSE` to restrict parameterizations without restricting the search, which is what the list form meant previously. A shape value of `TRUE` means "eligible, default shapes", and is how a categorical covariate is named (`list(WT = "power", SEX = TRUE)`) since a categorical takes no parameterization. A `var`-only rule makes every covariate eligible on that parameter alone; a rule naming neither `var` nor `covar` contradicts `fixCov = TRUE` and is an error. A character vector names no covariate, so `fixCov` does not apply and every covariate stays searchable.

`fixCov` is ignored when the model itself declares covariate effects: that already restricts the search (see `pinCovariates`) and the declaration is the more specific statement. The disagreement is reported in `$runInfo`, as is every covariate `fixCov` excludes. Categorical covariates always enter as indicators and take no shape, but `fixCov` still governs whether they are searched at all.

covCenterType

Statistic used to center a continuous covariate, `"median"` (default) or `"mean"`, computed over subjects rather than rows.

covCenter

Named numeric vector of centering values overriding `covCenterType` for those covariates, e.g. `c(WT = 70)`. Names are matched case-insensitively.

catCutoff

Minimum proportion of subjects a non-reference level must hold to get its own indicator. Rarer levels are lumped with the reference. Default `0.05`; `0` tests every level.

covSelectAlpha

Starting multiplier for the covariate-selection L0 penalty, ramped linearly from `covSelectAlpha` down to `1` over the `klWarmup` warmup iterations and held at `1` afterward (matching the reference implementation's `linspace(alpha, 1, kl_iter)`). Values `> 1` penalize covariate entry more heavily early in training; `1` disables the ramp.

covSelectSmooth

When `TRUE` (default) the covariate selection regresses the SAEM sufficient statistic – an exponential moving average of the posterior means, updated with the same gain as the M-step – rather than the current posterior means. This matches the reference implementation (Rohleff et al. 2025), which is the reason for the default. In practice it changes little: `gamma` is exactly 1 until `gammaIter`, so the statistic equals the posterior mean for most of a run and is averaged only over the closing tail. `FALSE` regresses the current posterior means.

gammaSeries

Decaying step-size series used once the smoothing phase starts (after `gammaIter`); the gain is 1 throughout the EM phase either way.

* `"reference"` (default): `1/(iter - gammaIter)`, the textbook Kuhn-Lavielle series the reference implementation uses. The first smoothing step is still a full replacement, and the decay follows. * `"saem"`: `1/(1 + iter - gammaIter)`, the CONTINUATION form saemControl() uses – nlmixr2est's SAEM builds its series so it continues rather than repeating a gain of 1, so the decay begins at `1/2`. Select this to match the step-size convention of the other nlmixr2 estimation methods rather than the reference.

sigma0Interp

How `sigma0` is turned into the encoder's initial posterior spread. The encoder head emits `logSigma` and forms `diag(L) = exp(logSigma)`, so `diag(L)` is the posterior standard deviation.

* `"sd"` (default): the bias is `log(sigma0)`, so the initial posterior SD is `sigma0` – what the argument says it is. * `"reference"`: the bias is `log(sigma0^2)`, matching the reference implementation, whose initial posterior SD is therefore `sigma0` SQUARED (`1e-6` rather than `1e-3` for the first neonatal dimension). The reference documents `sigma0` as a standard deviation, so this appears to be unintended there; it is offered only to reproduce its published behavior.

residOptimize

How the residual-error parameters are estimated.

Residual forms the optimizer estimates: `add`, `prop`, `add + prop`, `pow`, `lnorm`, and a `boxCox` or `yeoJohnson` lambda (bounded to `(-2, 2)`). For a transform-both-sides model the objective transforms `dv` only and carries the log-Jacobian, since `f` leaves the solve already on the transformed scale.

`nonMuTheta = "grad"` bypasses this entirely: the analytic outer gradient already carries a residual sigma and a transform lambda as its own directions, so those parameters are stepped by the gradient through Adam and the two-stage path never runs. Which converges better is model-dependent.

* `"moment"`: the closed-form moment estimator. For a model with a single additive error this is exactly the optimum (`sqrt(SSE/n)`); for any other error model it is either a different estimator or, for the forms with no closed form (`pow`, Box-Cox, Yeo-Johnson), no estimator at all – the parameter stays at its `ini()` value. There is no moment estimator for a log-likelihood (`ll()`) parameter either, so those also stay at `ini()`; use `"twoStage"` for such a model. * `"twoStage"` (default): block coordinate descent, as `npag`'s `residOptimize = "alternate"` does. Stage one optimizes the non-mu-referenced structural thetas with the residual parameters held, so it is driven by `(dv - f)`; stage two then holds those and optimizes the residual parameters alone against the extended least-squares objective `sum[(y-f)^2/r + log r]` over the CACHED `(y, f)` pairs. Because `f` is fixed by stage one, stage two needs no ODE re-solve – the same structure SAEM uses. On `theo_sd` this beats the moment estimator on both a pure-additive model (objective 131.79 vs 131.81) and a combined one (121.03 vs 122.47).

Which parameters stage two owns is decided per parameter: an error parameter, or one that no `d/dt()` right-hand side, initial condition or dosing modifier can reach. The second case is what a log-likelihood (`ll()`) or generalized endpoint needs – its residual-like parameters are plain thetas with no error row, and on the error-only rule stage two was empty for such a model, silently making `"twoStage"` behave like `"optimize"`. When no regressed theta qualifies (every one feeds the solve) stage two has nothing to do and `residOptimize` has no effect. * `"optimize"` (EXPERIMENTAL, diagnostic): a single JOINT solve over the structural and residual parameters together, against the full outer objective. Fine with one free residual parameter, but with `add` and `prop` both free it diverges – they are near-collinear, and routing the residual through the full outer objective lets the Laplace terms move with it at frozen etas (objective 320.7 against the moment estimator's 122.5). Retained for comparison; prefer `"twoStage"`.

residRhoend

Final trust-region radius (`rhoend`) of the bounded `bobyqa` that estimates the residual parameters – its convergence tolerance. `NULL` (default) derives it from `sigdig` (`10^(-sigdig)`), the same way every other optimizer tolerance in the package is derived, so `sigdig` stays the single knob that moves them together. Set it explicitly when the residual step should converge tighter than the rest: it runs with the ODE frozen, so tightening it is far cheaper than tightening `rhoend`, which also tightens the structural regression that re-solves per candidate.

omegaUpdate

How the population variances are updated in the covariate M-step. `"suffStat"` (default) follows the reference: `omega` is formed from the EMA sufficient statistics and ASSIGNED outright. `"blend"` is the historic behavior, blending the freshly computed `omega` with the previous value at the M-step gain (so it is smoothed twice). Applies to `omega` only.

Note this option reaches only ONE of the two omega M-steps. Which one runs is decided by `covariateSelection`: with `TRUE` the covariate M-step runs and honors `omegaUpdate`; with `FALSE` the plain closed-form M-step runs, whose variances are always raw posterior moments blended at the gain. A declared correlated block's OFF-diagonals always follow whichever estimator that branch's diagonal used – estimating the two halves of one block by different estimators need not even give a positive-definite result.

The two settings are the SAME update while the gain is 1, which it is throughout burn-in and the EM phase (assigning a value and blending it in with weight 1 are the same operation); they differ only once `gammaIter` decays the gain. A short run at default settings will show no difference.

`mStepObjective` does not enter the omega update at all – it scores the non-mu theta M-step. Omega has a closed-form EM update from the variational posterior either way.

The residual error estimate is still EMA-smoothed on the standard-deviation scale, where the reference smooths the residual sum of squares and takes the root afterwards – a known remaining difference. Matching it would need per-endpoint sufficient statistics plus an optimizer branch for the error models with no closed form (`add + prop`, `add + pow`, Box-Cox / Yeo-Johnson), as saemControl() does.

perNoCor

Fraction of the EM phase (`gammaIter` iterations) over which a declared correlated `omega` block is held at zero correlation, letting the variances settle before the correlations are estimated. This is saemControl()'s `perNoCor` rule (0.75 there as well); it has no effect on a model with no declared off-diagonals.

Held at ZERO, following saem, not at the `ini()` value: retaining an initial covariance while the variances shrink around it can leave the block non-positive-definite. A `fixed()` covariance is exempt – it is not being estimated, so it keeps its value through the hold and out the other side.

The fraction is of the EM phase, `min(gammaIter, iters)`, not of the whole run. That matters: the gain is 1 for `it <= gammaIter`, so the release point falls while the gain is still 1 and the correlations are estimable the moment they are unfrozen. (This is why no gain restart is needed here, whereas emviControl() – whose run has no separate unit-gain phase – has to restart the off-diagonal gain at release.)

A value greater than 1 is an ABSOLUTE iteration count rather than a fraction, and must be a whole number. Prefer the absolute form whenever a run may be resumed or reproduced at a different length: a fraction of a shorter run is a different schedule, not the same one truncated.

inputScale

Which observations the encoder-input centering and scaling are computed over. `"reference"` (default) matches the reference implementation, which takes the mean and SD across the whole padded observation matrix, so the zero padding of subjects with fewer observations enters both statistics. On a ragged dataset that is a materially different scale from `"observed"`, which uses only the observed values (on the neonatal case study the SD is 1582 against 506). Affects only the encoder's inputs, never the likelihood.

covSelectMethod

How the covariate M-step searches subsets. `"bnb"` is the exact branch-and-bound; it becomes impractical past a few dozen candidate covariates. `"l0learn"` has the `L0Learn` package propose supports, which are then scored and polished with the same exact objective – so the search is approximate but the scoring is not. `"auto"` (default) uses `"l0learn"` for a latent dimension with at least `covSelectMaxExact` candidate covariates and `"bnb"` otherwise. Set `covSelectMaxExact = Inf` to force the exact search everywhere.

covSelectMaxExact

Search size at or above which `covSelectMethod = "auto"` switches a latent dimension to `L0Learn` (default `17`, just above the measured wall-clock crossover of roughly 16 bits – see `tools/benchVaeCovSelect.R`, which finds the same crossover in bits whether a covariate carries one shape or two). Measured in bits of feasible-support space – `sum over covariates of log2(1 + shapes tried)` – after `pinCovariates` trimming, so it is the size of the search actually run. One shape per covariate costs exactly 1 bit, so with `shapes` set to a single shape this is a plain candidate count; two shape families of one covariate cost `log2(3)`, keeping the exact search's worst-case node budget the same either way. `Inf` forces the exact branch-and-bound everywhere.

bnbStrategy

Frontier discipline for the exact branch-and-bound covariate selection: `"lifo"` (default, last-in-first-out depth-first search), `"fifo"` (first-in-first-out) or `"lc"` (least cost / best-first). The solver is exact, so the selected covariates are identical for every strategy; only the search order (and thus efficiency) differs.

parEncoderBackward

Parallelize the encoder backward (gradient) pass over subjects. Defaults to `TRUE` unless `options(nlmixr2.identical = TRUE)` is set (which flips the default to `FALSE`); an explicit value here always wins. The encoder forward pass and the covariate branch-and-bound already run multi-threaded and are bit-identical to the serial run. The backward gradient is a continuous cross-subject sum, so parallelizing it (per-thread partials reduced in thread order) makes the result deterministic for a fixed number of `cores` but no longer bit-identical to the serial path: the per-step gradient differs at ~1e-12, which compounds through the iterative SGD/EM training to a small final difference (well below any estimation tolerance), and results may differ across different `cores`. When it is active (and `cores > 1`) a note is added to the fit's `$runInfo`. Set this to `FALSE` – or globally `options(nlmixr2.identical = TRUE)` – for bit-identical, fully reproducible results.

nonMuTheta

How to treat a structural population `theta` that has no random effect (is not mu-referenced) so it can still be estimated by the VAE (which only estimates parameters that occupy the latent space). For the eta-injection modes a small eta is injected so the parameter enters the latent space, and the reported fixed effect is `theta + mean(eta)` with the temporary eta dropped from the output model.

* `"regress"` (default, matching `saemControl(nonMuTheta=)`): no eta is injected; instead each such theta is estimated directly, re-optimized every M-step by a bounded `bobyqa` regression against the full FOCEi outer objective (bounds from the `ini()` lower/upper), blended with the M-step gain. `mStepObjective` selects which objective that regression targets. This recovers a no-random-effect population parameter without adding a spurious random effect. `nonMuEtaOmega` is unused in this mode. * `"grad"`: same target as `"regress"` but stepped with the EXACT analytic outer gradient (Almquist sensitivity equations, the machinery behind `foceiControl(fast=TRUE)`) instead of a derivative-free search: one augmented sensitivity solve per M-step replaces the bobyqa sweep. Both modes optimize the same full outer objective (with every mu-referenced theta held at its current M-step value), so this changes the optimizer, not the target. It is also the more natural fit for the method: the gradient is handed to the SAME Adam machinery that moves the encoder weights, so the parameter is learned alongside the rest of the model on a shared schedule (same gain, same KL warmup gate), whereas `"regress"` pauses each M-step to run a separate derivative-free optimizer to convergence and adopts its answer. This is NOT a speed option – it is measurably SLOWER than `"regress"` (on `theo_sd`, 1.47x with one non-mu theta and 1.13x with three; the gap narrows as the number grows, since bobyqa's cost scales in it and a single solve does not, but it does not close). Choose it for accuracy: the exact gradient lands closer to the maximum-likelihood value than the derivative-free search (`theo_sd` non-mu `tv`: 3.4294 vs 3.4324, against a FOCEi MLE of 3.4293). Applies to a conditionally Gaussian model and to a single non-Gaussian (`ll()`/generalized) endpoint, which differentiates the log-density directly. Falls back to `"regress"` when the model is out of analytic scope (`linCmt()`, IOV, `fo`, a multi-endpoint or censored `ll()` model, ...); `nonMuEtaOmega` is unused. * `"eta"`: inject the eta with an ESTIMATED omega (starting at `nonMuEtaOmega`); the typical value is estimated and appears in the iteration table. * `"fix"`: inject the eta with omega held FIXED at `nonMuEtaOmega` AND hold the typical-value theta fixed at its `ini()` value. Nothing about the parameter is estimated, so it is not shown in the iteration table (it is reported at its `ini()` value, marked fixed, with the injected eta dropped). * `"none"`: leave non-mu-referenced thetas frozen at their `ini()` value (the historic behavior).

nonMuEtaOmega

Variance of the eta injected for a non-mu-referenced theta (starting value for `nonMuTheta="eta"`, fixed value for `nonMuTheta="fix"`; unused for `"regress"`).

mStepObjective

Objective the non-mu-referenced theta M-step (`nonMuTheta = "regress"` or `"grad"`) is optimized against. It has no effect when there is no non-mu-referenced structural theta, and it never changes the encoder/ELBO training step or the covariate branch-and-bound criterion, both of which always follow the reference.

* `"outer"` (default): the full FOCEi outer objective – the frozen-eta joint likelihood PLUS the Laplace determinant, `0.5*log|Omega^-1|` and the DV-transform Jacobian. This is a deliberate deviation from Rohleff et al. (2025): it keeps the quantity being optimized equal to the objective the fit reports, and it is the functional the analytic outer gradient differentiates, so `nonMuTheta = "grad"` optimizes one target rather than stepping one and scoring another. * `"elbo"`: the reference behavior – the plain variational bound (frozen-eta joint likelihood, no Laplace term), matching the M-step in Rohleff et al. (2025). Use it to reproduce the reference implementation. The analytic outer gradient does not apply to this objective, so `nonMuTheta = "grad"` is downgraded to `"regress"` with a note in `$runInfo`.

The two objectives differ by terms that depend on the non-mu thetas through the eta Hessian, so they can land on different estimates, and – because those estimates feed the latent means the covariate search regresses on – on different covariate sets.

likelihood

Inner likelihood used for the objective, EBEs, and gradients, all run through the same FOCEi inner interface: `"focei"` (default, with eta-epsilon interaction), `"foce"` (no interaction, NONMEM FOCE with R frozen at the population prediction), `"focep"` (FOCE+, no interaction but R evaluated at the live conditional eta), or `"laplace"`.

objf

Which objective-function value is active for AIC/BIC/BICc. Both the linearization and importance-sampling -2LL are always computed and stored; this selects the default active one.

nIsSample

Number of importance-sampling draws for the IS -2LL.

returnVae

When `TRUE` return the raw VAE training object instead of the nlmixr2 fit.

print

Either a scalar print-frequency (`0` = suppress, `1` (default) = every evaluation, `N` = every Nth), OR a pre-built [iterPrintControl()] object. Equivalent to `iterPrintControl(every = print, ncol = printNcol, useColor = useColor)`.

useColor

Logical (or `NULL`) emit ANSI bold/color escapes in the iteration print. `NULL` (default) defers to [crayon::has_color()].

printNcol

Integer (or `NULL`) parameter columns per row before wrapping. `NULL` (default) uses `floor((getOption("width") - 23) / 12)`.

covMethod

Method for calculating the covariance at the VAE estimates, run through the FOCEi covariance step; the same choices as foceiControl(): "analytic" (default), "r,s", "r", "s", or "" to skip.

optExpression

Optimize the rxode2 expression to speed up calculation. By default this is turned on.

sumProd

Is a boolean indicating if the model should change multiplication to high precision multiplication and sums to high precision sums using the PreciseSums package. By default this is FALSE.

literalFix

boolean, substitute fixed population values as literals and re-adjust ui and parameter estimates after optimization; Default is `TRUE`.

literalFixRes

boolean, substitute fixed population values as literals and re-adjust ui and parameter estimates after optimization; Default is `TRUE`.

addProp

Type of additive-plus-proportional error: `"combined1"`, where standard deviations add: $$y = f + (a + b\times f^c) \times \varepsilon$$; or `"combined2"`, where variances add: $$y = f + \sqrt{a^2 + b^2\times f^{2\times c}} \times \varepsilon$$. Here y = observed, f = predicted, a = additive sd, b = proportional/power sd, c = power exponent (1 in the proportional case).

calcTables

This boolean is to determine if the foceiFit will calculate tables. By default this is TRUE

compress

Should the object have compressed items

adjObf

is a boolean to indicate if the objective function should be adjusted to be closer to NONMEM's default objective function. By default this is TRUE

ci

Confidence level for some tables. By default this is 0.95 or 95% confidence.

sigdig

Specifies the "significant digits" that the ODE solving requests. This is NULL by default, and while it is NULL it has no effect at all: rxSolve() uses the standard atol/rtol (and the standard sensitivity and steady-state tolerances). sigdig only changes a tolerance when you ask for it explicitly.

When it is supplied, the tolerances are derived with one solver-independent formula – the same for stiff, non-stiff and auto-switching solvers. The rtol exponent IS sigdig and atol sits three orders below it:

• rtol = 10^(-sigdig), atol = 10^(-sigdig-3)

• the sensitivity tolerances match the main solve, so rtolSens = rtol and atolSens = atol (gradients and covariances are built from them)

• the steady-state tolerances run one order looser than the corresponding main tolerance, so ssRtol = ssRtolSens = 10*rtol and ssAtol = ssAtolSens = 10*atol

Each of these is set only when you did not pass that tolerance yourself; a tolerance you supply always wins. Because they are resolved independently, an explicit atol/rtol overrides the main solve but does not propagate to the sensitivity or steady-state tolerances – set those directly if you need them changed too.

This mapping matches how nlmixr2est derives solver tolerances from its optimization sigdig, so a sigdig used for estimation and the same sigdig used for a plain rxSolve() mean the same thing. Note it is keyed to sigdig as a request for that many significant digits, and is looser than the atol/rtol defaults for small sigdig – at sigdig = 4 it gives rtol = 1e-4 against a default rtol = 1e-6. Raise sigdig, or set atol/rtol directly, when you want a tighter solve.

sigdigTable

Significant digits in the final output table. If not specified, then it matches the significant digits in the `sigdig` optimization algorithm. If `sigdig` is NULL, use 3.

rhoend

Final trust-region radius (`rhoend`) of the inner bounded `bobyqa` used by the non-mu / covariate regress M-step. `NULL` (default) derives it from `sigdig` (`10^(-sigdig)`, matching the optimizer convergence tolerance), or `1e-4` when `sigdig` is `NULL`.

stickyRecalcN

The number of bad ODE solves before reducing the atol/rtol for the rest of the problem.

maxOdeRecalc

Maximum number of times to reduce the ODE tolerances and try to resolve the system if there was a bad ODE solve.

odeRecalcFactor

The ODE recalculation factor when ODE solving goes bad, this is the factor the rtol/atol is reduced

outerStickyRecalcN

The number of bad analytic outer solves for a subject before its loosened tolerance is kept for the rest of the problem; the outer counterpart of `stickyRecalcN`.

outerMaxOdeRecalc

Maximum number of times to reduce the ODE tolerances for a single subject and retry when the analytic outer (augmented sensitivity) solve fails. Tracked separately from `maxOdeRecalc`, which governs the inner problem. A subject that solves after loosening still contributes an analytic gradient instead of dropping the whole gradient to finite differences.

outerOdeRecalcFactor

The factor the atol/rtol is loosened by on each analytic outer retry; the outer counterpart of `odeRecalcFactor`.

indTolRelax

When `TRUE` (default), only subjects whose ODE solve produced NaN/Inf have their tolerances relaxed, and the relaxed tolerance persists across optimizer calls (sticky). When `FALSE`, all subjects have their tolerances relaxed on each retry and tolerances are reset afterward.

eventSens

Controls how dosing/event-parameter (`alag`, `F`, `rate`, `dur`) sensitivities are computed for THETA/ETA gradients: `"jump"` (default) uses rxode2's analytic event sensitivities; `"fd"` uses the legacy finite-difference behavior.

rxControl

`rxode2` ODE solving options during fitting, created with `rxControl()`

...

Other arguments to control SAEM.

Value

vae control structure (class `vaeControl`)

Details

Covariate selection – MIQP vs. branch-and-bound. Per latent parameter the selection step minimizes the same L0/BIC objective `RSS_S/omega + log(N)*|S|` over subsets `S` of the candidate covariates (`RSS_S` is the residual sum of squares of the ordinary-least-squares fit on the intercept plus `S`). The reference implementation (Rohleff et al.) writes this as a Mixed-Integer Quadratic Program (MIQP) – binary include/exclude indicators with big-M constraints – and solves it with the commercial Gurobi solver through `cvxpy`. No MIQP-capable solver is freely available in R: Gurobi is commercial/licensed, and the open QP solvers on CRAN (e.g. `osqp`) are continuous-only and cannot represent the binary selection. A continuous convex relaxation (L1 / lasso) would be solvable but only approximates best subset.

This package instead solves the identical L0/BIC objective EXACTLY with a self-contained branch-and-bound: each candidate support's coefficients are the closed-form OLS fit and branches are pruned by a valid lower bound (the RSS of the OLS fit using all still-free covariates). It therefore returns the same optimum the MIQP would – no commercial dependency and no relaxation/accuracy loss – and scales to a few dozen covariates. The search is worst-case exponential in the number of covariates, but the pruning makes the practical (sparse) case fast (e.g. 32 candidate covariates in a fraction of a second).

Author

Matthew L. Fidler