Mixture Models

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

A mixture model treats each subject as belonging to one of \(K\) unobserved subpopulations (“classes”). Class \(k\) has its own set of class-specific typical values (and, optionally, its own \(\Omega\) / \(\Sigma\)), and a subject is assigned to a class only probabilistically, through covariate-dependent mixing probabilities. This is the ferx analogue of NONMEM’s $MIXTURE.

Model file

Declare the number of classes and the mixing rule in a [mixture] block, and select class-specific typical values with the reserved read-only MIXNUM index (1..=K) inside [individual_parameters]:

[parameters]
  theta TVCL1(1.0); theta TVCL2(3.0)   # class-specific clearance
  theta TVV(10.0)
  theta MIXL(0.0)                       # mixing logit intercept
  omega ETA_CL ~ 0.09
  sigma EPS ~ 0.04

[mixture]
  nsub = 2
  logit(1) = MIXL                       # p(1) = σ(logit);  p(2) = 1 − p(1)
  # omega(2) ETA_CL ~ 0.15              # optional per-class Ω override
  # sigma(2) EPS ~ 0.05                 # optional per-class Σ override

[individual_parameters]
  CL = if (MIXNUM == 1) TVCL1 * exp(ETA_CL) else TVCL2 * exp(ETA_CL)
  V  = TVV

[structural_model]
  pk one_cpt_iv(cl=CL, v=V)

[error_model]
  DV ~ proportional(EPS)

The mixing rule may be given as a logit (logit(k) = …) or a probability (p(k) = …) expression over thetas and covariates; the remaining class’s probability is the complement. Mixing expressions must not depend on any \(\eta\).

The two forms differ only in how \(p_{ik}\) is formed: logit(k) is passed through a softmax (reference class \(K\) has logit 0), while p(k) is used as the probability directly and the class weights are renormalised to sum to one. A constant p(1) = 0.3 is therefore identical to logit(1) = \ln(0.3/0.7). The softmax (logit) form has an analytic posterior-weighted outer gradient; the p(k) form does not, so a p(k) fit uses the finite-difference gradient.

Objective

Each subject’s marginal likelihood is the covariate-weighted mixture of its class-conditional likelihoods, and the population objective is the numerically stable log-sum-exp of the per-class FOCE/FOCEI negative log-likelihoods:

\[ p_{ik} = \operatorname{softmax}_k\big(g_k(\theta,\, \text{cov}_i)\big), \qquad L_i = \sum_{k=1}^{K} p_{ik}\, L_{ik}, \]

\[ \text{OFV} = -2 \sum_{i=1}^{N} \log \sum_{k=1}^{K} p_{ik}\, \exp(-\text{nll}_{ik}). \]

nll_ik is the ordinary per-subject FOCE/FOCEI negative log-likelihood evaluated with class-\(k\) parameters and a separate empirical Bayes solve per (subject × class). The mixing-logit coefficients are ordinary thetas and any per-class \(\Omega\)/\(\Sigma\) are estimated jointly. Estimation defaults to the derivative-free BOBYQA outer optimizer; an analytic posterior-weighted outer gradient is available, so a user-selected NLopt gradient optimizer (SLSQP / L-BFGS / MMA) is honoured, with an automatic finite-difference fallback for models outside analytic scope (e.g. MIXNUM-branched typical values).

Output

For a mixture fit, sdtab gains two per-subject diagnostics, matching the NONMEM reserved variables:

• MIXEST — the most-probable class \(\operatorname{argmax}_k \text{PMIX}_{ik}\) (1-based).
• PMIX_1 … PMIX_K — the posterior class-membership probabilities \(\text{PMIX}_{ik} \propto p_{ik}\exp(-\text{nll}_{ik})\), normalised over \(k\).

Every other per-subject diagnostic is evaluated in that subject’s MIXEST class: the reported \(\hat\eta\) are the winning class’s EBEs, and IPRED, PRED, IWRES, CWRES, EBE_OFV, and any [derived] / [output] column that branches on MIXNUM are all computed with MIXNUM set to that class. So a subject the fit assigns to class 2 gets class-2 typical values throughout its sdtab row, not the class-1 default.

