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).

Choosing `n_mh_steps`

n_mh_steps defaults to auto: the block count is sized from the dataset as

\[ \texttt{n\_mh\_steps} = \mathrm{clamp}\!\left(\left[\,2.5 \cdot \frac{n_{\text{obs}}}{n_{\text{subjects}} \cdot n_\eta}\,\right],\; 6,\; 20\right) \]

i.e. from the observations each random effect has to be informed by. The resolved count and where it came from are printed under verbose = true, and an explicit n_mh_steps = <n> overrides the rule entirely.

Two scope notes. With n_leapfrog > 0 the count is not ignored: HMC replaces the block kernel for the subjects it can serve, but n_mh_steps still governs the MH-fallback subjects and still sizes the componentwise sweep, which runs after an HMC proposal as well. And method = bayes reads the option but not the rule — under auto its η block keeps the historical fixed 20, because the rule is calibrated on SAEM quantities and, per the kernel split below, what makes a low count safe on dense data is SAEM’s componentwise sweep, which that sampler does not run. An explicit count applies to both.

The reason it is not a constant is the step-scale controller. On a sparse population it reaches its 40 % acceptance target at any proposal count — measured at 0.39–0.46 on cefepime, vancomycin, busulfan and pembrolizumab for every count from 2 to 20 — so each proposal is worth the same fraction of a move, mixing per iteration is linear in the count, and so is the cost. Twenty proposals then buy nothing that six do not. The one benchmark where extra proposals pay for themselves is the one where the controller cannot reach its target: the Emax PKPD model below realises 0.246, because its conditional $p(\eta_i \mid y_i)$ is far narrower than the $N(0,\Omega)$ prior the proposal is scaled by. Observations per subject per η is the cheapest proxy for that narrowing, which is why the rule reads it.

On that dense model it is specifically kernel 2 that earns the count, not kernel 1. Breaking the coupling experimentally (6 block proposals with the 10 componentwise sweeps that 20 implies, against 20 block proposals with the 3 sweeps that 6 implies, 6 seeds each) puts the sweeps arm on 20’s numbers — d_final 0.043 and cross-seed stability 0.045, against 20’s own 0.029 / 0.044 — and the block arm on 6’s — 0.102 / 0.138 against 0.091 / 0.132. Raising n_mh_steps is how a model asks for more sweeps today, and at equal quality the two spellings cost the same wall here (28.4 s vs 29.0 s), so the rule sizes the block count and lets the sweep count follow (#1466 tracks separating them).

Measured on ferx e7343462, 6 seeds, 150 + 250 iterations, IS −2 log L from the imp readout and d_final against a 300/700/40 long run of the same build:

dataset	obs/η	`auto`	IS −2 log L (auto vs 20)	`d_final` (auto vs 20)	SAEM-stage CPU
cefepime (458 subj, 3 η)	0.51	6	4296.3 vs 4294.9 (seed sd 1–2)	0.078 vs 0.075	−35 %
vancomycin (100 subj, 3 η)	0.91	6	1135.3 vs 1135.0 (sd 0.4–0.9)	0.037 vs 0.033	−41 %
pembrolizumab (303 subj, 4 η)	1.02	6	5983.8 vs 5981.2 (sd 4–10)	0.179 vs 0.115	−29 %
busulfan (600 subj, 3 η)	2.92	7	60044 vs 60045 (sd 3)	0.040 vs 0.041	−38…−40 %
Emax PKPD (100 subj, 2 η)	8.0	20	unchanged — the rule returns the old default		0 %

CPU is the SAEM stage’s user+sys, 3 seeds on an uncontended lane at --threads 2. Busulfan’s 7 is quoted as a bracket because the count it resolves to sits between two measured cells (6 → 4.40 s and 8 → 4.23 s, against 7.10 s at 20) whose costs differ by less than the run-to-run spread.

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).

Damped σ step. The residual σ that rides the numerical M-step gets its own capped step, $\gamma_k^\sigma = \min(\gamma_k,\,0.2,\,\gamma_k^\theta)$, applied to $\sigma^2$ rather than to $\log\sigma$ (issue #1445). Until then it was the one SAEM statistic with no stochastic approximation at all — assigned outright in both phases — so the reported σ was a single draw of the M-step maximiser rather than its average. The cap is on in both phases, which is what distinguishes it from the Ω step: Ω accepts the $\gamma=1$ assignment at the first convergence iteration because $(1/N)\sum_i \eta_i\eta_i^T$ over $N$ subjects is a well-determined statistic, while the σ maximiser of a minority variance component is not. The additive half of a combined(PROP, ADD) model on sparse data is exactly that: its single-draw conditional maximiser is boundary-heavy, so an assigned σ reports whichever value the last iteration happened to draw. Beyond $k-K_1 = 5$ the cap no longer binds and the full decaying $\gamma_k = 1/(k-K_1)$ takes over, which is what lets the SA estimate settle. Averaging on the variance scale is deliberate — it is the scale the scalar residual sufficient statistic already averages on, and a geometric mean of a sequence with mass near zero collapses where the arithmetic mean of $\sigma^2$ does not (measured: $\log\sigma$ averaging leaves the #1445 fixture’s $\sigma_{\text{add}}$ at 0.47 against a truth of 1.8, where the $\sigma^2$ average reaches 1.4). The third term is one-sided: σ never steps faster than θ, but the two are not always equal. They are equal during exploration whenever mstep_damping is at or below 0.2 (the iiv_on_ruv default of 0.03 is), and from the fifth convergence iteration on, where $1/(k-K_1) \le 0.2$. They differ over the four hand-off iterations $k-K_1 \in \{1,2,3,4\}$ — θ takes 1, ½, ⅓, ¼ while σ takes 0.2 — and during exploration for any mstep_damping above 0.2, where σ is the slower of the two. So no configuration, iiv_on_ruv included, is bit-identical to before; what is preserved is that σ is never given a larger step than the θ it shares a maximiser with.

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] \]

