Inter-Occasion Variability (IOV)

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

Inter-Occasion Variability (IOV) models the fact that a subject’s PK parameters may shift between occasions (treatment periods, study visits, dosing intervals). Unlike between-subject variability (BSV), which is fixed for a subject across the whole study, IOV is a random effect that is re-drawn for each occasion.

Concepts

Occasions

An occasion is a distinct time window during which the same kappa value applies. Occasions are identified by an integer column in the dataset (e.g. OCC). All rows — dose records and observation records — that share the same occasion index belong to one occasion.

Kappa parameters

Kappa (κ) is the IOV random effect, analogous to eta (η) for BSV. Each kappa is drawn independently per occasion:

\[ \kappa_{ik} \sim \mathcal{N}(0, \Omega_\text{IOV}) \]

where i indexes subjects and k indexes occasions. The IOV omega matrix Ω_IOV is estimated alongside the BSV omega.

Individual parameters

Kappas enter the individual-parameter expressions in the same way as etas:

CL = TVCL * exp(ETA_CL + KAPPA_CL)

At each occasion k, the effective individual CL is TVCL * exp(ETA_CL_i + KAPPA_CL_ik). The BSV eta captures stable between-subject differences; the kappa captures occasion-to-occasion fluctuation within a subject.

Option A: Diagonal IOV (independent kappas)

Each kappa is declared independently, giving a diagonal Ω_IOV:

[parameters]
  ...
  kappa KAPPA_CL ~ 0.05
  kappa KAPPA_V  ~ 0.03

The two kappas are uncorrelated across occasions. This is the most common formulation and corresponds to NONMEM’s IOV Option A.

Use FIX to hold a kappa variance constant during estimation:

kappa KAPPA_CL ~ 0.05 FIX

Option B: Block IOV (correlated kappas)

When occasion effects on different parameters are expected to covary, use block_kappa:

[parameters]
  ...
  block_kappa (KAPPA_CL, KAPPA_V) = [0.05, 0.01, 0.03]

Values are the lower triangle of Ω_IOV: Var(KAPPA_CL), Cov(KAPPA_CL, KAPPA_V), Var(KAPPA_V). This mirrors the block_omega syntax for BSV.

Mixed Option A + B

kappa and block_kappa declarations can be combined freely — uncorrelated kappas can coexist with a correlated block:

block_kappa (KAPPA_CL, KAPPA_V) = [0.05, 0.01, 0.03]
kappa KAPPA_KA ~ 0.10

A name may not appear in both a kappa and a block_kappa line.

Model File Setup

[parameters]
  theta TVCL(0.2, 0.001, 10.0)
  theta TVV(10.0, 0.1, 500.0)
  theta TVKA(1.5, 0.01, 50.0)

  omega ETA_CL ~ 0.09          # BSV
  omega ETA_V  ~ 0.04
  omega ETA_KA ~ 0.30

  kappa KAPPA_CL ~ 0.05        # IOV (Option A)

  sigma PROP_ERR ~ 0.02

[individual_parameters]
  CL = TVCL * exp(ETA_CL + KAPPA_CL)
  V  = TVV  * exp(ETA_V)
  KA = TVKA * exp(ETA_KA)

[structural_model]
  pk one_cpt_oral(cl=CL, v=V, ka=KA)

[error_model]
  DV ~ proportional(PROP_ERR)

[fit_options]
  method     = foce
  iov_column = OCC

See examples/warfarin_iov.ferx for a complete runnable example.

Algorithm Details

Inner loop

When a model has kappa declarations and the subject has occasion labels, the inner optimizer runs find_ebe_iov instead of the standard find_ebe. It jointly optimizes over:

\[ p = [\underbrace{\eta_1, \ldots, \eta_{n_\eta}}_{\text{BSV}},\ \underbrace{\kappa_{1,1}, \ldots, \kappa_{1,n_\kappa}}_{\text{occasion 1}},\ \ldots,\ \underbrace{\kappa_{K,1}, \ldots, \kappa_{K,n_\kappa}}_{\text{occasion K}}] \]