Validation against NONMEM

A two-class clearance mixture (tests/nonmem/mixture_iv.ctl and the slow-gated tests/mixture_nonmem.rs) was fit by both engines on the same 30-subject IV dataset — two subpopulations with CL ≈ 1 and CL ≈ 3, a constant mixing fraction, and \(\Omega\)/\(\Sigma\) fixed at the data-generating values (METHOD=1 INTER in NONMEM; FOCEI in ferx). NONMEM’s minimization was successful.

Quantity	NONMEM 7.5.1	ferx
OFV (without constant)	301.582	301.600
TVCL1 (class 1)	0.981	≈ 0.98
TVCL2 (class 2)	2.842	≈ 2.84
V	9.979	≈ 9.98
mixing fraction p(1)	0.472	≈ 0.47
MIXEST classification	—	30 / 30 subjects agree

The objectives agree to ~0.02 units and the per-subject MIXEST class assignment is identical for every subject.

Standard errors

A mixture fit runs the covariance step like any other fit: the finite-difference Hessian is built on the mixture objective (the K-fold log-sum-exp above), and the standard errors are \(\sqrt{\operatorname{diag}(R^{-1})}\). The mixing-fraction SE is reported on the scale the mixing form is written in — the coefficients of a logit(k) = … form carry a logit-scale SE, a p(k) = … probability an SE directly on the probability scale — because those coefficients are ordinary thetas.

Cross-checked against $COVARIANCE MATRIX=R on the same model (tests/nonmem/mixture_iv_cov.ctl); MATRIX=R is the pure Hessian-inverse estimator ferx computes (NONMEM’s default sandwich \(R^{-1}SR^{-1}\) differs by the usual 10–25 % between the two estimators):

SE	NONMEM (MATRIX=R)	ferx
TVCL1	0.1056	0.1067
TVCL2	0.2543	0.2559
V	0.1861	0.1860
p(1) (delta-method)	0.1026	0.1030

All four agree to < 1 %. The p(1) row compares NONMEM’s direct \(P(1)=\theta_4\) SE to ferx’s logit-scale SE mapped through \(\operatorname{SE}(p) = p(1-p)\operatorname{SE}(\text{logit})\).

Inter-occasion variability under a mixture