Closed-form EM update for mu-referenced thetas

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

\[ g(\theta_j) \leftarrow g(\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$, and $g$ is the mu-scale link. For a mu-referenced model where $P_i = g^{-1}(g(\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$ — and it is link-independent: the same shift is the maximiser for $g = \log$, $g = \mathrm{id}$ and $g = \mathrm{logit}$ alike. After the update the etas are re-centred by the same shift so they remain deviations from the new $g(\theta_j)$.

What decides eligibility is therefore not the link but whether the packed theta — the quantity the optimiser actually steps — is $g(\theta)$. ferx packs a theta as $\log \theta$ when its lower bound is non-negative and as $\theta$ otherwise, so:

Form	mu scale $g(\theta)$	Closed-form EM update
`P = THETA * exp(ETA)`, `exp(log(THETA) + ETA)`	$\log \theta$	yes (theta is log-packed)
`P = inv_logit(THETA + ETA)`	$\theta$	yes — a logit-scale theta has a negative lower bound, so it is packed as-is (#918)
`P = inv_logit(logit(THETA) + ETA)`	$\mathrm{logit}\,\theta$	no — no packing has that scale
`P = THETA + ETA`	$\theta$	no — routed to NLopt

The two excluded forms are still mu-referenced (they still centre the inner/E-step search); they just take the numeric M-step. This matters most for bounded $(0,1)$ parameters — bioavailability, pathway fractions, $E_{max}$ fractions — in models that also carry IIV on the residual error, where the numeric M-step works from noticeably noisier conditional samples: declare the theta on the logit scale (theta LOGIT_F(-0.477, -10.0, 10.0)) and write F = inv_logit(LOGIT_F + ETA_F).

The packing rule cuts both ways, and a mismatch is reported rather than silently routed: a lognormal theta declared with a negative lower bound is packed on the identity scale, and a logit-scale theta declared with a non-negative lower bound (theta LOGIT_F(0.5, 0.0, 5.0)) is packed on the log scale. Neither packed value is the mu, so the closed form does not apply; SAEM and IMP/IMPMAP estimate such a theta through the numerical M-step and emit an advisory naming it, with the bound to change. Under mu_referencing = false SAEM runs no closed-form shift at all and suppresses the advisory; single-population IMP/IMPMAP is different — its shift is EM-mandatory, so it runs (and advises) independently of that option, as the IMP/IMPMAP page describes.

A theta that anchors more than one eta (F1 = inv_logit(LOGIT_F + ETA_F1) and F2 = inv_logit(LOGIT_F + ETA_F2), or the lognormal CL = TVP*exp(ETA_CL) / V = TVP*exp(ETA_V)) is also routed to the numerical M-step with an advisory: the closed form shifts a packed theta by one eta mean and pins it, and no single such shift is the joint maximiser. Give each eta its own typical value to get the closed-form update.

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

Covariate mu-references: typical values that read several thetas

A typical value that reads two or more thetas has no single anchor:

CL = (TVCL + (CRCL - 90) * TH_CRCL) * exp(ETA_CL)   # additive covariate
CL = TVCL * (WT / 70) ^ TH_WT * exp(ETA_CL)           # estimated exponent
F  = inv_logit(LOGIT_F + TH_SEX * SEX + ETA_F)        # logit-scale covariate

Before #619 every theta in such a value sat on the numerical M-step, which holds each subject’s sampled η fixed. That is the wrong coordinate system for a covariate slope: once the MH sampler has let η absorb the covariate-correlated part of the between-subject variation, the slope sees almost no gradient, and the two drift together — on the fluconazole model of #619 the renal gradient reached its bound at 0.

ferx now records such a value as a covariate mu-reference (all its thetas, one eta) and gives it the mu-referencing M-step in its general form. Holding each subject’s $\phi_i = g(A_i(\theta)) + \eta_i$ fixed, the group’s thetas are re-fitted to the population of individual values:

\[ \hat\theta = \arg\min_\theta \; \sum_i \tfrac12\,(\phi_i - g(A_i(\theta)))^\top \Omega^{-1} (\phi_i - g(A_i(\theta))) \;+\; \sum_i -\log p(y_i \mid \phi_i, \theta) \]

with $A_i$ the eta-free typical value at subject $i$’s covariates and $g$ the link (log or the logit-scale identity). The data term drops out — and the step is an exact nonlinear least-squares fit (Gauss–Newton) — when three things hold: every covariate the value reads is constant within each subject, no estimated theta of the group is read anywhere else in the model, and no other individual parameter reads the group’s eta. Otherwise the term stays live and the sum is minimised numerically (BOBYQA), still in the $\phi$-frozen coordinates. A covariate that varies within a subject leaves a $\theta$-dependence through the within-subject ratio $A_i(t)/A_i(t_0)$; a theta the model reuses — CL = (TVCL + TH_X*WT) * exp(ETA_CL) alongside V = TVV + TH_X — leaves one through that other parameter, since preserving $\phi_{CL}$ pins CL but not V; an eta the model reuses — the same CL alongside V = TVV * exp(0.5*ETA_CL) — leaves one through the re-centring, which holds $\phi_{CL}$ fixed but moves V. ferx detects both kinds of reuse and says so in a fit warning naming the theta or the eta. The new estimate is blended with the same SA step $\gamma_k$ as the single-anchor shift, and each η is re-centred by the change in its own mu so $\phi_i$ is unchanged. This is NONMEM’s MU_1 = LOG(THETA(1) + (CRCL-90)*THETA(2)) — a MU that is nonlinear in THETA — done automatically.

A covariate mu-reference is used by SAEM and IMP/IMPMAP only; mu_refs (inner-loop centring, suggest_start, reporting) is unchanged, so FOCE/FOCEI/Laplace fits are byte-identical with and without one. It is declined — with a warning naming the thetas, which then stay on the numerical M-step as before — when one of its thetas is also another eta’s single anchor, when the eta’s initial IIV is below 1e-3, when every theta of the group is FIXed, under a mixture model, and under mu_referencing = false. A group step also stands down for one iteration when the typical value is not finite for some subject at the current θ — an additive form such as TVCL + (CRCL-90)*TH_CRCL can go ≤ 0 for a low-covariate subject — and the fit reports how many iterations that happened on, since the thetas were then moved by the numerical M-step (or, under IMP/IMPMAP, left to it) rather than by the group. A value that reads TIME, MIXNUM, a second eta, or a local variable that is assigned conditionally is not recorded at all.

Scalar residual-error update

For a single, free additive or proportional residual SD in a single-population, uncensored Gaussian model whose free structural thetas are all log-mu-referenced, SAEM also averages the residual sufficient statistic instead of replacing the SD with the latest MCMC draw’s conditional optimum. For additive error it stores

\[ S_r \leftarrow (1-\gamma_k)S_r + \gamma_k\sum_{ij}(y_{ij}-f_{ij})^2, \qquad \sigma = \sqrt{S_r/n_{\mathrm{obs}}}. \]

For proportional error, the summand is divided by $f_{ij}^2$. This reduces otherwise persistent final-draw Monte Carlo noise in the residual SD. It is an internal exact M-step for that restricted likelihood, not a new fit option. Combined or multi-endpoint error, block_sigma, M3 censoring, IOV, iiv_on_ruv, FREM, custom residual magnitude, log-transform-both-sides, mixtures, fixed residual SDs, and models with a free numerical structural theta continue to use the general NLopt M-step.

The NLopt M-step (BOBYQA)

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.

Thetas the closed-form update cannot reach

Two model shapes leave a theta on the numerical M-step. The first is slower; the second can bias the estimate outright.

Non-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 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. (A typical value that reads several thetas, such as an additive covariate model, is not in this class since #619 — see covariate mu-references.) 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)), or, for a parameter bounded to $(0,1)$, the logit-scale form P = inv_logit(LOGIT_P + ETA_P).

A theta with no ETA at all rides a different, noisier update

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 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 each iteration’s MCMC η draw and assigns the maximiser. That is a valid stochastic-EM update — its fixed point is the marginal optimum whenever the E-step samples the right conditional — but it is noisier than the closed-form shift, it runs only every third exploration iteration, and on a poorly mixing chain it can settle away from the marginal optimum (see the FREM case below).

Before 0.4.x (#1415) this channel was broken in a way that read as “biased”: the NLopt solve behind it started from a first design a quarter of the bound range wide (×12–×26 of the current value on the default bounds) and stopped on a 1e-4 relative tolerance, so it returned a few percent of a step; and #1011 had then blended even that at 3 % per exploration M-step. The result was a theta that travelled a small fraction of the way to the optimum and stopped, on every model with such a theta — allometric and maturation exponents, covariate slopes, inter-compartmental clearances. Both are fixed: the solve starts local and converges, and the exploration cap is off by default (it is kept only for iiv_on_ruv models — see the next section).

Measured on the three ferx-testdata models of #1415, importance-sampled −2 log L at the final estimates (lower is better; the IS readout, not the FOCE approximation the fit YAML reports — see below):

model (no-ETA thetas)	before	after	reference
melphalan (`TVQ`, `CL_CRCL`)	2705.8	2684.0	NONMEM SAEM 2682.4
clofarabine (`TVQ`, `TVV2`, `THALF`)	17778.5	17764.3	NONMEM SAEM 17734.6
thiotepa (nine of eleven thetas)	6047.8	5827.3	IS at the FOCEI optimum 5823.6

On thiotepa every no-ETA theta had ended within a few percent of its initial value; after the fix they land where FOCEI lands. The in-repo anchor is covmuref_power with the weight covariate routed through an if block so the #619 group detection does not apply and TH_WT rides this channel: NONMEM SAEM 0.921, ferx 0.332 before (its 0.3 start) and 0.961 after.

Two structural remedies remain better than this channel where they apply:

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.
Write a covariate model in a form #619 can see — a covariate read from the data directly, not through a conditionally assigned local — so the slope joins the typical value’s mu-reference group.

Otherwise cross-check against FOCEI or IMPMAP, which do not depend on this channel.

Averaging a maximiser is not maximising the average

The numerical M-step re-maximises the frozen-η objective against one draw and then adopts that maximiser. A maximiser is a nonlinear function of the draw, so what the run converges to is E[θ*(η)] rather than the θ that maximises E[Q(θ, η)] — SAEM’s actual target. The difference is a Jensen-type bias that grows with the dispersion of the draws, which is why it gets worse when the E-step samples better (#1458).

It is invisible in the trace: the objective is fine, and only a comparison against a long-run reference or a FOCEI fit shows it. Two opt-in alternatives attack it: one removes the Jensen term by never forming a maximiser, the other shrinks it by averaging more draws into the one it forms. Both are off by default, and a fit that sets neither is bit-identical to before.

The two alternatives

mstep_solver = score_sa

Solves the score equation E[∇Q(x, η)] = 0 by stochastic approximation, preconditioned by the SA-averaged expected information (Kuhn & Lavielle 2005):

I_k = (1 − γ_k)·I_{k−1} + γ_k·I(x_k, η_k)
x_{k+1} = x_k − γ_k·(I_k + λ·diag I_k)⁻¹·∇Q(x_k, η_k)

A fixed point needs E[∇Q(x*, η)] = 0, which is the stationarity condition of the averaged objective, so no Jensen term arises. σ takes the same Newton direction as every other coordinate, but the step gives it a target that damp_mstep_sigma_variance then moves it a γ_σ = min(γ, 0.2, γ_mstep) fraction of the way to, on the variance scale — the same #1445 policy the derivative-free arms use, not a second copy of it. That is a smaller gain on the step, not the blend stacked on top of it (which really would move σ at γ²). One gradient sweep is cheaper than the mstep_maxiter·(n+1) population passes the derivative-free solve costs, so this M-step runs every iteration rather than every third during exploration.

mstep_draws = K

Averages the derivative-free M-step objective over this iteration’s η draw and the previous K − 1 iterations’ draws, so the maximiser NLopt returns is the maximiser of an average. The bias falls roughly as 1/K — it shrinks, it does not go away — at K times the M-step’s objective evaluations, and only to the extent that consecutive draws decorrelate, which under the default sticky random-walk E-step they largely do not. Ignored under score_sa.

Both are restricted to the plain Gaussian residual scope the expected information has a closed form on: no IOV, no mixture, no M3, no time-to-event endpoint, no block_sigma residual correlation, no [covariate_nn] θ, no residual magnitude, no FREM. A model outside it keeps the historical solver and says so by name in the fit’s warnings.

Measured

On the busulfan benchmark of #1458 (ferx-testdata/busulfan_dominic, 600 subjects, analytic 2-cpt IV, block_omega(ETA_CL, ETA_V1) + omega ETA_V2, combined error, TVQ with no ETA), 6 seeds at the default 150/250/20 schedule, on the importance-sampled −2 log L — the neutral readout, since a SAEM fit has no objective of its own. B6qfix is the same model with TVQ held FIX, which removes this channel and nothing else:

	`TVQ` free	`TVQ` FIX	what a free no-ETA θ costs
default (`bobyqa`)	60044.87 ± 3.29	60035.12 ± 2.08	+9.75
`mstep_solver = score_sa`	60034.59 ± 5.27	60031.22 ± 3.35	+3.37
`mstep_draws = 2`	60043.87 ± 3.77	60034.62 ± 3.45	+9.25
`mstep_draws = 3`	60045.61 ± 2.36	60033.99 ± 3.50	+11.62

score_sa removes about two-thirds of the penalty, at the same cost (6.40 against 6.47 CPU-seconds: one gradient sweep every iteration in place of a derivative-free solve every third). mstep_draws does not move it at all here and costs 53 % (K = 2) to 89 % (K = 3) more CPU — consecutive draws of the default random-walk E-step are too correlated for a two- or three-draw average to be worth anything, which is the regime #1458 predicted it would fail in. It is kept because that argument is about the kernel, not about the estimator, and a near-i.i.d. kernel changes it.

On the 100-subject subset the three arms are a tie (7531.37 ± 0.37 default against 7531.14 ± 0.93 for score_sa), with score_sa’s seed spread 2.5× wider.

Why `score_sa` is not the default

It was a candidate in #1449 and it is not a close call on either side, which is why the answer is worth recording.

For. Re-measured at today’s shipped defaults (n_mh_steps = auto) over all eight SAEM benchmarks, 6 seeds each, single-threaded, it is tie-or-better than bobyqa on the importance-sampled −2 log L on every one — best on busulfan (−5.14), busulfan-TVQ-FIX (−5.58) and the thiotepa ODE bench (−11.07) — and it costs 25–54 % less CPU (cefepime 5.36 → 2.45 s, thiotepa 215 → 94 s), because one gradient sweep per iteration replaces a derivative-free solve whose cost is mstep_maxiter·(n+1) population passes.

Against — the σ half, fixed in #1480. As #1458 shipped it, score_sa stepped σ at the shared γ and in packed (log σ) units, which is #1445’s two defects at once: during exploration γ = 1, so that is a full Newton step to one draw’s score root with no Robbins-Monro averaging, and a log-scale blend is a geometric mean of the σ sequence. On the sparse combined-error fixture of #1445 — 300 subjects with a median of one observation each, truth combined(prop = 0.13, add = 1.8) — it scattered the minority additive σ. Measured seeds 1 / 2 / 3, ci-test, aarch64:

	`ADD_ERR` (truth 1.8)	cross-seed sd
`bobyqa` (the default)	1.33 / 1.83 / 2.17	0.42
`score_sa` as #1458 shipped it	0.84 / 0.52 / 1.13	0.31
`score_sa` with the γ_σ schedule only	0.87 / 0.99 / 1.20	0.17
`score_sa` after #1480	1.90 / 1.86 / 2.11	0.14

The third row is why the fix is both halves: the schedule alone still leaves seed 1 at 0.87, under the factor-of-two gate of 0.9. The objective could not see any of it — the FOCE readout differs by under 2 units across all those fits, and the benchmark suite’s −2 log L prefers score_sa throughout. This is the same channel #1445 is about, a variance component identified only where σ_add and σ_prop·f are comparable, and it is why tests/saem_combined_error.rs scores it against simulation truth rather than against another estimator.

Against — the θ half, still open (#1480). On the #619 covmuref_power anchor, whose TH_WT carries no ETA and therefore moves only through this channel, score_sa at the default 150/250 schedule reaches 0.746 against NONMEM SAEM’s 0.921 — outside the 0.12 window the default solver satisfies. It is slow rather than wrong: at 300/700 it lands at 0.9845 (|Δ| 0.063) while bobyqa at that same schedule overshoots to 1.1485 (|Δ| 0.227). So the window is a property of the estimator and the schedule together, and score_sa needs a longer run on a θ that starts far from its optimum. tests/saem_covariate_mu_ref.rs keeps the anchor at 300/700 live and gates the default-schedule arm on the open issue.

A default has to be safe for the parameter the likelihood barely constrains and has to converge in the schedule it ships with, so score_sa stays opt-in. It is worth reaching for on a model with a free θ that has no ETA, with a longer schedule than the default if that θ starts far away.

One further thing #1449 measured: score_sa’s score/information sweep does not reuse the cached event schedule, so on a model whose [individual_parameters] read TIME or a time-varying covariate it rebuilds two schedules per subject per iteration. That is an unrealised saving on top of the CPU numbers above, not a regression against them.

Two things the numbers do not say. First, d_final against a long-run reference gets worse under score_sa (TVQ moves from −0.04 to −0.18 log units away from it) — but that reference is the same estimator run longer, so it carries the same bias, and it cannot arbitrate a change that removes it. Second, a FOCEI fit of the same model is not the oracle either: its own importance-sampled −2 log L is 60046.43, worse than every SAEM arm above, so it locates TVQ (7.60) no better than they do.

Neither is a substitute for attaching an ETA

A θ that carries a random effect gets the exact closed-form Robbins-Monro shift and never enters this channel at all. Both options above make the numerical channel unbiased in the limit; the closed form is exact at every iteration. Where the model can be written so #619’s covariate mu-references apply, that remains the better fix.

The exploration cap (`mstep_damping`) is a hold, not a cure

mstep_damping blends the numerical M-step result in as θ ← θ + γ_θ·(θ* − θ) during exploration, with γ_θ capped at the value you set, and averages it at γ = 1/(k−k1) in convergence. It is off by default (1.0, plain assignment in both phases) for every model except one with iiv_on_ruv, which keeps the 0.03 cap — the reason is measured below:

[fit_options]
  method        = saem
  mstep_damping = 0.03   # opt in: hold the no-ETA thetas near their start

It was introduced by #1011 at a default of 0.03 after a FREM iiv_on_ruv model (475 subjects, 12 etas) whose absorption-split fraction TVFRD1 (D1 = MAT*(1-TVFRD1), KA = 1/(MAT*TVFRD1), no eta of its own) drifted from 0.383 to 0.039 under the undamped update, while IMP (0.311), IMPMAP (0.318) and NONMEM IMP (0.394) agree. The capped fit ended at 0.29–0.36, which read as “landing on the marginal optimum”. #1415 measured what the cap does, with the M-step solver repaired:

It holds the theta where it starts. On the same reprex: 0.383 → 0.361, and from a deliberately wrong 0.2 start, 0.2 → 0.191. Its good answer was the model’s initial estimate, which that auto-generated model had taken from a FOCEI fit.
It freezes every other no-ETA theta the same way, because the cap divides the number of EM steps the exploration phase amounts to (50 M-steps × cap), and a covariate slope confounded with an eta needs all of them. covmuref_power TH_WT (NONMEM 0.921, start 0.3): cap 0.03 → 0.323, 0.1 → 0.340, 0.3 → 0.403, off → 0.961. Thiotepa (IS −2 log L, optimum 5824): 0.03 → 5928, 0.1 → 5897, 0.3 → 5860, off → 5827.
The convergence-phase Robbins-Monro average hurts too. γ = 1/(k−k1) is the right average for a statistic that is stationary around the optimum, but the sequence of η-frozen maximisers is an EM trajectory, so averaging it weights the early, still-crawling iterates most (TH_WT 0.637 with the average, 0.961 without).

The FREM drift itself is real and is not fixed: undamped, that reprex settles near TVFRD1 0.18 from either start on the default schedule, keeps drifting on a 3× longer one (0.03), and 200 MH steps per iteration do not change it — so it is not plain E-step under-mixing. The importance-sampled objective cannot rank the two points on a 12-eta model (ESS/K = 0.001), but the FOCEI objective can: the drifted fit is 47 units worse than the held one (7492.7 against 7445.3), so on this model the hold is the better answer.

The drift is the iiv_on_ruv coupling, not FREM and not the no-ETA channel as such: the same model with its iiv_on_ruv line and ETA_RUV omega removed lands undamped at TVFRD1 0.414, TVMAT 2.680, TVV 137.6 — on NONMEM IMP’s 0.394 / 2.680 / 133.8 — where the 0.03 cap leaves it at 0.342 / 2.228 / 147.5. Under iiv_on_ruv the E-step is modified (η_RUV is re-centred into σ every iteration) and σ shares the numerical M-step with the no-ETA thetas; that is where the pathology lives, and it needs its own fix (#1421). So the default is keyed to that shape: a model with iiv_on_ruv keeps the 0.03 cap, every other model runs undamped, and a value you set wins either way. FIX is the honest spelling of the same hold when the value is known.

The option only acts on the channel it is about: 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 is left on the undamped update, and setting mstep_damping on a model where no theta trips the gate warns that it has no effect rather than accepting it silently. Sigma takes the same step as theta when the cap applies, because the numerical M-step is a single NLopt problem over [θ; σ] and a partial step on one component alone would be inconsistent. Under mu_referencing = false a theta written in a mu-referenceable form is not this channel either, so such a fit keeps the undamped update and the option reports no effect.

Mixture models are excluded outright: a MIXNUM-switched typical value goes through the same numerical M-step, but its class typical values must separate from a common start, and damping that excursion stalls the separation. For the same reason the mixture arm also keeps the pre-#1415 M-step solver configuration: its class typical values come from a hard class draw, and the NONMEM-anchored mixture seed sweep was calibrated on the partial step that configuration returns — a converged per-class solve moves TVCL2 from the NONMEM MLE 2.84 to 3.23. An exact M-step on a hard class draw is a different estimator, so that channel changes only with its own anchor.

Comparing a SAEM fit to FOCEI: read the importance-sampled objective

The ofv a SAEM fit reports is the FOCE approximation at the final estimates — SAEM has no objective of its own — and until 0.4.x (#1415) a method chain such as method = [saem, imp] computed it without the η–ε interaction, while method = saem alone and every FOCEI fit computed it with. On the thiotepa model that is a 175-unit difference at the same estimates, and it was the bulk of the “SAEM is 341 units worse than FOCEI” that #1415 reported. A chain now keeps the interaction flag unless its final stage is foce, as the single-method form always did. When you compare estimators, compare the importance-sampled −2 log L that an imp_eval_only stage prints (IS done. −2 log L = …) or FitResult::importance_sampling, evaluated at each fit’s estimates.

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

The per-subject step scales $\delta_i$ (MH) and leapfrog step sizes (HMC) are adapted toward a target acceptance rate. Which rule does the adapting is set by [fit_options] scale_adaptation.

Target acceptance rates: 40% for the primary block kernel (MH or HMC’s own 65% when n_leapfrog > 0), 44% for the componentwise kernel. The block kernel and the componentwise kernel carry independent per-subject scales, and each η’s componentwise scale is adapted separately.

`scale_adaptation = interval`

Every adapt_interval iterations, from that window’s acceptance rate:

above target: increase $\delta_i$ by 10% (up to 5.0)
at or below target: decrease $\delta_i$ by 10% (down to 0.01, or 1e-6 for a componentwise scale)

Know its reach. The factor is fixed, so the total movement available is 1.1^n up or 0.9^n down over the n times the rule fires — and on the default 150 + 250 = 400 iterations at adapt_interval = 50 that is 8 times, i.e. at most ≈2.1× up or ≈0.43× down however far from target the chain is. A model whose optimal step is an order of magnitude below the 0.3 starting value never reaches it. examples/warfarin_saem.ferx is such a model: on the default settings its combined acceptance collapses to 1.6% by iteration 26 and sits at 2–4% for the remaining ~375 iterations against the 40% target (issue #1444).

`scale_adaptation = robbins_monro` (default)

A Robbins–Monro step taken every iteration, from that iteration’s own rate:

\[\log \delta_i \leftarrow \log \delta_i + c \, k^{-0.6} \, (\hat{a}_k - a^\ast)\]

with $c = 1$, $k$ the 1-based iteration, $\hat{a}_k$ the realised rate and $a^\ast$ the target. The correction is proportional to the discrepancy, so the reach is unlimited, and the $k^{-0.6}$ decay is the diminishing-adaptation condition that keeps an adapting chain valid (Roberts & Rosenthal 2007). The clamps are the same as the interval rule’s.

It is the default since #1449, and it was not always a free win. What ferx measures directly is the acceptance rate, and on examples/warfarin_saem.ferx at its shipped settings (150 + 250 iterations, n_mh_steps = 20, seed 12345) the combined rate over the last 100 iterations goes from 3.90% under interval to 42.05% under robbins_monro, against the 40% target. That much is reproduced by tests/saem_scale_adaptation.rs.

Whether a rate at target makes the estimates better is a separate question, and one this page does not answer from its own runs. The SAEM benchmarking in issue #1444 reported, over 6 seeds on real datasets, that Robbins–Monro roughly halves warfarin’s final-estimate error against a long-run reference, changes little on models whose baseline scale is already near target, and regressed pembrolizumab (303 subjects, 4 η) — a model whose acceptance was already above target — by +4.7 on the importance-sampled −2 log L. That regression is what kept the rule opt-in.

It is gone at today’s defaults, and the reason is worth knowing: those measurements were taken at a fixed n_mh_steps = 20, and #1468 made the proposal count auto, which resolves to 6 on that model. Re-measured at the shipped defaults over 6 seeds (#1449), Robbins–Monro improves pembrolizumab by −4.7, and the only bench it costs anything on is vancomycin (+0.98 against a 0.94 cross-seed sd) — which the dead band below removes. #1444’s own data already pointed this way: at n_mh_steps = 8 it won on that bench too. These numbers have no NONMEM counterpart (NONMEM’s SAEM exposes no comparable step-scale knob).

scale_adaptation = interval restores the historical rule exactly.

The κ (IOV) scales stay on the interval rule under both settings: the per-occasion kernel makes one proposal per occasion per iteration, so a single iteration’s rate for one subject is a ratio of small integers, and the Robbins–Monro arm has not been measured there.

`scale_deadband` — only correct a chain that is off target

[fit_options] scale_deadband = <lo>,<hi>, default 0.15,0.60, makes the Robbins–Monro step conditional: it fires only when that iteration’s acceptance is outside [lo, hi], and inside the band the scale is left bit-for-bit alone (there is no interval bump underneath it — the interval rule is off for the whole run under robbins_monro). The band is one absolute acceptance window shared by the block kernel (target 40%) and each componentwise coordinate (44%). It is ignored under scale_adaptation = interval, which has no per-iteration step to skip.

It exists for the failure the previous section names: the models Robbins–Monro was written for sit at 2–4%, and the one it regressed was already above target and got pulled down onto it.

What it does, measured, is not quite what it was built for. Over the 8-benchmark, 6-seed suite of #1449 the band turns the unbanded rule’s one loss into a win (vancomycin +0.98 → −0.13 on the importance-sampled −2 log L, against a 0.94 cross-seed sd) and keeps most of its gains (pembrolizumab −4.69 → −3.74, warfarin’s rescue intact: tail acceptance 0.041 under interval, 0.422 unbanded, 0.391 banded). But it does not leave a near-target chain alone: the realised tail acceptance is lower under the band than without it on every bench measured, because of the fixed-point shift described below. The improvement is real; the mechanism is not “the controller declined to act”.

Two things to know before setting one:

A wide band freezes the scales; it does not disable the dead band. The rule fires only outside the band, so a band wider than the spread of the per-iteration rate means no step ever runs. scale_deadband = none (or an empty band, lo == hi) is the spelling that restores plain Robbins–Monro, and a band covering all of [0, 1] is rejected with that message rather than accepted.
Setting it from Rust is validated too. FitOptions::saem_scale_deadband is public, so a band can reach the sampler without passing the model-file parser. An empty, reversed, out-of-range or everything-covering band is dropped there — the fit runs the unbanded Robbins-Monro rule and says so in its warnings — so the two front ends cannot produce different numerical behaviour for the same option, and neither can freeze the step scales.
A band moves the fixed point. A per-iteration rate is n_accepted / n_proposals, so it is quantised to multiples of 1/n_mh_steps (block) or 1/n_cw_sweeps (componentwise), and the band censors the small deviations while keeping the tails. Unless the band is symmetric about the target and the rate’s distribution is symmetric, the surviving steps do not cancel at the target: the scale settles where the censored deviations balance, which is not where the unbanded rule settles. With the shipped band and n_mh_steps = 6 that equilibrium is near 33% rather than 40%.

Acceptance-rate diagnostics

Two warnings land in FitResult.warnings (never on stderr), most severe first:

Not mixing — the cumulative post-burn-in acceptance is below 1%. The sampled ETAs barely moved, so the M-step ran on degenerate sufficient statistics and Ω/σ are unreliable. The realised tail rate is also available as FitResult::saem_mh_accept_tail, so a script can check mixing without matching on the warning text.
Persistently far from target — the acceptance over the last 100 post-burn-in iterations is outside 10–80%. The band is deliberately wide: a fit at 0.25 or 0.60 is unremarkable, and this is not a range of “good” rates but the range outside which the step scales have demonstrably failed to find their target. The message carries the realised tail rate, the target, the window length and the cumulative rate, so you do not have to re-run with optimizer_trace to see how far off it was. It is the tail rather than the whole-run average that is reported, because a run that mixes badly early and recovers averages to a number describing neither half; and it stays quiet for runs with fewer than 25 post-burn-in iterations, where the realised rate describes the starting scale rather than a failure to reach target.

Post-SAEM Finalization

After the SAEM iterations complete:

EBE Refinement: Run the standard FOCE inner loop (BFGS optimization) warm-started from the SAEM ETAs to obtain final empirical Bayes estimates
Combined-error additive-collapse report: for combined(PROP, ADD) residual-error models, when SAEM has left a free additive component ADD at or below 1e-3 — the detection band just above its exp(-8) ≈ 3.35e-4 optimizer lower bound — the fit carries a warning saying so. It is a report, not a repair — ferx once ran a FOCEI marginal-likelihood polish here and adopted it if it lowered the OFV, and that silent second fit was removed in favour of the warning. The channel it was papering over is fixed at source: the σ M-step now takes a damped σ step instead of assigning a single draw’s maximiser (issue #1445), so the collapse it guarded against no longer happens on a combined() model whose additive term is identified. The warning remains, because an additive term that genuinely is not identified in the data will still land on the bound and you should see that: check it yourself (RSE, σ-correlation, profile likelihood), and on sparse data prefer the importance-sampling −2LL over the Laplace OFV when judging whether the additive term is real.
FOCE OFV: Compute the objective function using the FOCE/Laplace approximation, so AIC and BIC are directly comparable with FOCE results
Covariance Step: Optionally compute standard errors via finite-difference Hessian (same method as FOCE)
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.)

Validation of the logit closed-form M-step (#918)

Parallel dual first-order absorption (fast KA1 / slow KA2) with the pathway fraction carrying logit-scale IIV — FR1 = 1/(1 + exp(-(LOGIT_FR1 + ETA_FR1))). 60 subjects × 16 samples, simulated from the model itself (nonmem_anchor/simulate_logit_fraction_data.py); the NONMEM side writes the same relationship with explicit mu syntax (MU_3 = THETA(3), nonmem_anchor/logit_fraction_saem.ctl, METHOD=SAEM). Both start at FR1 = 0.4, the mirror image of the truth.

Quantity	truth	ferx SAEM (before #918)	ferx SAEM	NONMEM SAEM
`LOGIT_FR1`	0.4055	0.2462	0.3793	0.4219
`FR1` = `inv_logit(LOGIT_FR1)`	0.600	0.561	0.594	0.604
`ω²(ETA_FR1)` (logit scale)	0.25	0.326	0.186	0.183
`CL`	5.0	5.35	5.31	5.29
`V`	50.0	51.4	54.6	55.7

Before the fix the fraction’s typical value went through the numeric M-step and stalled short of the truth, with ω²(ETA_FR1) inflating to absorb the bias — the same θ / η-mean confounding the closed-form step exists to remove. The slow-gated tests/saem_logit_mu_ref.rs pins the comparison.

Validation of the covariate mu-reference M-step (#619)

Two matched fixtures, 60 subjects × 7 samples each, simulated from the model under test (nonmem_anchor/simulate_covmuref_data.py; covariates constant within each subject, so the exact Gauss–Newton engine runs). The NONMEM side is METHOD=SAEM with the mu written out: the nonlinear MU_1 = LOG(THETA(1) + (CRCL-90)*THETA(2)) for the additive renal gradient, the linear MU_1 = LOG(THETA(1)) + THETA(2)*LOG(WT/70) for the allometric exponent (nonmem_anchor/covmuref_{additive,power}_saem.ctl). Both engines start at the same values, deliberately away from the truth.

Quantity	truth	start	ferx SAEM (before #619)	ferx SAEM	NONMEM SAEM
additive: `TVCL`	5.0	4.0	3.990	4.921	4.917
`TH_CRCL`	0.05	0.02	0.0198	0.04179	0.04178
`TVV`	50	40	48.09	48.09	48.09
`ω²(ETA_CL)`	0.09	0.1	0.083	0.0403	0.0401
`ω²(ETA_V)`	0.04	0.1	0.0360	0.0358	0.0360
power: `TVCL`	5.0	4.0	4.919	4.787	4.788
`TH_WT`	0.75	0.3	0.332	0.9221	0.9213
`TVV`	50	40	49.51	49.50	49.52
`ω²(ETA_CL)`	0.09	0.1	0.0975	0.0701	0.0697
`ω²(ETA_V)`	0.04	0.1	0.0459	0.0460	0.0459

Before #619 the covariate theta never left its start in either shape — the eta-frozen numerical M-step has no leverage on it once the sampled η have absorbed the covariate structure, and the TVCL/TH_CRCL pair reported a correlation of 1.00 with RSEs in the tens of thousands of percent. With the group step ferx and NONMEM agree to three or four figures on every parameter (TH_CRCL to 8e-6), including the shape NONMEM handles most efficiently. tests/saem_covariate_mu_ref.rs (slow-gated) pins both comparisons and an IMP twin.

The time-varying engine has no NONMEM SAEM comparator (NONMEM’s MU must be constant within a subject). On the fluconazole model that motivated #619 (CL = (TVCL_NR + (CKDCRCYS-90)*THETA_CLCKD) * exp(ETA_CL), CKDCRCYS changing within 27 of 31 subjects, SAEM seed 123) the renal gradient moves from a frozen 1.96 (start 2.0) to 1.18 and TVCL_NR from 153.9 (start 150) to 109.0, against NONMEM FOCE’s 1.41 and 122; the remaining distance to NONMEM’s OFV on that dataset predates this change and is tracked with the FOCEI comparison of the same model.

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    = auto     # block-MH steps/subject/iteration, or `auto` (see above)
  n_leapfrog    = 0        # Set > 0 (e.g. 3) to use HMC instead of MH
  adapt_interval = 50      # Step-size adaptation frequency (interval rule)
  scale_adaptation = interval  # or robbins_monro
  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 above what auto resolves to (e.g. 30-50) for better mixing in the E-step on hard surfaces. The rule caps at 20, the count calibrated to escape the basin trap observed on Emax PKPD with stressful initial values (see below), and it reaches that cap only at 8 observations per subject per η; an ODE-with-Form-C readout may need more proposals than either 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.

This benchmark has 16 observations per subject against 2 η, so the auto rule resolves it to the same 20: the count that closed this trap is unchanged on the model it was calibrated on (#1459). Re-measuring it on e7343462 at 6 seeds, the boundary collapse itself no longer reproduces at any count from 2 to 20 — E0 lands at 1.83–2.02 against a truth of 2 and GAMMA at 0.89–1.08, never at a bound — but what 20 still buys is reproducibility along the EMAX/EC50 ridge: final-estimate distance 0.029 and cross-seed stability 0.044 at 20, against 0.076–0.173 and 0.090–0.160 at 4–8. That ridge is a property of the data, not of the sampler (FOCEI on the same dataset lands at EMAX 20.0 / EC50 1.78 where SAEM lands at 53.1 / 7.28), so report both endpoints and check several seeds on any Emax/Hill fit.

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

Try n_mh_steps = 50 (above the auto cap of 20)
Warm-start from a FOCEI fit (method = [focei, saem])
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.

The additive term of a `combined()` error model

Before ferx #1445 the additive component of combined(PROP, ADD) could be driven onto its optimizer lower bound on sparse data where NONMEM’s own SAEM estimates a clearly non-zero value, and it got worse with more iterations. The cause was the σ half of the numerical M-step being assigned outright rather than averaged: the additive term of a combined model is often a minority variance component, identified only where $\sigma_{\text{add}}$ and $\sigma_{\text{prop}}\cdot f$ are comparable, and the conditional maximiser of a weakly identified variance component from a single η draw sits near zero a large fraction of the time. The damped σ step makes the estimate the average of that sequence rather than its last draw.

Anchored against NONMEM 7.5.1 METHOD=SAEM on the same 300-subject, median-one-observation dataset (nonmem_anchor/saem_sparse_combined_saem.ctl and data/saem_sparse_combined.csv, which is the fixture tests/saem_combined_error.rs scores; results in nonmem_anchor/results/saem_sparse_combined_saem.ext, read from the .ext because a SAEM+IMP chain truncates its .lst):

	additive SD (truth 1.8)	proportional SD (truth 0.13)
NONMEM SAEM, `ISAMPLE=10`	1.5425	0.1370
NONMEM SAEM, `ISAMPLE=2` (its default)	1.6463	0.1491
ferx SAEM after #1445, 8 seeds	1.42 [1.11, 1.75]	0.139
ferx SAEM before #1445, 3 seeds	0.011 / 0.0025 / 0.011	0.135
ferx FOCEI	3.220	0.1200

ferx SAEM now agrees with NONMEM SAEM on the additive term to 8% of the simulation truth and on the proportional term to 1.5%, where before it sat four orders of magnitude below it. (NONMEM’s ISAMPLE=2 default collapses the two weakly identified Ω diagonals on this design — 8.8e-4 and 0.032 against truths 0.15 and 0.69 — which is why the ISAMPLE=10 arm is the headline; the σ estimate moves only from 1.65 to 1.54 between them.) On a real dataset, ferx-testdata/cefepime_jordan run64, NONMEM FOCEI estimates the additive RUV at 1.550 (nm/run64.ext, THETA(6)); ferx SAEM reports 1.16 [0.97, 1.34] after the change and ~1e-3 before it.

Note what the objective does not say here: on the simulated fixture the importance-sampled −2 log L is 3183.1 after against 3183.4 before, a 0.2 difference that is pure noise. That is not a defect in the comparison, it is the whole problem restated — a minority variance component is weakly identified, the likelihood is nearly flat in it, and an estimator that reports a single draw of a flat direction will wander to the boundary while barely moving the objective. Rank the two by the parameter against a known truth and an external engine, not by the likelihood they are nearly tied on.

The price is that σ now lags: it is a Robbins-Monro average over the convergence phase, so a run that is too short reports a σ still on its way down from the starting guess. If σ has clearly not settled, lengthen the run (n_exploration first, so averaging starts from an equilibrated chain) rather than reading the last value. If ADD still lands on its bound, it is genuinely not identified by the data — check it (RSE, σ-correlation, profile likelihood) rather than assuming the estimator.

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, or scale_adaptation = robbins_monro when the acceptance-rate warning says the scales never reached their target

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.