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. All rows — dose records and observation records — that share the same occasion index belong to one occasion. Occasions can be supplied two ways:

• From the data — an integer column (e.g. OCC) named via iov_column in [fit_options].
• From the model — a rule (iov_occasion = dose or iov_occasion = time(...)) that derives the occasion of each row from its timeline, so you do not have to precompute an occasion column. See Declaring occasions in the model.

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.

Structural-model scope

IOV composes with every structural model, including the analytic absorption closed forms — pk one_cpt_transit / two_cpt_transit / one_cpt_ig / two_cpt_ig. Because a per-dose closed-form superposition cannot carry a dose’s drug across an occasion boundary where the disposition changes, a subject with IOV on such a model is transparently rerouted to the model’s exact transit() / igd() ODE twin, which integrates the cross-occasion carryover exactly (the same reroute used for time-varying covariates). This is automatic — no model change is needed — and covered on the Absorption page. A closed-form absorption model that has no ODE twin (one carrying a user [odes] / [scaling] / [initial_conditions] block, a role that shadows a reserved F/lagtime slot, or one whose twin could not be built — see When a model loses its ODE twin) still rejects IOV with an actionable message.

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.

Sample-size-weighted IOV (weight = ...)

Maturity: experimental (#1031).

The arm-level random effect of a longitudinal model-based meta-analysis — the between-treatment-arm variability (BTAV) — is not drawn from one variance. It scales with the number of subjects behind the arm:

\[ \kappa_{ik} \sim \mathcal{N}\!\left(0, \frac{\gamma^2}{N_{ik}}\right) \]

A 400-subject arm’s mean wanders a quarter as far as a 25-subject arm’s. Declare that where the rest of the variance structure is declared:

[parameters]
  kappa KAPPA_EMAX ~ 2.0 (sd) weight = NARM

[individual_parameters]
  LEMAX = LEMAX0 + log(OR_ABATA) * ABATA + ETA_EMAX + KAPPA_EMAX

[covariates]
  NARM continuous     # number of subjects in this arm

weight = W means κ ~ N(0, Ω_IOV / W). The engine applies the scaling by rewriting every reference to the kappa as KAPPA / sqrt(W) — precisely the hand-written form the modifier replaces — so the objective, the analytic sensitivities, the EBEs and every estimator see a model they already support.

Why declare it rather than write it

KAPPA_EMAX / sqrt(NARM) in the middle of a structural equation works, but it puts a variance-structure decision inside a structural expression:

• / NARM instead of / sqrt(NARM) produces a plausible wrong answer rather than an error;
• nothing in the model summary marks KAPPA_EMAX as weighted — it reports as an ordinary IOV term;
• the term can only be scaled where it is used, so the same divisor has to be repeated in every expression the kappa enters.

What is reported

The estimate stays the unweighted \(\gamma^2\) — the quantity a published MBMA reports — and the fit prints the effective SD at the median arm alongside it:

--- KAPPA (IOV) Estimates ---
  KAPPA_EMAX           = 4.145382  SE = 0.412000
                         weight = NARM  →  SD = 0.1440 at NARM = 200.0000 (κ ~ N(0, KAPPA_EMAX/NARM))

A raw \(\gamma\) of 2.0 on a logit scale reads as alarming until it is divided; \(\gamma/\sqrt{200} \approx 0.14\) is the between-arm SD a reader is looking for. FitResult carries both as kappa_weights (the source text) and kappa_weight_typical (the median arm size).

Rules and diagnostics

• The weight may reference covariates, FIXed thetas and TIME — not random effects, and not estimated thetas. A weight is the known precision of an occasion: a random effect there would make the variance of κ a function of κ, and an estimated theta would move the weight during the fit, so the positivity check below (evaluated at the initial θ) could not certify it.
• A weighted kappa may be referenced only in [individual_parameters]. Elsewhere it would read as the unweighted κ, so the parser rejects it — and it must be referenced somewhere in that block, since a weight that scales nothing would still be reported as if it had been applied.
• block_kappa cannot carry a weight: the divisor would have to apply to the block’s covariances too for the matrix to stay a valid one. Declare a weighted kappa as a diagonal kappa.
• A weight that is zero, negative or non-finite at any record is a fit-time error (E_KAPPA_WEIGHT_NONPOSITIVE) — the usual cause is a blank arm-size cell. simulate() enforces the same check. A weight that changes within an occasion warns (W_KAPPA_WEIGHT_VARIES_WITHIN_OCCASION): κ is drawn once per occasion, so Ω_IOV/W is ambiguous when W moves underneath it.

Validation against NONMEM

NONMEM has no weight = on an $OMEGA, so its equivalent model writes the scaling into $PK by hand — which is exactly what the declaration desugars into. Both engines were run on the warfarin IOV dataset with an arm-size column (NARM = 400 on occasion 1, 25 on occasion 2), 1-cpt oral, FOCE, same starting values (tests/nonmem/kappa_weight_iov.{ctl,lst,ext}):

Quantity	NONMEM 7.6.0	ferx
OFV (without constant)	207.455	205.683
TVCL	0.31727	0.32464
TVV	8.4018	8.4445
TVKA	2.6180	2.4467
ω²(CL)	0.59721	0.61281
ω²(KA)	1.0571	0.89725
γ² (Ω_IOV)	2.0013	1.8402
σ² (proportional)	0.037588	0.040884

ferx reaches a lower objective from the same start on a surface where κ is weakly identified (it enters divided by √400 and √25). The residual difference is the pre-existing cross-engine FOCE/IOV gap, not the weighting — evaluated at identical fixed parameters it is smaller with the weight than without:

Parameters	Model	NONMEM	ferx	Δ
inits, γ² = 0.01	weighted κ/√NARM	466.779	466.720	0.06
inits, γ² = 0.01	plain κ	332.265	331.620	0.65
NM MLE, γ² = 2.00	weighted κ/√NARM	207.455	205.943	1.51
NM MLE, γ² = 2.00	plain κ	238.889	238.524	0.37

The declaration itself is exact: weight = NARM and the hand-written KAPPA_CL / sqrt(NARM) score bit-identically in ferx under FOCE, FOCEI, Laplace and SAEM (tests/weighted_kappa_iov.rs), and the analytic ∂f/∂κ matches finite differences of the production predictor to tolerance.

This is the random-effect sibling of the residual weight = ... modifier, and the two share syntax and validation.

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.

Declaring occasions in the model

Instead of carrying an occasion column in the data, you can have ferx derive the occasion partition from each subject’s timeline with the iov_occasion key in [fit_options]. The kappa declarations are unchanged — only where the occasion labels come from differs.

iov_occasion value	Meaning
column (default)	Occasions come from the iov_column dataset column.
dose	Each administration time begins a new occasion. An observation takes the occasion of the most recent dose at or before its time (dose-before-observation at equal TIME). Doses that share a time (e.g. a simultaneous depot + central bolus) share one occasion. ADDL-expanded doses each start their own occasion.
time(t₁, t₂, …)	Time-window breakpoints. The strictly-increasing edges define bins: time(24, 48) → occasion 0 for TIME < 24, occasion 1 for 24 ≤ TIME < 48, occasion 2 for TIME ≥ 48. The first occasion runs from -∞ to the first breakpoint and the last runs from the last breakpoint to +∞, so you specify only the interior boundaries.

For example, a two-period crossover with a washout at 24 h needs no OCC column:

[fit_options]
  method       = foce
  iov_occasion = time(24)

This single breakpoint defines two occasions: TIME < 24 (occasion 0) and TIME ≥ 24 (occasion 1). You give only the interior boundaries — the breakpoint list is not a list of occasion start times. In particular, do not add a leading 0: occasion 0 already starts at -∞, so time(0, 24, 48) is a mistake (it would carve out an empty occasion for TIME < 0 and shift every real occasion up by one). ferx rejects a 0 breakpoint with an error pointing you at the correct time(24, 48) form. Likewise there is no trailing edge for the last occasion — it runs to +∞.

Precedence. If both iov_occasion (set to dose/time(...)) and iov_column are given, the model-side rule wins and a warning notes that the column was ignored.

Clock. Both rules bucket observations and doses on ferx’s internal event clock — the same monotonic timeline the solver integrates on. For reset-stacked subjects (EVID=3/4, where the data-file TIME restarts each period) that clock keeps advancing across the reset, so time(...) breakpoints separate the periods correctly and each dose stays in the same occasion as its co-timed observations.

Identifiability. A rule that yields only one occasion for every subject (e.g. time(...) breakpoints outside the sampling window, or dose on single-dose subjects) leaves kappa confounded with between-subject variability, so fit() errors — set breakpoints inside the sampling window, use dose on a multiple-dose design, or supply an iov_column. On the other end, dose on a long ADDL/steady-state train produces one occasion per expanded dose; this can balloon the per-subject inner problem, so ferx warns when the derived occasion count is large — prefer time(...) windows there. If the model declares no kappa, iov_occasion has no effect and a warning says so.

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.

Simulation

simulate() is occasion-aware for a kappa model: it draws one independent κ ~ N(0, Ω_IOV) per occasion group, so a VPC from an IOV model carries the inter-occasion spread the fit estimated (drawing a single κ for the whole subject would silently under-disperse it).

That draw needs Ω_IOV in the parameters you hand to simulate(). The model’s own parameters (model.default_params) always carry it, but a caller that rebuilds parameters from a fit must thread the fitted IOV covariance through as well — estimation::uncertainty_samples::fitted_params_from_result() does this and is the supported way to go from a FitResult to a simulate-ready ModelParameters. Parameters with a missing omega_iov are rejected: simulate_with_options() / simulate_with_options_diag() return an Err naming the omission, and the Vec-returning simulate() / simulate_with_seed() panic with the same message rather than silently emitting occasion-free draws (issue #1019).

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.
• Overlapping-occasion absorption tails use the current occasion’s kappa (built-in input-rate forcings only — transit(), igd(), weibull(), first_order()). Each integration segment evaluates its absorption shape parameters (e.g. MTT, MAT, KA) from that segment’s occasion snapshot, so if a dose from an earlier occasion is still absorbing when a new occasion begins, its remaining input rate is computed with the new occasion’s parameters. This is exact for pure IIV and whenever the inter-occasion interval exceeds the absorption window (the usual washout case); only genuinely overlapping-occasion absorption tails are approximated. Bioavailability F, lagtime, and zero_order(dur) windows are read per-dose from each dose’s own occasion snapshot and stay occasion-correct even under overlap.
• 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 from data	OCC in $INPUT; referenced in $PK	iov_column = OCC in [fit_options]
Occasion from model	IF(TIME.LT.24) OCC=1 logic in $PK	iov_occasion = dose or time(...) in [fit_options]
Kappa in PK	ETA(n) with SAME block	KAPPA_xxx in [individual_parameters]