A mixture may carry inter-occasion variability (a kappa on a structural parameter) — the two features compose (#985). Each class’s per-subject inner solve estimates the per-occasion κ̂ under that class’s typical values, and the FOCE/FOCEI marginal L_ik is the κ-augmented occasion likelihood, so the same L_i = Σ_k p_ik L_ik mixture is formed over IOV-aware class likelihoods. The IOV Ω is shared across classes (it sits in the base parameter block, ahead of any per-class Ω/Σ override), matching NONMEM’s $OMEGA BLOCK(1) … SAME occasion structure. The analytic outer gradient does not yet emit κ-slot derivatives for a mixture, so an IOV mixture optimises against a central-difference outer gradient; the covariance step runs as usual on the K-fold mixture objective.

NONMEM cross-check (tests/nonmem/mixture_iv_iov.ctl)

Two-class bimodal-CL model with a per-occasion κ on CL (30 subjects × 2 occasions), IIV/IOV Ω and Σ FIXed at the data-generating values in both engines, METHOD=1 INTER:

Quantity	NONMEM 7.5.1	ferx
OFV (without constant)	470.814	470.930
TVCL1 (class 1)	1.057	1.069
TVCL2 (class 2)	3.330	3.275
V	10.259	10.253
mixing fraction p(1)	0.526	0.506
MIXEST classification	—	30 / 30 subjects agree

The objectives agree to ~0.12 units and every subject’s most-probable class matches. Standard errors are not cross-checked for this model: NONMEM’s covariance step terminates with rounding errors and substitutes the R matrix (R MATRIX SUBSTITUTED: YES, no printable SE block) — a known NONMEM fragility for mixture + IOV covariance — whereas ferx’s covariance step completes. The mixture SE machinery itself is anchored to $COVARIANCE MATRIX=R in the non-IOV standard-errors comparison above.

SAEM under a mixture

Mixtures can also be estimated with SAEM (method = saem, #985). Where FOCE/FOCEI marginalise the latent class analytically (the log-sum-exp above), SAEM samples it: each E-step draws the class indicator \(z_i \sim \operatorname{Categorical}(\mathrm{PMIX}_i)\) from the current per-subject posterior \(\mathrm{PMIX}_{ik} \propto p_{ik}\,e^{-\mathrm{nll}_{ik}}\), then runs the usual η-MCMC within the drawn class (its Ω/Σ and its MIXNUM branch). The M-step partitions the sufficient statistics by sampled class:

• Class-switched typical values (if MIXNUM == k …) are estimated from each class’s own members — the ordinary θ/σ M-step run under each subject’s class gives \(\partial\text{pred}/\partial\text{TVCL}_2 = 0\) for a class-1 subject, so no special handling is needed. They are not mu-referenced under SAEM; see Mu-referencing under a mixture below.
• Class-shared typical values (V = TVV * exp(ETA_V), the same theta in every class) do take the closed-form mu-referencing shift \(\log\theta \mathrel{+}= \gamma\,\operatorname{mean}_i(\eta_i)\), exactly as in a single-population SAEM fit (#996).
• Per-class Ω overrides (omega(k)) use per-class stochastic-approximation diagonal statistics.
• Mixing coefficients are updated from a small optimisation of the SA-averaged expected complete-data mixing log-likelihood \(\sum_i \sum_k \bar r_{ik}\,\ln p_{ik}(\theta_\text{mix})\), where \(\bar r_{ik}\) is the Robbins–Monro average of the sampled class indicators. This is exact for constant mixing (it recovers \(p_k = \operatorname{mean}_i \bar r_{ik}\)) and for covariate-dependent logit mixing (a weighted logistic fit) alike.

The reported OFV, per-subject MIXEST/PMIX, and standard errors are computed from the K-fold mixture marginal at the SAEM optimum, exactly as for FOCEI. Per-class σ overrides (sigma(k)) are the one piece not yet re-estimated under SAEM — they are held at their initial values (with a warning), and every class-k subject is scored under that held σ, so the free base σ is estimated only from the classes that actually share it. Being held, the overrides report a standard error of 0, as FIXed parameters do; route the model to FOCEI if they must be fit. HMC and the conditional-distribution pass are disabled for a mixture (they are class-unaware); the MH E-step is used.

Two further SAEM-specific rules follow from the separated M-step:

• A mixing theta marked FIX is honoured — it is excluded from the mixing M-step and stays at its declared value, so a known mixture fraction can be pinned.
• A theta that drives both the mixing expression and a structural typical value is rejected with an error rather than silently double-owned (SAEM estimates the mixing coefficients and the structural thetas in separate M-steps). Split it into two thetas, or fit with FOCE/FOCEI, whose joint marginal handles the shared parameter.

Mu-referencing under a mixture

A MIXNUM-switched typical value pairs one η with K class thetas, so the classical pooled mu-referencing move \(\log\theta \mathrel{+}= \gamma\,\operatorname{mean}_i(\eta_i)\) (one theta per η) does not apply to it directly. ferx detects the class-aware form at parse time — a chain of MIXNUM == k branches whose arms are all THETA_k * exp(ETA) (or exp(log(THETA_k) + ETA)) on the same ETA — and resolves one anchor theta per class. Anything else (mixed anchors, a different ETA per arm, an additive arm, a non-MIXNUM condition, an out-of-range class) falls back to the numerical M-step, and the parser says so in a warning rather than degrading silently. The warning is only raised for expressions that actually carry an ETA: a class-switched typical value written without IIV (V = if (MIXNUM == 1) TVV1 else TVV2) or an intermediate MIXNUM flag variable is not a missed mu-reference and is left alone.

Each anchor theta belongs to one η. If two η are written against the same typical value, only the first keeps the closed-form shift — applying both would move that θ twice per iteration — and the other η’s θ go to the numerical / weighted M-step.

What each estimator does with that anchor set differs, and the difference is empirical:

Estimator	Class-shared θ	Class-switched θ
SAEM	closed-form pooled shift	numerical M-step
IMP / IMPMAP (estimating)	closed-form shift	responsibility-weighted closed-form shift

SAEM deliberately keeps the switched θ numerical. It draws one hard class per subject each E-step, so the per-class mean \(\operatorname{mean}_{i:c_i=k}(\eta_i)\) is a classification-EM statistic: the class boundary is re-drawn from the \(\theta_k\) that the shift just moved, and the two feed back. Measured on the two-class anchor (tests/nonmem/mixture_iv.csv, seed 20250818) the class-aware variant lands at TVCL2 = 3.02 with OFV 305.0, against 2.75 / 302.2 for the numerical M-step — worse on every coordinate (NONMEM SAEM: 2.735; FOCEI optimum: 2.842). A Rao-Blackwellised soft-responsibility variant converges to the same wrong point, so the bias is in the hard-class sufficient statistic rather than the weighting. IMP/IMPMAP importance-samples η within every class and weights by the responsibilities, so it has no hard assignment and no such feedback.

Setting mu_referencing = false disables the class-aware shift under both estimators and restores the pre-#996 behaviour (every θ through the numerical / weighted M-step). Note that for a single-population IMP/IMPMAP fit the closed-form shift is EM-mandatory and stays on regardless of that option; only the mixture path is switchable.

Two further conditions turn an individual θ’s shift off, on either path, each with a warning naming the parameter:

• A θ declared with a negative lower bound is packed on the identity scale, so log θ += mean(η) would apply θ += mean(η) where the closed form means θ *= exp(mean(η)). Give the parameter a non-negative lower bound to keep the closed-form update.
• A class that receives no responsibility mass in an iteration has no per-class η mean to shift with. Its θ is handed back to the weighted M-step for that iteration rather than held fixed, and the fit warns at the end; a class that collapses this way usually means the mixture is over-specified.

When the class-aware shift is off for a class θ — by any of these routes, or by mu_referencing = false — that θ is estimated by the importance-weighted M-step alone, which is biased for weakly-identified parameters (#406). IMP/IMPMAP fits say so in a warning; cross-check such a fit against FOCEI.

Under IOV, the per-occasion κ is sampled inside each subject’s drawn class — under that class’s MIXNUM branch and its Ω/Σ — so a class-switched typical value is never absorbed into κ.

NONMEM cross-check (tests/nonmem/mixture_iv_saem.ctl)

The two-class constant-mixing model above, estimated with $EST METHOD=SAEM (NBURN=1000 NITER=1000 ISAMPLE=2) followed by a METHOD=IMP EONLY=1 objective-evaluation pass, Ω/Σ FIXed:

Quantity	NONMEM SAEM	ferx SAEM
OFV (IMP evaluation)	300.87	302.17
TVCL1 (class 1)	1.002	0.982
TVCL2 (class 2)	2.735	2.746
V	9.993	10.089
mixing fraction p(1)	0.4712	0.4670
MIXEST classification	—	29 / 30 agree (1 borderline)

Estimates agree to ≤ 3 %. The single MIXEST disagreement is a subject whose posterior sits almost exactly on the class boundary, which NONMEM SAEM and ferx SAEM split marginally differently. Covariate-dependent mixing under SAEM is exercised end-to-end in tests/mixture_nonmem.rs.

Bayesian estimation under a mixture

Mixtures can also be fit with full-MCMC Bayesian estimation (method = bayes, #985). The latent class is handled Rao-Blackwellised — marginalised out rather than sampled — so the sampler still sees a smooth continuous target: each subject’s likelihood contribution is the K-class marginal \(-\log \sum_k p_{ik}\,\exp(-\mathrm{nll}_{ik})\). Because every Gibbs block evaluates the per-subject likelihood through this one marginal, the whole Gibbs-within-MH sampler becomes mixture-aware in one place:

• the η block samples \(\eta_i\) against the class-marginal posterior (a chol(Ω)-preconditioned random walk; HMC is disabled for mixtures because its analytic gradient is single-class);
• the (θ, σ) block — which also carries the mixing thetas — samples them against the marginal, so the mixing coefficients (constant or covariate-dependent logit mixing) are sampled with no extra machinery;
• the Ω block draws the (class-shared) Ω from its conjugate full conditional;
• the κ block (inter-occasion variability, #985) samples the per-occasion \(\kappa_{ik}\) against the same marginal, so κ targets the class-marginal posterior rather than class 1’s.

Ω/Σ are shared across classes (the classes differ only through the MIXNUM-switched typical values); per-class Ω/Σ overrides (omega(k)/sigma(k)) are rejected with a clear error — fit those with FOCE/FOCEI. The reported OFV, per-subject MIXEST/PMIX, and EBEs are the K-fold marginal values at the posterior mean.

Priors, label switching, and the reported OFV

The population priors are weakly informative normals, \(\mathcal{N}(u_0, 10^2)\) on the unconstrained scale (log θ where θ is positive, log σ), centred on the initial estimates \(u_0\) from the model file. Any [parameters] bounds you declare are not enforced as a uniform prior inside the (θ, σ) Metropolis block — only the mu-referenced conjugate θ move clamps to them — so a mixing logit is kept in range by the prior, not by its declared bounds.

Because the classes are identified only through those priors, a chain can swap labels (class 1 ↔︎ class 2) and the pooled posterior mean would then average across two modes and estimate neither class. The signature is a high split-R̂ on the class-switched θ; ferx emits a warning naming label switching when R̂ exceeds its convergence threshold on a mixture. Separating the classes in the initial estimates is the practical fix.

The OFV reported for a Bayes mixture run is the Laplace K-fold marginal at the posterior mean (not the posterior-mean joint NLL ×2 that a non-mixture Bayes run reports), and — like every mixture marginal — it is the FOCE or FOCEI form depending on interaction. In a method chain (method = [saem, bayes]) that flag follows the chain’s last estimator, so compare mixture OFVs only across runs with the same interaction setting.

Under diffuse priors the Bayes posterior mean concentrates at the maximum- likelihood optimum, so it recovers the same estimates the FOCEI/SAEM anchors pin (observed on tests/nonmem/mixture_iv.csv: posterior mean TVCL1 ≈ 0.94, TVCL2 ≈ 2.78, V ≈ 10.0, p(1) ≈ 0.47, marginal OFV ≈ 301.8, vs the FOCEI optimum 0.98 / 2.84 / 9.98 / 0.472 / 301.58; the two class means sit slightly inside the MLE pair, the direction partial label mixing pulls them). That shared optimum — already NONMEM-anchored by the FOCEI run above — is the Bayes cross-check: a direct NONMEM METHOD=BAYES reference is not used because NONMEM’s BAYES sampler aborts in burn-in on this model (it insists on Gibbs-sampling Ω, which is FIXed here), a known NONMEM $MIX+BAYES fragility. The same mixture-marginal machinery is anchored directly against NONMEM SAEM above.

IMP objective evaluation under a mixture

Importance sampling evaluates the class-marginal likelihood for a mixture (method = imp with imp_eval_only, NONMEM METHOD=IMP EONLY=1, #985). For each subject and each class it runs the class-k MAP inner solve (under the MIXNUM guard), importance-samples the class-conditional marginal \(L_{ik}\), and combines them exactly as FOCEI does:

\[ -2\log L = -2\sum_i \log \sum_k p_{ik}\,L_{ik}. \]

This is the Monte-Carlo marginal likelihood — the number to compare against NONMEM’s #OBJV for IMP EONLY — and it is what a chained fit reports on FitResult.importance_sampling. The reported per-subject ESS is the worst ESS among the classes that subject actually loads on (classes with \(\mathrm{PMIX}_{ik}\approx 0\) are excluded, since a poor proposal there barely enters the marginal), so imp_low_ess_threshold flags mixture subjects the same way it does single-population ones. The usual workflow is to estimate with SAEM (or FOCE/FOCEI) and append an evaluation stage:

[fit_options]
  method = [saem, imp]
  imp_eval_only = true

On tests/nonmem/mixture_iv.csv this recovers the NONMEM IMP EONLY objective to well under a Monte-Carlo standard error: ferx −2 log L = 300.82 ± 0.06 vs NONMEM 300.87 (mixture_iv_saem.ctl’s second $EST), with ess_min/K = 0.82 and ess_med/K = 0.91.

Estimating IMP / IMPMAP

Importance sampling also estimates a mixture (method = imp / impmap without imp_eval_only) via a class-partitioned MCEM. Each iteration importance- samples every subject’s η within each class (under the MIXNUM guard), forms the deterministic class responsibilities \(\mathrm{PMIX}_{ik} \propto p_{ik}L_{ik}\), then runs responsibility-weighted M-steps: the mixing coefficients from the responsibilities, the class-shared Ω from \(\sum_i\sum_k \mathrm{PMIX}_{ik}\,\widehat{\eta\eta^{\top}}_{ik}/N\), and the class-switched typical values from a responsibility- and importance-weighted observation M-step (each class’s samples scored under its own MIXNUM). Ω/σ are class-shared (per-class overrides are rejected); IOV / FREM / SDE are not supported (as for the single-population MCEM).

The class-switched typical values are mu-referenced per class (#996). Because the MCEM importance-samples η within every class, a per-class posterior mean \(\bar\eta_{ic}\) exists directly, and the EM move is the responsibility-weighted

\[\log\theta_k \mathrel{+}= \frac{\sum_i r_{ic}\,\bar\eta_{ic}}{\sum_i r_{ic}}, \qquad r_{ic} = \mathrm{PMIX}_{ic}\]

summed over the classes \(\theta_k\) serves. Those θ are pinned out of the weighted M-step, which then fits only σ and the remaining θ. A class-shared typical value repeats the same theta in every class slot, so it collapses to the pooled single-population shift. See Mu-referencing under a mixture for the detection rules and for why SAEM does not do the same.

This is what closes the old gap: on the two-class anchor, estimating IMPMAP now recovers TVCL1 = 0.952, TVCL2 = 2.821, V = 10.008, p(1) = 0.472, against 0.949 / 2.819 / 10.010 / 0.473 from NONMEM IMPMAP with MU_1 assigned inside the MIXNUM branch (tests/nonmem/mixture_iv_impmap_mu.ctl, NM 7.6.0) — better than 0.5 % on every coordinate. With mu_referencing = false the same run falls back to the weighted M-step alone and returns TVCL1 = 0.998, TVCL2 = 2.597, which is the pre-#996 behaviour and shows how much of the accuracy the shift buys. Note that the EM estimators settle slightly below the FOCEI optimum (0.98 / 2.84) on TVCL1/TVCL2; NONMEM’s mu-referenced IMPMAP does the same.

A θ whose paired ω is negligible (ω < 1e-3) carries no mean-shift information, so it is routed to the weighted M-step with a warning, as in the single-population MCEM; so is a θ declared with a negative lower bound (identity-packed, so the additive log-scale shift is not its EM optimum). The fit also names any estimated θ that carries no ETA at all — the mixing thetas are excluded, since they are owned by the mixing M-step.

The per-subject proposal machinery is the same as the single-population MCEM: each class’s draws get their own ISCALE pilot search (the guard against importance-weight collapse), and impmap_mceta runs its multi-start MAP within every class. Two IMPMAP options are not wired for mixtures and warn when set: impmap_trace (no iteration trace is produced) and impmap_auto / imp_auto (the sample count stays fixed).

Not yet supported

The following error clearly rather than silently mis-fitting:

• Per-class σ overrides (sigma(k)) under SAEM (held at their initial values, with a warning), and per-class Ω/σ overrides under Bayes and estimating IMP/IMPMAP (rejected) — use FOCE/FOCEI to fit those.
• IOV / FREM / SDE under the mixture IMP paths (objective evaluation and estimation) — use FOCE/FOCEI.
• A theta shared between the [mixture] mixing expression and a structural typical value, under SAEM and under estimating IMP/IMPMAP. Both run the mixing M-step separately from the θ/σ M-step, so a shared theta would be double-owned; FOCE/FOCEI’s joint marginal handles it.