Skip to contents

Control for the variational inference methods emvi and fbvi in nlmixr2

Source: R/emviControl.R
emviControl.Rd

Variational-inference NLME estimation in the style of Kucukelbir et al. (2017): the latent variables are transformed to an unconstrained real coordinate space, a Gaussian variational family is posited there, and the ELBO is maximized by stochastic gradient ascent using the reparameterization trick.

Usage

emviControl(
  seed = 42L,
  iters = 300L,
  nMc = 1L,
  viFamily = c("fullRank", "meanField"),
  pointEstimate = NULL,
  optim = c("advi", "adam"),
  adaptEta = TRUE,
  perNoCor = 0.75,
  etaCandidates = c(0.01, 0.025, 0.05, 0.1, 0.25),
  tau = 1,
  alpha = 0.1,
  tol = NULL,
  evalElbo = 100L,
  klWarmup = 0L,
  temperInit = 10,
  likelihood = c("focei", "foce", "focep", "laplace"),
  returnVi = FALSE,
  resume = NULL,
  print = 1L,
  useColor = NULL,
  printNcol = NULL,
  covMethod = c("vi", "analytic", "r,s", "r", "s", ""),
  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,
  stickyRecalcN = 4,
  maxOdeRecalc = 5,
  odeRecalcFactor = 10^(0.5),
  indTolRelax = TRUE,
  eventSens = c("jump", "fd"),
  rxControl = NULL,
  ...
)

Arguments

seed

Random seed for the variational optimization (reparameterization sampling); default 42. The Monte-Carlo gradient is stochastic, so a fixed seed makes every fit reproducible. Reparameterization noise is drawn from a counter-based stream keyed by the global iteration index, so a shorter run is a bit-for-bit prefix of a longer one and results are independent of the number of cores.

iters

Total number of stochastic gradient ascent iterations.

nMc

