SAEM

Maturity: beta — see Feature Maturity for what this means.

Stochastic Approximation Expectation-Maximization (SAEM) is an alternative estimation method that uses MCMC sampling for random effects instead of MAP optimization. It is more robust to local minima and can handle complex random effect structures, including models with inter-occasion variability (IOV).

Algorithm Overview

SAEM replaces the deterministic inner loop of FOCE with stochastic sampling, following the Monolix convention with a two-phase step-size schedule.

References

  • Delyon, Lavielle, Moulines (1999). Convergence of a stochastic approximation version of the EM algorithm. Annals of Statistics, 94–128.
  • Kuhn & Lavielle (2004). Coupling a stochastic approximation version of EM with an MCMC procedure. ESAIM: Probability and Statistics 8:115–131.

Two-Phase Schedule

Phase 1: Exploration (iterations 1 to K1)

Step size \(\gamma_k = 1\). The algorithm explores the parameter space rapidly, with the sufficient statistics being fully replaced each iteration. This allows fast movement toward the basin of the MLE.

Default: 150 iterations.

Phase 2: Convergence (iterations K1+1 to K1+K2)

Step size \(\gamma_k = 1/(k - K_1)\). The algorithm performs a decreasing-weight average, which guarantees almost-sure convergence to the MLE under regularity conditions.

Default: 250 iterations.

Per-Iteration Steps

Each SAEM iteration consists of:

1. E-Step: Sampling

For each subject, sample from the conditional distribution of random effects:

\[ p(\eta_i | y_i, \theta, \Omega, \sigma) \]

Two samplers are available:

Metropolis-Hastings (default, n_leapfrog = 0)

Two MH kernels run per subject per SAEM iteration (the mixture of Kuhn & Lavielle 2004):