The joint negative log-posterior is:

\[ -\log p(p \mid y_i) = \frac{1}{2}\left[ \eta^T \Omega^{-1} \eta + \log|\Omega| + \sum_{k=1}^{K} \kappa_k^T \Omega_\text{IOV}^{-1} \kappa_k + K \log|\Omega_\text{IOV}| + \sum_{j} \left(\frac{(y_{ij} - f_{ij}^{(k_j)})^2}{V_{ij}} + \log V_{ij}\right) \right] \]

where $k_j$ is the occasion of observation j and $f_{ij}^{(k)} = f(\theta, \eta_i, \kappa_{ik})$.

Optimization uses BFGS. Within the analytical 1-/2-/3-compartment scope the inner-loop H-matrix (the Jacobian ∂f/∂[η, κ] over the stacked random effects) is taken from the exact analytic second-order-dual sensitivity provider rather than finite differences. User-[odes] IOV models within the event-driven analytic scope (a compiled RHS, bolus/finite-infusion dosing, resets and EVID=2 breakpoints, dose-only occasions, and a per-subject stacked axis count up to 96) use the same exact analytic stacked-η inner gradient. Models outside both scopes — steady-state IOV doses, IIV-on-residual-error (iiv_on_ruv), FREM, an estimated lagtime, or a per-subject axis count above the cap — keep the finite-difference Jacobian.

Outer loop

The IOV omega is packed after the sigma block in the optimizer vector. For Option A (diagonal Ω_IOV) only the log-diagonal entries are appended; for Option B the full Cholesky lower triangle is appended (log-diagonals plus off-diagonals as-is), mirroring the BSV omega packing:

\[ x = [\log\theta,\ \text{chol}(\Omega_\text{BSV}),\ \log\sigma,\ \text{chol}(\Omega_\text{IOV})] \]

The per-subject FOCE objective is a proper augmented marginal: the per-occasion kappas are integrated out through the linearised marginal exactly like the BSV etas. The random-effect vector is $b = [\eta,\ \kappa_1, \ldots, \kappa_K]$ with block-diagonal prior covariance

\[ \Sigma_b = \text{blockdiag}(\Omega_\text{BSV},\ \Omega_\text{IOV},\ \ldots,\ \Omega_\text{IOV}) \]

(K copies of Ω_IOV), and the H-matrix is augmented with kappa columns, $H_\text{full} = [\,\partial f/\partial\eta \mid \partial f/\partial\kappa_1 \mid \cdots \mid \partial f/\partial\kappa_K\,]$. A kappa column can affect later observations through cross-occasion carryover, including when the kappa belongs to a dose-only occasion. The objective is then the ordinary Sheiner–Beal form

\[ \text{NLL}_i = \tfrac{1}{2}\left[(y - f_0)^T \tilde{R}^{-1} (y - f_0) + \log|\tilde{R}|\right],\qquad \tilde{R} = H_\text{full}\,\Sigma_b\,H_\text{full}^T + R \]

with no separate kappa prior added (it is already folded into $\tilde{R}$). This reduces exactly to the BSV-only FOCE objective when there are no kappas.