Number of Monte-Carlo samples used to approximate the ELBO gradient at each iteration (the paper's `M`; typically 1-10).

viFamily

Variational family in the unconstrained space. `"fullRank"` (default) uses a block full-rank Gaussian: a dense `neta x neta` Cholesky factor per subject plus a dense block over the population vector (mean-field across blocks). `"meanField"` uses a fully factorized (diagonal) Gaussian. Mean-field is faster but is known to underestimate marginal variances.

pointEstimate

Which of the two methods to run, normally left at its `NULL` default so it follows `est`: `est="emvi"` implies `TRUE` and `est="fbvi"` implies `FALSE`. `est` wins over a contradicting value, and says so – it has to, because re-estimating a fit with the other method pipes the completed fit's control forward. `TRUE` runs the variational-EM hybrid: the variational posterior covers the per-subject etas only, and the population parameters (thetas / omega / residual error) are point estimates maximized by the ELBO gradient; output semantics match FOCEi/SAEM. `FALSE` runs full Bayes: the variational posterior also covers the unconstrained population vector, with flat priors.

Two things about "flat" are worth being explicit about, because they define the prior rather than merely describe the implementation. (1) A BOUNDED theta is fitted on its unconstrained scale, and the log-determinant of that constraining transform IS added to the full-Bayes objective, so the flat prior is flat on the NATURAL parameter, as in Stan. It is deliberately NOT added when `pointEstimate=TRUE`: a maximum-likelihood estimate has to stay invariant to reparameterization, which is why Stan's own `optimize` defaults to `jacobian=0`. (2) The between-subject variances are carried as per-eta LOG-variances and no Jacobian is applied to them, so the prior is flat on `log(omega)` – the conventional weakly-informative choice for a scale parameter, but a choice, not an accident: it is not flat on `omega`.

These point estimates maximize the ELBO, NOT the likelihood, and the difference is not merely cosmetic. Since `ELBO = log p(y|theta) - KL(q || p(eta|y,theta))`, any dependence of that KL on `theta` displaces the maximizer from the MLE – "variational maximum likelihood". For VARIANCE components the displacement has a known direction: a variational family that understates posterior spread makes the omega M-step, `Omega = mean_i(mu_i mu_i' + Sigma_i)`, inherit that understatement, so between-subject variability is biased DOWNWARD. The bias is worst for `viFamily="meanField"`, which cannot represent within-subject posterior correlation at all; `"fullRank"` can, which is why it is the default. Structural (typical-value) parameters are far less affected. If the between-subject variances are themselves the quantity of interest, prefer `"fullRank"` and cross-check against `est="focei"` or `est="saem"`.

optim

Stochastic optimizer. `"advi"` (default) uses the adaptive step-size sequence from the ADVI paper (Eqs 10-11) – the value keeps that name because the step-size rule really is the published one, even though the surrounding method is not; `"adam"` uses Adam.

adaptEta

When `TRUE` (default) adaptively choose the step-size scale `eta` by a short search over `etaCandidates` before the main loop; when `FALSE` use a fixed `eta` (the first `etaCandidates` entry).

perNoCor

Fraction of the run over which a declared correlated `omega` block is held at zero correlation, letting the population 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. When the correlations are released the off-diagonal gain restarts rather than continuing the decayed schedule, so they are still estimable at the point they are unfrozen. This is needed here and not in vaeControl(): there the fraction is of the EM phase, where the gain is still 1 at release, whereas this run has no such phase and the fraction lands in the decayed schedule.

A value greater than 1 is an ABSOLUTE iteration count rather than a fraction, and must be a whole number. The resolved absolute iteration is stored with the fit and reused by `resume=`, so a resumed run releases the correlations at the same global iteration a single long run would – recomputing the fraction from the resumed call's `iters` would re-apply a hold the original run had passed. Only `est="emvi"` estimates the block at all; `est="fbvi"` carries `omega` as per-eta log-variances and errors on a correlated model rather than dropping the off-diagonal.

etaCandidates

Candidate step-size scales searched when `adaptEta` is `TRUE`. The default is narrower and smaller-valued than the paper's `c(0.01, 0.1, 1, 10, 100)` because these gradients come from FOCEi sensitivities on the model's own scale rather than from AD through a Stan program, so the useful step sizes sit lower; the paper's grid can be passed verbatim if wanted. Each candidate costs `min(iters, 75)` iterations (a diverging one aborts early and is cheap), so the search is a substantial share of a fit – widen it deliberately. When the search selects the largest or smallest candidate the grid may be the binding constraint, and the fit says so in `$runInfo`; `$etaScores` reports the per-candidate scores behind the choice.

tau

Stabilizing constant `tau > 0` in the step-size denominator (paper Eq 10); the step-size is insensitive to it.

alpha

Weighting `alpha` in (0, 1) of new vs old gradient information in the step-size memory recursion (paper Eq 11).

tol

Convergence tolerance on the relative change in the ELBO: the loop stops early once the change falls below this. `NULL` (default) derives it from `sigdig` as `10^(-sigdig)`, the same rule `saemControl()` and `foceiControl()` use for their optimizer tolerances, so it tightens with `sigdig` instead of staying pinned. `0` disables early stopping (run all `iters`). Because the per-iteration ELBO is an `nMc`-sample Monte-Carlo estimate and therefore noisy, the test compares the MEAN over the last `evalElbo` iterations against the mean over the window before it, rather than consecutive iterations. The `adaptEta` step-size search never stops early on this criterion (its scorer reads a short run as divergence).

evalElbo

Window length, in iterations, for the `tol` convergence test. Stan's ADVI re-evaluates the ELBO every 100 iterations with fresh draws; averaging the draws already taken is the cheaper equivalent. Shrunk automatically on a short run so `iters` well below `evalElbo` can still trigger the check.

klWarmup

Number of iterations of PRIOR TEMPERING (0, the default, disables it). During the warm-up the population prior is inflated by a factor ramping geometrically from `temperInit` down to 1, which down-weights the prior term of the ELBO and keeps the per-subject variational posterior from collapsing before it is informative. This is the variational analogue of saemControl()'s `perSa` simulated-annealing phase (true simulated annealing does not transfer: variational inference has no MCMC kernel to keep wide).

It CHANGES THE OBJECTIVE FUNCTION MID-RUN. Early iterations maximize a tempered surrogate rather than the ELBO, so the convergence theory for the ELBO does not cover the warm-up, the reported ELBO trace is not comparable across the boundary, and the `tol` early-stopping test is suppressed until tempering ends. The `adaptEta` step-size search also scores candidates on the untempered objective. Off by default for those reasons; turn it on for a model where the variational scale collapses early.

temperInit

Initial prior inflation factor for `klWarmup` tempering (default 10); ignored when `klWarmup = 0`.

likelihood

Inner likelihood used for the per-subject objective and gradient, run through the FOCEi inner interface: `"focei"` (default), `"foce"`, `"focep"`, or `"laplace"`.

returnVi

When `TRUE` return the raw variational optimization object instead of the nlmixr2 fit.

resume

Optional warm-resume state: a previous `emvi`/`fbvi` fit (or its `$env$viState`). The optimization continues from that state for `iters` more iterations, bit-for-bit identical to a single fresh run of the combined length (the counter-based RNG is keyed by the global iteration index).

That equivalence requires every schedule point to be an ABSOLUTE iteration. A FRACTIONAL `perNoCor` cannot provide it, and not because of any bookkeeping that could be fixed: `perNoCor = 0.75` of one 120-iteration run releases the correlations at iteration 90, while 0.75 of a first 60-iteration leg releases at 45. Those are different schedules, and the resulting correlation estimates genuinely differ. Pin the schedule (`perNoCor = 90`) whenever a fit may be resumed; the resolved value is then stored with the fit and reused.

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 covariance. In this discussion, R is the Hessian matrix of the objective function. The S matrix is the sum of each individual's gradient cross-product (evaluated at the individual empirical Bayes estimates).

"sa" (default) Use the stochastic-approximation Fisher Information Matrix. After estimation, a dedicated covariance phase (nSaCov iterations) holds the parameters at the converged estimate and keeps resimulating the individual parameters, Monte-Carlo averaging the Louis observed-information integrand into a converged FIM decoupled from the cooling schedule (the approach used by Monolix; Kuhn & Lavielle 2005). Always includes every estimated population parameter (theta, the Omega diagonal variances, and residual).

"analytic" Compute the FOCEI analytic observed-information covariance at the converged SAEM estimates. When the model is out of analytic-covariance scope (e.g. linCmt(), a non-normal likelihood, or a non-SD IOV parameterization) or the result is not positive definite, it falls back to the linearized Fisher information (linFim) with a message.

"linFim" Use the Linearized Fisher Information Matrix to calculate the covariance.

"fim" Use the Fisher Information Matrix accumulated during SAEM estimation to calculate the covariance. Like sa it inverts the observed information to a full theta + Omega diagonal + residual covariance, but uses the (noisier) estimation-phase matrix rather than a dedicated cov phase.

For both fim and sa the simulation-based Fisher information covers the structural theta, the Omega diagonal variances, and additive residual error. Off-diagonal Omega covariances and proportional/combined residual error are not estimated reliably by the simulation FIM (the complete-data correction is unstable when between-subject variability dominates the residual), so those variance-block standard errors are spliced in from the linearized FIM (linFim).

"r,s" Uses the sandwich matrix to calculate the covariance, that is: \(R^-1 \times S \times R^-1\)

"r" Uses the Hessian matrix to calculate the covariance as \(2\times R^-1\)

"s" Uses the crossproduct matrix to calculate the covariance as \(4\times S^-1\)

"" Does not calculate the covariance step.

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.

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

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

variational-inference control structure (class `emviControl`)

Details

One control serves both methods, which differ only in what the variational posterior covers:

  • `est="emvi"` – variational EM. The variational family covers the per-subject etas only; the population parameters (thetas, omega, residual error) are point estimates updated by an M-step, so the output semantics match FOCEi/SAEM.

  • `est="fbvi"` – full-Bayes variational inference. The variational posterior additionally covers the unconstrained population vector under flat priors, which is the closer analogue of the published algorithm.

`pointEstimate` is the switch between them and defaults to whichever the chosen `est` implies. `est` decides: a control that contradicts it is overridden, with a message saying so.

Neither method is the published ADVI algorithm, which is why neither is named for it. Two deviations matter. First, even `fbvi` parameterizes the between-subject omega as per-eta log-variances rather than freely, so a full-rank family does not yield the joint posterior covariance over the whole parameter vector that published ADVI gives. Second, the gradient of the log-joint comes from the FOCEi forward sensitivities (inner per-subject eta gradient plus the outer population sensitivity contraction) rather than from automatic differentiation – there is no AD in this implementation, which is why "advi" would misname it. Results should therefore not be read as reproducing a Stan `vb()` fit except on models where the two objectives coincide. The whole optimization loop runs in C++.

Author

Matthew L. Fidler