Kernel 1 — block proposal: run n_mh_steps symmetric random-walk steps \(\eta_{\text{prop}} = \eta_{\text{current}} + \delta_i \cdot L \cdot z\), with \(z \sim N(0, I)\) and \(L = \text{chol}(\Omega)\). This mixes the joint scale efficiently. For FREM models the covariate ETAs are near-deterministic (their pseudo-observations pin them to a posterior SD of \(\approx\sqrt{\text{EPSCOV}} \ll \sqrt{\Omega_{jj}}\)), so a full-scale joint proposal for them is rejected every time and the joint acceptance collapses to 0%. Each covariate coordinate is therefore damped in the block move by a factor \(\min(1, \sqrt{\text{EPSCOV}}/\sqrt{\Omega_{jj}})\), so the joint step still explores the correlated PK block while the componentwise sweep handles the covariate coordinates themselves. Ordinary (non-FREM) ETAs use a damping factor of exactly 1, i.e. the plain \(\text{chol}(\Omega)\,z\) move (issue #895).

Kernel 2 — componentwise sweep: for multi-η models, follow the block move with \(\max(2, \lfloor \texttt{n\_mh\_steps} / n_\eta \rfloor)\) sweeps that perturb one coordinate at a time, \(\eta_j' = \eta_j + \delta_i^{\text{cw}} \cdot \sqrt{\Omega_{jj}} \cdot z\), holding the others fixed. Because the block proposal is shaped by \(\text{chol}(\Omega)\), once \(\Omega\) drifts toward a high correlation the block move can only travel along that near-degenerate direction; the single-draw Ω M-step then feeds the correlation back into \(\Omega\), and during the \(\gamma=1\) exploration phase this compounds into a runaway collapse toward a near rank-1 \(\Omega\) (every off-diagonal correlation \(\to \pm 1\), one variance \(\to 0\)). A componentwise move can always shift a single η independently of \(\Omega\)’s off-diagonals, so the sampled draws are not forced collinear. The kernel is skipped for single-η models (no off-diagonal to decorrelate).

Both kernels are symmetric in \(\eta\), so the proposal density cancels and the acceptance log-ratio is the difference of individual_nll values, which encodes the prior \(N(0, \Omega)\) plus the observation likelihood.

Acceptance: \(\min(1, \exp(\text{NLL}_{\text{current}} - \text{NLL}_{\text{prop}}))\). Target acceptance rate: 40%.

Reported acceptance is combined. The mh_accept_rate column in the optimizer trace (and the verbose banner) is the combined block + componentwise rate. Reporting the block kernel alone is misleading for FREM-scale \(\Omega\): as above, the block move reads 0% acceptance even when the componentwise sweep is mixing the chain perfectly well. If the combined post-burn-in rate stays below 1%, the E-step genuinely is not mixing and SAEM appends a warning to FitResult.warnings — the sampled ETAs never moved, so the \(\Omega/\sigma\) estimates are unreliable (issue #895).

HMC (Hamiltonian Monte Carlo, n_leapfrog > 0)

When n_leapfrog is set to a positive integer (e.g. 3), one HMC proposal replaces the n_mh_steps MH proposals per subject per iteration. HMC uses the gradient of the individual NLL to make longer, more directed moves through the posterior:

\[ H(\eta, p) = \text{NLL}(\eta) + \tfrac{1}{2} \|p\|^2 \]

with momentum \(p \sim N(0, I)\) and a standard velocity Störmer-Verlet (leapfrog) integrator. Acceptance is on \(\Delta H\), targeting ~65% acceptance.

Requirements: an analytical PK model (no ODE). The HMC gradient is the exact analytic Dual2 η-gradient (the same one FOCEI uses). A warning is emitted if n_leapfrog > 0 but the model is out of the analytic provider’s scope.

Per-subject fallback to MH: hmc_step silently falls back to MH for a subject when any of the following conditions hold: - The model uses an ODE ([odes] block present) - The model has no analytical PK path (no tv_fn — pure ODE-only models) - The Ω matrix has a non-finite log-determinant (degenerate variance) - The model is outside the analytic gradient’s scope (time-varying covariates, oral infusion, SS+reset, or LTBS). A η-dependent expression obs_scale is served analytically (since #486), so on a closed-form model it takes the HMC path — only LTBS combined with it falls back to MH.

In a single run with n_leapfrog > 0, different subjects can therefore use different samplers. The acceptance rate reported in verbose output and the optimizer trace is an aggregate across all subjects; in mixed HMC/MH runs the target (65% for HMC, 40% for MH) may not be meaningful for the aggregate. The n_mh_steps option governs the number of proposals for MH-fallback subjects even when n_leapfrog > 0.

The E-step sampling is parallelized across subjects using Rayon.

The startup banner reports the resolved E-step kernel on a sampler: line, e.g. sampler: Metropolis-Hastings random walk or sampler: HMC (3 leapfrog steps, Dual2 analytic gradients). If n_leapfrog > 0 but HMC is unavailable (an ODE model, or one outside the analytic gradient’s scope), the line says so and reflects the MH fallback. Because SAEM is sampling-based rather than gradient-driven, it does not print the gradient: line that FOCE/FOCEI use — the gradient_method option only governs the inner EBE/Hessian step (used for diagnostics, and consumed by a following imp stage), not the SAEM iterations themselves.

2. Stochastic Approximation Update

Update the sufficient statistic for \(\Omega\):

\[ S_2 \leftarrow (1 - \gamma_k^\Omega) \cdot S_2 + \gamma_k^\Omega \cdot \frac{1}{N} \sum_{i=1}^{N} \eta_i \eta_i^T \]

Damped Ω step. The θ/σ M-step uses the full \(\gamma_k\) (1.0 during exploration), but the Ω sufficient statistic uses a capped step \(\gamma_k^\Omega = \min(\gamma_k, 0.1)\). With the full \(\gamma=1\), Ω would be overwritten every exploration iteration by a single warm-started, not-yet-equilibrated MCMC draw; for a correlated block that one snapshot is biased toward the chain’s current correlation, and the bias feeds back through \(\mathrm{chol}(\Omega)\) into the next proposal — the same rank-1 runaway the componentwise kernel guards against, here attacked from the M-step side. Capping the Ω learning rate during exploration averages those draws (Robbins-Monro) and breaks the feedback while θ still moves at full speed. The cap applies during exploration only; in the convergence phase it is lifted and Ω uses the full decaying \(\gamma_k = 1/(k-K_1)\) — the same Robbins-Monro schedule as θ — so the SA estimate settles correctly (by then the chain is equilibrated, so the single-draw overwrite risk no longer applies).

3. M-Step for Omega (Closed Form)

\[ \Omega_k = S_2 \]

with structurally-zero entries (cross-block off-diagonals, standalone-vs-block off-diagonals, and all off-diagonals in a fully-diagonal Ω) zeroed out — the SA accumulator \((1/N) \sum_i \eta_i \eta_i^T\) is dense by construction, but the model declares which entries are free parameters. Without this projection the chain feeds spurious sampling correlations into the next iteration’s MH proposal Cholesky and Ω drifts toward rank-deficiency. Free diagonals are then floored at 1e-6 to keep them away from zero (which would collapse the per-eta proposal scale δ·chol(Ω)); this floor does not by itself guarantee positive-definiteness of a full block Ω.

Burn-in. This M-step is suppressed for the first omega_burnin iterations (default 20, clamped to n_exploration): Ω is held at its starting value while the MH chain warms up. The MH proposal scale is \(\delta_i \cdot \mathrm{chol}(\Omega)\), so Ω and the sampler are coupled. On sparse data (few observations per subject) a cold-start chain (η = 0, only n_mh_steps proposals) produces a tiny \((1/N) \sum_i \eta_i \eta_i^T\); with \(\gamma_1 = 1\) the M-step would install that as Ω on iteration 1, shrinking the proposal, which keeps the chain near zero — a self-reinforcing collapse that dumps between-subject variability into the residual error. The SA statistic \(S_2\) is still refreshed each burn-in iteration (at the damped rate \(\gamma_k^\Omega\), so it is a running average of the warming chain), so the first Ω update after burn-in reflects the warmed-up chain rather than the cold-start spread. Set omega_burnin = 0 to disable the burn-in; note that the damped Ω step above now also guards this cold-start collapse on its own (it is a strict generalisation — continuous rather than for the first omega_burnin iterations only), so disabling the burn-in no longer reproduces the collapse by itself.

4. M-Step for Theta and Sigma (Optimization)

Minimize the conditional observation negative log-likelihood with ETAs held fixed:

\[ \sum_{i=1}^{N} \sum_{j=1}^{n_i} \left[ \frac{1}{2} \log V_{ij} + \frac{1}{2} \frac{(y_{ij} - f_{ij})^2}{V_{ij}} \right] \]

When mu_referencing = true (the default), ferx detects lognormal parameters from the [individual_parameters] block and applies the closed-form EM update for those thetas instead of running NLopt:

\[ \log \theta_j \leftarrow \log \theta_j + \gamma_k \cdot \overline{\eta_j} \]

where \(\overline{\eta_j} = (1/N) \sum_i \eta_{i,j}\) is the empirical mean of the post-MH random effects for the eta paired with \(\theta_j\). For a log-mu-referenced model where \(\log P_i = \log \theta + \eta_i\) with \(\eta_i \sim N(0, \omega^2)\), this is exactly the M-step that maximises the complete-data log-likelihood, scaled by the SA step size \(\gamma_k\). After the update the etas are re-centred by the same shift so they remain deviations from the new \(\log \theta_j\).

NLopt still runs for any remaining thetas (non-mu-referenced) and for sigma — the closed-form-updated thetas are pinned at their new values for the NLopt call.

The NLopt M-step uses BOBYQA (derivative-free trust-region with quadratic interpolation). The earlier gradient-based SLSQP path was found to lock onto one side of the Emax-Hill identifiability ridge on the dense-Emax PKPD benchmark (under-estimating EMAX by ~40% at virtually identical OFV); BOBYQA’s quadratic trust-region exploration lands much closer to truth and ~40% faster on that benchmark (no FD-gradient eval per parameter), while remaining numerically equivalent (ΔOFV < 0.1) on simpler PK-only models.

When mu_referencing = false, the full NLopt M-step runs for all thetas as before.

WarningNon-mu-referenced parameters slow SAEM convergence

Whenever the estimation chain runs SAEM, ferx warns if any individual parameter’s random effect is not mu-referenced — for example CL = TVCL + ETA_CL (or any non-loglinear form) instead of CL = TVCL * exp(ETA_CL). Such parameters cannot use the closed-form EM update above and fall back to the slower NLopt M-step, which can strongly slow convergence. The warning fires independently of the mu_referencing fit option: turning the option off does not remove the risk — it removes the mu-centering that mitigates it — so the diagnostic is most warranted in exactly that case. Prefer log-mu-referenced forms (P = TV * exp(ETA)) where the model allows.

ImportantA theta with no ETA at all is worse than slow — it can drift

The warning above covers a parameter whose eta ferx could not mu-reference. A fixed-effect-only theta — one that carries no eta anywhere in [individual_parameters] — is a distinct and more damaging case, and SAEM warns about it separately:

SAEM: estimated parameter(s) [TVFRD1] have NO associated ETA, so they are not
mu-referenced and are moved only by the η-frozen numerical M-step. ...

Such a theta never receives the closed-form log θ += γ·mean(η) shift, which is an exact Robbins-Monro average. Its only channel is the numerical M-step, which re-maximises the conditional observation likelihood against a single MCMC η draw. Assigning that maximiser outright is argmax of one draw rather than the SA average of E[argmax], so the estimate carries a Monte-Carlo bias that does not decay with iteration count, and it can settle a long way from the marginal optimum — dragging correlated typical values with it.

Measured on a FREM iiv_on_ruv model (475 subjects, 12 etas) whose absorption fraction TVFRD1 splits a zero-order duration and a first-order rate (D1 = MAT*(1-TVFRD1), KA = 1/(MAT*TVFRD1)) with no eta of its own. The marginal −2logL optimum for TVFRD1 on this dataset is ≈ 0.29:

Run TVCL TVV TVMAT TVFRD1
NONMEM IMP (reference) 7.060 133.8 2.680 0.394
ferx IMP 7.068 135.5 2.751 0.311
ferx IMPMAP 7.046 134.2 2.780 0.318
ferx SAEM, before the #1011 damping 7.169 143.2 3.020 0.039
ferx SAEM, before, started at the IMPMAP solution 7.046 136.4 3.191 0.065
ferx SAEM, current 7.123 139.5 2.686 0.290
ferx SAEM, FRD1 = TVFRD1*exp(ETA_FRD1), ω² = 0.01 7.181 142.5 2.646 0.313
ferx SAEM, TVFRD1 held FIX at 0.394 7.025 131.9 2.770 (fixed)

Restarting SAEM at the right answer did not rescue the old behaviour — the M-step walked the parameter back down — so this was never a starting-value problem.

The remedy depends on what the parameter is. A typical value should be put in a mu-referenceable form (P = TVP * exp(ETA_P), with a small — optionally FIX — omega; ferx applies the mu-referencing automatically). A parameter that has no eta to give — a covariate coefficient, an allometric exponent, a structural constant — should be held FIX, or the fit cross-checked against FOCEI/IMPMAP.

The advisory deliberately stays quiet about two kinds of theta, because “has no ETA” would be false for them:

  • Mixing coefficients of a mixture model. The numerical θ/σ M-step never moves them — they do not enter the residual/eta likelihood, and the mixing M-step overwrites them from the class responsibilities — and a mixing expression is forbidden from depending on an eta in the first place.
  • MIXNUM-switched class typical values, and log-mu-referenced thetas dropped to the identity scale by a negative lower bound. These carry an eta in the model text and are on the numerical M-step for a reason SAEM already states in its own warning; repeating “add an ETA” on top of that would be both false and contradictory.

What SAEM does about it now

The numerical θ/σ M-step result is blended in as a Robbins-Monro step θ ← θ + γ_θ·(θ* − θ) rather than assigned outright, with γ_θ capped at 0.03 during exploration and following the full decaying γ = 1/(k−k1) in convergence. This is the θ-side counterpart of the cap SAEM already applied to the Ω sufficient statistic, and it exists for the same reason: a single un-equilibrated MCMC draw must not be able to carry a parameter away.

The cap is the mstep_damping fit option, and defaults to 0.03:

[fit_options]
  method        = saem
  mstep_damping = 0.01   # smaller damps harder; 1.0 turns the damping off

Values must lie in (0, 1]. Smaller damps harder — on the model below, 0.01 moves TVFRD1 to 0.355 and 0.003 to 0.357, against NONMEM’s 0.394 — but tighter caps also slow legitimate exploration, and here they pull TVV off (129.9 at 0.01 against NONMEM’s 133.8). mstep_damping = 1.0 disables the damping entirely, in both phases, restoring the pre-0.3.x behaviour exactly; that is the escape hatch if a model fitted better without it.

1.0 is therefore a sentinel, not a cap value, and the option is deliberately discontinuous there: mstep_damping = 0.999 still buys the full convergence schedule γ = 1/(k−k1), while 1.0 buys none of it. There is no setting that means “cap exploration a little and leave convergence alone” — the option only ever governs the exploration cap — and “off” has to restore the whole pre-damping trajectory to be a useful escape hatch. (Lifting only the exploration cap left the model above at TVFRD1 0.047 rather than its true undamped 0.039.)

The damping is applied only to the channel the bias lives in: a fit is damped when NLopt is left estimating a theta that is not mu-referenced, which is exactly the situation the advisory above warns about. A theta that is FIX, that is pinned out by the mu-reference shift, or that is mu-referenceable at all (so the exact closed-form log θ += γ·mean(η) update is available to it) is left on the undamped update, and a fit whose every estimated theta is one of those is unaffected down to the last bit — and setting mstep_damping on such a model warns that it has no effect rather than accepting it silently.

Sigma is damped with the same γ_θ, on the same gate, even though the gate is a statement about theta. That is not an oversight: the numerical M-step is a single NLopt problem over the concatenated [θ; σ] vector, so θ* and σ* are one joint maximiser of the same frozen-η likelihood. Blending theta at γ_θ while assigning sigma outright would take an inconsistent partial step — sigma would jump to the value that is optimal for the theta the M-step wanted, while theta is held near the theta it had, leaving sigma to absorb the very misfit the damping just stopped theta from fixing. So a model that trips the gate has its sigma damped too, including when the theta that trips it is a covariate coefficient or an allometric exponent rather than a typical value.

One case is deliberately left out of scope: under mu_referencing = false every theta goes through the numerical M-step, but a theta written in a mu-referenceable form is not the fixed-effect-only channel the damping was calibrated against, so such a fit keeps the undamped update and mstep_damping reports that it had no effect (naming that reason).

“Mu-referenced” here means ferx detected the pairing — the TVX * exp(ETA_X) and logit forms. A parameter that carries an eta in a form ferx cannot pair with a single theta (an additive X = TVX + ETA_X, or a covariate model whose typical value is not log-linear in one theta) counts as un-referenced and is damped, which is the intended side to err on: that theta has no closed-form shift either, so the numerical M-step really is its only mover.

Mixture models are also excluded. A MIXNUM-switched typical value goes through the same numerical M-step, but it is solving a different problem — the class typical values must separate from a common start, and the class assignments only settle once they have. Damping that excursion stalls the separation, so mixtures keep the undamped update.

Averaging the M-step objective over the last K η draws instead — the obvious alternative — is both weaker and far more expensive: it reaches TVFRD1 0.081 at K = 50 and 0.130 at K = 150 (15× the wall time), because consecutive SAEM draws are heavily autocorrelated. Damping the iterate accumulates over every iteration rather than a K-window, which is why it wins on both axes.

Damping reduces the Monte-Carlo bias but does not remove it, so the advisory still fires and the two structural remedies remain the better answer:

  • Attach an eta — P = TVP * exp(ETA_P) with a small omega (optionally FIX). ferx then mu-references it automatically and the exact closed-form update takes over.
  • Hold it FIX at a value estimated by another method. Every theta then lands within 3% of NONMEM, TVV included.

Otherwise estimate the model with FOCEI or IMPMAP, which do not depend on this channel.

The number of NLopt evaluations saved is stored in FitResult::saem_mu_ref_m_step_evals_saved, accumulated across SAEM iterations as 2 × mstep_maxiter × n_mu_ref_pairs per outer step (one finite-difference probe pair per pinned mu-ref dimension, capped at mstep_maxiter NLopt gradient requests). The field is None when mu-referencing is off or method ≠ SAEM.

When n_leapfrog > 0, FitResult::saem_n_subjects_hmc records how many subjects used HMC at least once during the E-step (the remainder used MH fallback). The field is None for MH-only runs. The fit YAML also emits saem_n_subjects_hmc and saem_n_subjects_mh when the field is Some.

5. Adaptive Step Sizes

Every adapt_interval iterations, the per-subject step sizes \(\delta_i\) (MH) or leapfrog step sizes (HMC) are adjusted based on acceptance rate:

  • If acceptance rate exceeds the target: increase \(\delta_i\) by 10% (up to 5.0)
  • If acceptance rate falls below the target: decrease \(\delta_i\) by 10% (down to 0.01)

The block kernel and the componentwise kernel carry independent per-subject scales (the componentwise scale is adapted toward the higher 1-D optimum), and the kappa MH (IOV) scale is adapted the same way.

Target acceptance rates: 40% for the block MH kernel, 44% for the componentwise kernel, 65% for HMC.

Post-SAEM Finalization

After the SAEM iterations complete:

  1. EBE Refinement: Run the standard FOCE inner loop (BFGS optimization) warm-started from the SAEM ETAs to obtain final empirical Bayes estimates
  2. Combined-error additive-collapse repair: for combined(PROP, ADD) residual-error models, only when SAEM has driven the additive component ADD onto its lower bound, ferx runs a final FOCEI marginal-likelihood polish from the SAEM estimates (retrying from the model’s initial parameters if ADD is still pinned) and adopts it only if it lowers the marginal OFV. Fits whose ADD converged to a healthy non-zero value are left untouched. This guards against point-η SAEM M-steps overfitting the low-concentration tail and collapsing ADD when the marginal likelihood identifies a non-zero additive term. It is a safety net, not a substitute for checking ADD identifiability yourself (RSE, σ-correlation, profile likelihood); on sparse data prefer the importance-sampling −2LL over the Laplace OFV when judging whether the additive term is real.
  3. FOCE OFV: Compute the objective function using the FOCE/Laplace approximation, so AIC and BIC are directly comparable with FOCE results
  4. Covariance Step: Optionally compute standard errors via finite-difference Hessian (same method as FOCE)
  5. Diagnostics: Compute PRED, IPRED, CWRES, IWRES for each subject

For sparsely-sampled data where the Laplace OFV is biased, you can append an importance-sampling stage that estimates −2 log L by Monte Carlo:

method = [saem, imp]

See Importance Sampling (IMP).

Conditional Distribution (conditional mode vs. distribution)

The post-SAEM finalization above produces the conditional mode of each subject’s random effects — the empirical Bayes estimate (EBE), the single most probable \( _i \). That is a point estimate. SAEM’s MCMC E-step is, however, already sampling each subject’s full conditional distribution \( p(_i y_i; ) \); the mode discards everything but its peak.

Set conddist = true to run an opt-in post-fit pass that characterises that distribution. With the population parameters fixed at their converged values, the same MH kernels (block, componentwise, and the per-occasion kappa kernel for IOV) are re-run per subject — warm-started at the EBE mode — and the draws are accumulated rather than discarded. The pass reports, per subject:

  • the conditional mean \( [_i y_i] \),
  • the conditional SD \( (_i y_i) \),
  • optionally the raw draws (conddist_keep_samples = true), and
  • a distribution-based η-shrinkage, \( 1 - _i(_i)/\), reported alongside the usual mode-based shrinkage_eta.

This mirrors the conditional-mode vs. conditional-distribution distinction in saemix (map.saemix vs conddist.saemix; Comets, Lavenu & Lavielle, J. Stat. Soft. 80(3), 2017) and Monolix (the “Conditional Mode” vs “Conditional Distribution” tasks). Why prefer the distribution for diagnostics: EBEs and conditional means are shrunk toward the population, so η–covariate and η–η relationships built on them can be hidden or fabricated; samples from the conditional distribution are not shrinkage-biased.

Results are exposed on FitResult.cond_dist and written by the CLI to {model}-conddist.csv (ID, ETA, COND_MEAN, COND_SD, COND_MODE), with the raw draws in {model}-conddist-samples.csv when retained.

Validation against saemix and NONMEM

On the bundled warfarin data (10 subjects, 1-cpt oral, log-normal CL/V/KA, proportional error), all three engines converge to identical population parameters (TVCL 0.1327, TVV 7.737, TVKA 0.811), and ferx’s per-subject conditional distribution agrees with both references to within Monte-Carlo noise.

vs saemix (conddist.saemix):

Quantity (per-subject η) corr max|diff| RMSE
conditional mean 1.0000 0.0029 0.0011
conditional SD 0.9913 0.0024 0.0009
mode / MAP 1.0000 0.0002 0.0001

vs NONMEM ($EST METHOD=SAEM then METHOD=IMP EONLY=1; conditional moments read from the .phi file — PHI(k) − log θ_k is the η conditional mean, sqrt(PHC(k,k)) the conditional SD):

Quantity (per-subject η) corr max|diff| RMSE
conditional mean 1.0000 0.0012 0.0004
conditional SD 0.9987 0.0017 0.0004

(ferx conddist_nsamp = 2000. Comparison scripts: tests/reference/saem_conddist/bench_saemix_conddist.R, tests/reference/saem_conddist/parse_nm_conddist.R.)

Inter-Occasion Variability (IOV)

method = saem supports models with kappa declarations (n_kappa > 0). IOV is handled with a per-occasion Gibbs Metropolis-Hastings step interleaved with the standard eta MH:

  • E-step: After sampling \(\eta\), one MH proposal is made for each occasion’s \(\kappa_k\). Both the eta and kappa samplers target the correct conditional distributions: \(p(\eta | \kappa, \theta, \text{data})\) and \(p(\kappa_k | \eta, \theta, \text{data})\) respectively.
  • SA update: \(S_2^{\text{iov}} \leftarrow (1 - \gamma_k) S_2^{\text{iov}} + \gamma_k \cdot \frac{1}{N_{\text{occ}}} \sum_i \sum_k \kappa_{ik} \kappa_{ik}^T\)
  • M-step: \(\Omega_{\text{iov}} = S_2^{\text{iov}}\) (analytic update, same structure as the BSV omega M-step).

No additional configuration is required; method = saem works for both BSV-only and IOV models.

Models with no random effects are rejected

SAEM is an EM algorithm over the random effects. A model that declares none (n_eta = 0 — see no omega at all) has an empty E-step: there is nothing to sample and the stochastic approximation has nothing to average, leaving an M-step that is plain maximum likelihood on θ/σ. Left to run, it terminates on its own criterion slightly short of the optimum — on the one_cpt_iv zero-Ω anchor it reports converged at OFV −269.599 where FOCEI reaches −269.637.

Rather than return a quietly-wrong answer, method = saem errors at n_eta = 0:

method = saem requires at least one random effect (n_eta = 0). SAEM is an EM over
the random effects, so with none declared its E-step is empty and it only
approaches the objective that FOCE/FOCEI minimise exactly. Use method = foce,
focei, or laplace for a fixed-effects-only (naive-pooled) model.

The check runs in the same validation pass as ferx check, and fires if saem appears anywhere in a methods = [...] chain — so methods = [saem, focei] fails up front rather than after the SAEM stage has already run.

imp, impmap and bayes are rejected at n_eta = 0 in the same pass, under the code E_METHOD_NO_RANDOM_EFFECTS (#1007) — each integrates over the random effects, and with none declared the marginal likelihood collapses to the observation likelihood FOCE/FOCEI already minimise exactly. Like the SAEM guard this fires anywhere in a chain, so methods = [focei, imp] fails up front instead of after the FOCEI stage has run. Each of the three still refuses at run time as well, which is what a direct fit() caller bypassing ferx check sees.

Configuration

[fit_options]
  method        = saem
  n_exploration = 150      # Phase 1 iterations
  n_convergence = 250      # Phase 2 iterations
  n_mh_steps    = 20       # block-MH steps per subject per iteration (ignored when n_leapfrog > 0)
  n_leapfrog    = 0        # Set > 0 (e.g. 3) to use HMC instead of MH
  adapt_interval = 50      # Step-size adaptation frequency
  omega_burnin  = 20       # Iterations to hold Ω fixed while the chain warms up
  seed          = 12345    # RNG seed for reproducibility
  covariance    = true     # Compute standard errors

  # Conditional-distribution pass (opt-in; off by default)
  conddist              = true    # Estimate p(η_i | y_i) per subject after the fit
  conddist_nsamp        = 2000    # Retained MCMC draws per subject
  conddist_burnin       = 500     # Burn-in draws discarded before accumulation
  conddist_keep_samples = false   # Retain the raw draws (writes -conddist-samples.csv)

Tuning Guide

Not Converging

  • Increase n_exploration (e.g., 300) to give more time for basin finding
  • Increase n_convergence (e.g., 500) for a longer averaging window
  • Raise n_mh_steps further (e.g. 30-50) for better mixing in the E-step on hard surfaces — the default 20 is calibrated to escape the basin trap observed on Emax PKPD with stressful initial values (see below) and to drive the componentwise kernel that prevents block-Ω collapse, but ODE-with-Form-C readouts may need more proposals to fully decorrelate samples between M-step calls.

PD-curve thetas collapse on cold start (Emax / sigmoid readouts)

A failure mode specific to models that read population thetas through a Form C [scaling] block (e.g. y[CMT=N] = E0 + EMAX * effect^GAMMA / ...) and score them via a per-CMT additive [error_model]: from stressful initial values (e.g. 1.5× truth) the M-step can lock the PD-curve thetas into a degenerate basin where E0 → 0, EMAX and EC50 blow up, and GAMMA collapses below 1. The likelihood at the bad basin is only modestly worse than at truth (~150 OFV units on a 100-subject benchmark), so SAEM doesn’t back out on its own.

The underlying cause is MCMC sample correlation: with an early default n_mh_steps = 3 the chain didn’t decorrelate enough between SAEM outer iterations, so the single-draw stochastic M-step received sticky correlated ETAs that biased the population-θ update toward the basin. The default was raised to 10 and then to 20 (alongside the componentwise kernel and the damped Ω step), which resolves this reliably across seeds at modest extra wall on the affected model and ~0% on simpler PK-only models.

If you still see this signature (E0 hitting its lower bound, EMAX large, EC50 large) on a related Emax/Hill model:

  1. Try n_mh_steps = 20 or n_mh_steps = 50
  2. Warm-start from a FOCEI fit (method = [focei, saem])
  3. Run with several seeds and keep the lowest OFV

Ω Collapses / Residual Error Inflates

On sparse data (few observations per subject) the variance components can collapse toward zero on the first iterations while the residual error absorbs the between-subject variability (e.g. tiny omega with a large additive sigma). The default omega_burnin = 20 and the damped Ω SA step guard against this by keeping Ω near its starting value while the chain warms up. If it still occurs, raise omega_burnin (e.g. 40) and/or n_mh_steps so the chain reaches a representative spread before Ω is first estimated, then polish with method = [saem, focei].

Block Ω correlations near ±1 (rank-1 collapse)

For a block_omega (correlated) random-effects block, a faulty E-step/M-step coupling can drive every off-diagonal correlation toward ±1 while one variance collapses toward zero — a near rank-1 Ω that FOCEI on the same data does not show. The mechanism is the block proposal (preconditioned by chol(Ω)) plus the single-draw Ω M-step feeding correlation back into Ω during the γ=1 exploration phase. ferx guards this by default with two mechanisms — the componentwise MH kernel and the damped Ω SA step described above — so on a poorly-identified 2-cpt model the SAEM Ω now matches the FOCEI/NONMEM estimate (e.g. corr(CL,V1) ≈ 0.67, corr(V1,V2) ≈ 0.4) across seeds instead of collapsing to ≈0.99. If you still see inflated block correlations, raise n_mh_steps (this also raises the componentwise sweep count n_mh_steps / n_eta) and run several seeds.

IIV on residual error (iiv_on_ruv): σ × ω_RUV ridge

With iiv_on_ruv = ETA, the residual is \(Y = f + \varepsilon \exp(\eta_{\text{RUV}})\), so the residual variance is \(\sigma^2 \exp(2\eta_{\text{RUV}})\). This form is invariant to \(\eta_{\text{RUV}} \to \eta_{\text{RUV}} - c,\ \sigma \to \sigma\,e^{c}\) — \(\sigma\) and a shift of \(\eta_{\text{RUV}}\) are a degenerate pair. Unlike a mu-referenced structural \(\eta\) (whose mean is folded into its typical value \(\theta\) each M-step), \(\eta_{\text{RUV}}\) has no typical-value \(\theta\), so nothing absorbs its mean: it drifts along that degenerate direction, and because \(\omega_{\text{RUV}}\) is estimated as \(\overline{\eta_{\text{RUV}}^2}\) the drift injects a spurious \(\text{mean}^2\) term that pumps \(\sigma\) and \(\omega_{\text{RUV}}\) up together — the runaway reported in issue #895 (\(\omega_{\text{RUV}}\) toward ~49, \(\sigma\) toward its e⁵ ceiling), worst on FREM models with an extreme \(\Omega\)-diagonal scale range.

ferx fixes this by treating \(\sigma\) as the residual-scale typical value: each iteration \(\eta_{\text{RUV}}\) is re-centred to zero mean, absorbing the shift into \(\sigma\) (\(\eta_{\text{RUV}} \mathrel{-}= \overline{\eta_{\text{RUV}}}\), every RUV-scaled \(\sigma \mathrel{\times}= e^{\overline{\eta_{\text{RUV}}}}\)). Each subject’s residual variance is exactly unchanged, so the likelihood is untouched, but \(E[\eta_{\text{RUV}}]=0\) is restored and \(\omega_{\text{RUV}}\) measures the true variance. On the 475-subject FREM reprex this converges to \(\omega_{\text{RUV}} \approx 0.28\), \(\sigma \approx 0.20\) from both a too-small and a too-large \(\sigma\) start (NONMEM: 0.28 / 0.18), where it previously ran to \(\omega_{\text{RUV}} \approx 49\). Re-centring is applied only when every RUV-scaled \(\sigma\) component is free to absorb the shift — this keeps each subject’s residual variance exactly unchanged. If any residual \(\sigma\) component is FIXed (e.g. a fixed additive term of a combined error), scaling only the free component would leave the fixed part uncompensated, so re-centring is skipped and the growth caps below carry the load (a FIXed \(\sigma\) already partially pins the mean).

Two growth caps remain as belt-and-braces backstops, both no-ops on a well-behaved fit:

  • σ growth cap — each free RUV-scaled residual \(\sigma\) is capped at \(\approx 20\times\) its scale (a post-M-step clamp). A \(\sigma\) whose own user upper bound is already tighter carries no separate cap — NLopt enforces that bound directly.
  • ω_RUV growth cap — the RUV \(\Omega\) variance is capped at \(\approx 20\times\) its scale, applied as a correlation-preserving rescale of the RUV row/column so a block \(\Omega\) stays positive-definite. A FIXed off-diagonal covariance with the RUV eta is left exactly as declared.

The reference “scale” is the larger of the starting value and the data-informed value the fit reaches by the end of the exploration phase, so a run started from a \(\sigma\)/\(\omega_{\text{RUV}}\) guess many-fold below the truth is not spuriously clamped. If either cap ever binds (it should not, now that the drift is removed), SAEM warns that the split is weakly identified and suggests fixing \(\sigma\).

Slow Convergence

  • Decrease n_exploration and n_convergence if parameters stabilize early
  • Use adapt_interval = 25 for faster step-size adaptation

Reproducibility

  • Always set seed for reproducible results
  • Different seeds will produce slightly different estimates due to the stochastic nature of the algorithm

Output

The SAEM iteration progress is printed to stderr:

SAEM: 10 subjects, 3 ETAs, 400 total iter (150 explore + 250 converge)
  SAEM iter    1/400 [explore] γ=1.000  condNLL=95.244
  SAEM iter   50/400 [explore] γ=1.000  condNLL=56.705
  SAEM iter  150/400 [explore] γ=1.000  condNLL=46.071
  SAEM iter  200/400 [converge] γ=0.020  condNLL=36.799
  SAEM iter  400/400 [converge] γ=0.004  condNLL=38.096
SAEM iterations complete. Computing final EBEs and OFV...
SAEM completed. Final OFV = ...
Running covariance step...

The Final OFV = ... line is printed before the covariance step starts (#893). SAEM only learns its OFV at the very end (the final FOCE approximation), and the covariance matrix is often the most expensive part of the run — so seeing the OFV first lets you interrupt (Ctrl-C) when the OFV already rules the run out, instead of waiting for a covariance matrix you won’t use.

Why γ (gamma) is shown

\(\gamma_k\) is the stochastic-approximation step size, not a model quantity — it is intrinsic to the SAEM algorithm and tells you which phase the run is in and how aggressively the estimates are still moving:

  • Exploration ([explore], \(k \le K_1\)): \(\gamma_k = 1\). Each iteration fully replaces the running sufficient statistics, so the chain roams freely toward the basin of the MLE.
  • Convergence ([converge], \(k > K_1\)): \(\gamma_k = 1/(k - K_1)\) decays toward zero. Updates shrink into a decreasing-weight average that damps the Monte-Carlo noise so the estimates settle. The printed value (e.g. 0.020, 0.004) is exactly this decaying weight applied to the Ω, θ, and σ updates each iteration.

Why condNLL and not OFV

During the iterations ferx prints condNLL, the conditional (joint) negative log-likelihood summed over subjects, evaluated at the current MH/HMC-sampled etas:

\[ \text{condNLL} = \sum_{i=1}^{N} \text{NLL}(\eta_i^{\text{sampled}}) \]

This is a cheap per-iteration progress signal. It is not the marginal objective function value (OFV): it is evaluated at one stochastic draw of the random effects rather than integrated over their distribution, so unlike the FOCE/FOCEI outer-loop OFV it is noisy, will not decrease monotonically, and is not comparable across runs for model selection.

The true marginal OFV (\(-2 \log L\) via the Laplace approximation, directly comparable with FOCE for AIC/BIC) is expensive and is therefore computed only once, after the iterations finish — this is the Final OFV = ... line. See Post-SAEM Finalization.

(NLL = negative log-likelihood, i.e. \(-\log L\). NONMEM’s OFV = -2 \log L is essentially 2 × NLL plus a constant.)

condNLL should generally decrease during the exploration phase and stabilize during convergence.