Why this matters (issue #101). An earlier implementation linearised only the BSV etas and added an explicit kappa MAP penalty $\tfrac{1}{2}\sum_k(\kappa_k^T\Omega_\text{IOV}^{-1}\kappa_k + \log|\Omega_\text{IOV}|)$. That penalty omits the kappa Laplace determinant, so Ω_IOV was unidentified and the FOCE OFV disagreed with SAEM and NONMEM. The augmented marginal above fixes this; all estimation paths now share one consistent objective.

Outer-loop gradient (EBE re-convergence for IOV)

The IOV variance components — especially Ω_IOV — are weakly identified, and their gradient is dominated by the EBE response. Raising Ω_IOV un-shrinks the per-occasion kappas and improves the data fit, an effect a fixed-EBE gradient cannot see, so a gradient that holds the EBEs fixed leaves Ω_IOV pinned at its initial value (while derivative-free methods, which re-solve the EBEs at each trial point, move it freely).

For IOV models within the analytical 1-/2-/3-compartment scope, under a gradient-based outer optimizer (lbfgs, bfgs, slsqp), the outer gradient is now the exact analytic marginal gradient (Almquist et al. 2015) assembled over the stacked random-effect vector [η_bsv, κ₁, …, κ_K] with the block-diagonal prior Ω_bsv ⊕ K·Ω_iov. It carries the EBE response (Eq. 46) in closed form on every θ/Ω/σ/Ω_iov block — including the κ-variance — so it sees the un-shrinkage effect directly and converges Ω_IOV without re-solving the inner loop. Both FOCEI (Laplace marginal) and FOCE (Sheiner–Beal linearized marginal) have exact closed-form gradients here, and the cross-occasion dose carryover is differentiated exactly through the second-order-dual event-driven walk. Validated against reconverged finite differences (~1e-6 on every packed parameter) and against NONMEM on the warfarin IOV model.

User-[odes] IOV models get the same exact analytic outer gradient when they are in the event-driven analytic scope (a compiled RHS, bolus or finite-infusion dosing with or without time-varying covariates, EVID=2 breakpoints, resets, dose-only occasions, and a per-subject stacked axis count up to 96), assembled by the same second-order-dual walk fed per-occasion parameters. The outer gradient is assembled per subject: subjects in analytic scope use the exact gradient, and the few outside it (steady-state IOV doses, an estimated lagtime, iiv_on_ruv, FREM, an above-cap axis count, or a non-finite analytic component) fall back to central finite differences that re-converge that subject’s inner loop at each perturbed point (subject_reconverged_fd_gradient_iov) — so one out-of-scope subject no longer forces the whole population onto finite differences. The default BOBYQA optimizer is derivative-free and uses neither path. Setting reconverge_gradient_interval = 1 forces the reconverged-FD gradient even within the analytical scope, as a documented escape hatch.

In addition, on the SLSQP path IOV models auto-enable per-coordinate scaling. SLSQP applies a uniform gradient cap that rescales the whole gradient by its largest (theta) component; with the disparate-magnitude IOV parameters (block-diagonal omega + kappa) that cap otherwise starves the omega/Ω_IOV step and leaves the variances pinned at their initial values. Presenting O(1) scaled coordinates removes the starvation, so a pure FOCEI/SLSQP fit moves the variance components off their initial values and descends toward the marginal minimum from a cold start. It does not always reach it: the final descent is driven by a noisy re-converged FD gradient on a weakly-identified variance surface, and SLSQP can terminate a few OFV units high in a platform-dependent way (see Limitations). The default BOBYQA optimizer and the methods = [saem, focei] chain reach the true minimum platform-independently. (This scaling is scoped to IOV + SLSQP: it is unnecessary for BFGS, which optimizes in unscaled space, and counter-productive for MMA — and the global scale_params default stays off to preserve the issue-#99 non-IOV behaviour.)

Covariance step

When covariance = true, standard errors for the IOV omega diagonal (or its Cholesky elements) are computed via the finite-difference Hessian at convergence, using the same procedure as for BSV omega.

Output

The following IOV-related fields are populated on FitResult and printed by print_results(). The fit YAML gains an omega_iov: section listing the kappa variances and (for Option B) the off-diagonal covariances/correlations.

Field	Description
`omega_iov`	Estimated IOV covariance matrix (`Option<DMatrix<f64>>`)
`kappa_names`	Names of kappa parameters in declaration order
`kappa_fixed`	Per-kappa `FIX` flags
`se_kappa`	Standard errors for IOV omega elements (requires `covariance = true`)
`ebe_kappas`	Per-subject, per-occasion kappa EBEs. `ebe_kappas[i][k]` is the kappa vector for subject i, occasion k

When the dataset has occasion labels the sdtab CSV gains an OCC column (one row per observation, the subject’s occasion index for that row).

Note: shrinkage_kappa exists as a placeholder field on FitResult but is not yet computed — it is always returned as an empty Vec. Use ebe_kappas directly if you need a manual shrinkage diagnostic.

Limitations

Cross-occasion dose carryover. Predictions are computed by pk::predict_iov, which builds per-event PK parameters carrying each event’s occasion kappa and propagates the compartment amounts continuously across occasion boundaries (via the event-driven solver). A dose from an earlier occasion is therefore eliminated through a later occasion with the later occasion’s clearance — matching NONMEM’s integration model (this replaced the earlier “Option A” superposition, which used a single clearance for the whole dose history and biased no-washout designs; issue #104). On the warfarin IOV reference, ferx’s population prediction now matches NONMEM’s to 5 significant figures on every row. A ~17 OFV-unit gap versus NONMEM remains at matched parameters, but it is not a prediction difference — it is a cross-engine FOCEI marginal/EBE difference for this weakly-identified multi-random-effect model (and NONMEM’s reference itself did not converge cleanly here). See tests/warfarin_iov_nonmem.rs.
Local-optimizer convergence from a far-off cold start. Pure FOCEI (SLSQP via the auto-scaling above, or BFGS) reaches the minimum from a sensible start, but the IOV variance surface is ill-conditioned and weakly identified, so a far-off start (e.g. a residual-error init off by more than ~2×) can still stall. Set a realistic sigma init (a too-small one drives early omega inflation). The most robust workflow remains SAEM, or a methods = [saem, focei] chain (SAEM finds the basin, FOCEI polishes to a strictly better marginal optimum). See tests/iov_convergence.rs.
Platform-dependent pure-SLSQP cold-start termination. Even from the model’s cold default start, a pure FOCEI/SLSQP fit on warfarin_iov reaches a different point on different architectures: macOS arm64 lands at ≈307.8 (the minimum), while Linux x86_64 deterministically stalls at ≈314.7, ~7 OFV units high (issue #160). The descent direction near the minimum comes from a re-converged central-FD gradient whose floating-point summation order differs across architectures, and SLSQP’s stagnation guard latches at whichever flat-looking point it reaches first. This affects only an explicit optimizer = slsqp; the default BOBYQA (derivative-free) and the SAEM→FOCEI chain reach ≈307.8 on both platforms. Prefer them for IOV fits — reserve pure SLSQP for a warm start already in the basin.
ODE Form C output ([scaling] y = <expr>) cannot reference KAPPA_* directly under IOV — the readout is evaluated once per observation with a single eta, so a direct kappa reference would see kappa = 0. The parser rejects it with a clear error (issue #107); put the occasion dependence in the structural parameters (e.g. CL) and reference those in the readout instead. The PK dynamics themselves are always occasion-correct.
Inner-loop and outer gradients come from the exact analytic sensitivity provider within scope for IOV models — both the analytical 1-/2-/3-compartment models and user-[odes] IOV models with a compiled RHS, event-driven dosing, dose-only occasions, and a per-subject stacked axis count up to 96 (see Inner loop and Outer-loop gradient). Steady-state IOV doses, iiv_on_ruv, FREM, an estimated lagtime, and above-cap axis counts still fall back to finite differences (per subject for the outer gradient).
The occasion column must contain non-negative integers. Non-integer or negative values are silently treated as missing (OCC = 0) and a single summary warning is emitted.
Per-kappa shrinkage is not yet computed (see Output note above).
Kappa EBEs are not yet emitted as columns in the sdtab CSV; access them via FitResult.ebe_kappas.

Comparison to NONMEM

Concept	NONMEM	ferx-core
Diagonal IOV	`OMEGA BLOCK(n) SAME` with one kappa per block	`kappa NAME ~ variance` (Option A)
Correlated IOV	`OMEGA BLOCK(n)` shared across occasions	`block_kappa (N1, N2) = [...]` (Option B)
Occasion column	`OCC` in `$INPUT`; referenced in `$PK`	`iov_column = OCC` in `[fit_options]`
Kappa in PK	`ETA(n)` with `SAME` block	`KAPPA_xxx` in `[individual_parameters]`