Advanced usage

Independent settings, each needing a paragraph more than a reference page gives it. The sections have nothing else in common and stand alone.

Full covariance vs variance NLL

Each study’s V selects a branch of the NLL, auto-detected:

V structure	Auto-detected `method`	NLL cost
Non-diagonal matrix	`"cov"`	O(n_t³) Cholesky solve
Diagonal matrix	`"var"`	O(n_t) element-wise
Plain vector (variances only)	`"var"`	O(n_t) element-wise

Override by setting method explicitly in the study list:

# Force full covariance NLL even for a diagonal V
study_cov <- c(study, list(method = "cov"))

# Diagonal approximation when you only have marginal variances
study_var <- list(
  E      = E,
  V      = diag(diag(V)),   # drop off-diagonal entries
  n      = n,
  times  = times,
  ev     = rxode2::et(amt = 100),
  method = "var"
)

Use "cov" when you have a full covariance matrix and expect the temporal correlation to inform the estimates; "var" when only marginal variances exist — a published table of means and SDs — or when runtime matters more.

Gradient modes

Every control takes a grad argument. admControl() names its analytical mode "sens"; adfoControl(), adghControl() and adirmcControl() call theirs "analytical". All four also take "fd" and "none".

"fd" is a central difference. Forward differencing was removed in 0.4.1: it measured 10² to 10⁴ times less accurate than the analytic gradient, and the one solve per parameter it saved did not pay for a gradient the optimizer struggles to descend. Nor is the step a constant — it is measured per parameter by the Shi (2021) procedure, with grad_h the fallback when it cannot be.

`grad =`	What it does	Notes
`"sens"` (`admc`) / `"analytical"` (the other three)	Sensitivity equations, contracted in closed form	The default everywhere. Needs an ODE or `linCmt()` model; `adfo` gets its structural thetas from the order-2 sensitivity model, and falls back to FD only if that cannot be built
`"fd"`	Central finite difference of the full NLL	`2 n_p` evaluations per step. Also the automatic fallback where the sensitivity model cannot be built, with a warning
`"none"`	BOBYQA, derivative-free	No gradient; `adfo`’s default up to 0.4.0

The GH objective is noise-free – deterministic quadrature, no MC draws – so its analytical gradient is exact and the Hessian well-conditioned. Keep the default.

For adfo, every mode still uses the sensitivity model inside each NLL evaluation to get J in one rxSolve; grad only controls how the optimizer gradient is formed.

# Analytical gradient (default; fastest for ODE models)
fit_sens <- nlmixr2(pk_model, admData(), est = "admc",
  control = admControl(studies = list(s = study), grad = "sens",
                       n_sim = 5000L, seed = 1L))

# Derivative-free — ignores maxeval; use nloptr stopping criteria instead
fit_bobyqa <- nlmixr2(pk_model, admData(), est = "admc",
  control = admControl(studies = list(s = study), grad = "none",
                       n_sim = 5000L, seed = 1L))

Mu-referencing and sensitivity equations

Mu-referencing is the nlmixr2 convention of writing each structural parameter as a fixed effect plus a random effect under a known back-transformation:

cl <- exp(tcl + eta.cl)   # mu-referenced: tcl paired with eta.cl

admixr2 uses this pairing to classify parameters into three cases, each handled differently by grad = "sens":

Paired (mu-referenced) parameters — rxode2 augments the ODE system with sensitivity equations d(pred)/d(eta_i). Since tcl and eta.cl enter additively on the log scale, the NLL gradient with respect to tcl equals the one with respect to eta.cl, and a single solve recovers both.

Non-mu-referenced parameters — written separately (cl <- exp(tcl) * exp(eta.cl)), the eta-keyed sensitivities do not carry d(pred)/d(tcl). admixr2 therefore emits its own first-order sensitivity model over an explicit direction set — one direction per random effect, one per unpaired theta — compiled with eventSens = "jump" so dosing-modifier (f/lag/rate/dur) sensitivities come too. The unpaired thetas then get an analytical gradient from that same solve, no finite differences, with a message saying so. This is the direction-set scheme nlmixr2est’s fast-focei uses, at first order, cross-validated against its inner model to ~1e-13. A targeted CRN finite-difference solve is the fallback when the augmented model cannot be built.

Parameters without a random effect — one with no corresponding eta (v2 <- exp(tv2), IIV dropped from V2) is the same unpaired case: its own direction in that sensitivity model, and an exact analytical gradient. (adfo needs one thing more: V_pred = J Ω Jᵀ + Σ depends on the second derivative dJ/dθ, which comes from the order-2 sensitivity model. Finite differences are the fallback when that cannot be built, not the rule.)

IRMC kappa correction

IRMC shifts its importance weights by moving the proposal mean when a structural parameter changes. For a mu-referenced parameter that shift is a closed-form function of the back-transformed ratio.

A theta with no paired random effect (ka <- exp(tka), no eta.ka) has no such weight path: changing tka in the inner loop moves f(theta, 0) and nothing else follows. The kappa correction adds kappa = f(theta_cand, 0) - f(theta_outer, 0) to the predicted mean, anchoring the inner NLL to the true prediction at the candidate. Kappa is zero when every theta is mu-referenced.

How f(theta_cand, 0) gets evaluated is the choice:

"exact" (default) re-solves it at each inner NLL evaluation — one extra rxSolve per inner step.
"linearized" precomputes J = df/d(theta) once per outer iteration from a small FD batch and takes kappa_fn(theta_cand) ≈ f0 + J %*% (theta_cand - theta0) — arithmetic only, no rxSolve at all.

The approximation holds while inner steps stay small against the outer box constraint, which is typical once phases converge, so prefer "linearized" wherever a solve is expensive.

adirmcControl(..., kappa_method = "exact")       # default
adirmcControl(..., kappa_method = "linearized")  # faster for complex ODE models

Parallel restarts

Multi-restart fitting guards against local optima. workers > 1 runs the restarts in parallel, over a pool admixr2 manages itself:

fit_par <- nlmixr2(
  pk_model, admData(), est = "admc",
  control = admControl(
    studies    = list(examplomycin = study),
    n_sim      = 5000L,
    n_restarts = 4L,
    workers    = 4L,      # background worker processes (mirai daemons)
    cores      = 8L,      # total rxSolve OpenMP threads across workers
    restart_sd = 0.3,     # SD of log-scale perturbation from the starting point
    seed       = 1L
  )
)

cores is split automatically — floor(total_cores / n_workers) each, remainder cyclically, so 13 cores over 4 workers gives 4, 3, 3, 3. Workers stop after the restart phase, freeing every core for the covariance Hessian. Call admStopWorkers() by hand if a fit is interrupted before that cleanup.

Quasi-random sampling

The sampling argument controls how eta samples are drawn:

`sampling =`	Method
`"sobol"`	Sobol sequence (default)
`"halton"`	Halton sequence
`"torus"`	Kronecker torus
`"lhs"`	Latin hypercube
`"rnorm"`	Plain normal

Sobol typically needs 2–5× fewer samples than plain normal draws for the same NLL variance; at small n_sim, "lhs" covers better.

Parameter uncertainty

covMethod = "r" takes a numerical Hessian after optimisation and reports standard errors in print(fit). The default is covMethod = "r,s" (below), the sandwich correction of that same Hessian; ask for "r" to get the uncorrected ones. SEs cover the structural, residual-error and omega parameters, omega on its natural scale, with the Hessian rotated off the Cholesky factor the optimiser works in. A larger cov_n_sim cuts MC noise in the Hessian:

fit_cov <- nlmixr2(
  pk_model, admData(), est = "admc",
  control = admControl(
    studies     = list(examplomycin = study),
    n_sim       = 5000L,
    cov_n_sim   = 10000L,   # more samples → lower Hessian noise
    cov_h_outer = 2.5e-3,   # outer FD step scale (default: eps^(1/5) ≈ 7.4e-4)
    covMethod   = "r",
    seed        = 1L
  )
)

Skip uncertainty altogether during early development, when only the point estimates matter:

admControl(..., covMethod = "none")

If the Hessian is non-positive-definite (SEs printed as NA), increase cov_h_outer or cov_n_sim.

`covMethod = "r,s"`: scoring the summary against its own sampling law

The aggregate objective is the exact log-likelihood of n iid draws from N(yt, Vt), which assumes each subject’s observation vector is multivariate normal. It is not: y_i = f(theta, b_i) + eps_i is nonlinear in the random effect, so the marginal is a mixture. The point estimates do not pay for this — the score has expectation zero at the true parameters under any weight, so the fit stays consistent. The reported uncertainty does, in two ways that differ in kind: Cov(V_ij, V_kl) is mis-sized by the excess kurtosis, and Cov(ybar, vech V) is assumed zero where a real correlation of 0.3–0.6 sits. Sample mean and sample covariance are exactly independent for a multivariate normal and for nothing else.

covMethod = "r,s" scores (ybar, vech V) against its own asymptotic law and reports the sandwich H^-1 J H^-1, where H is the same Hessian "r" inverts and J is the variance of the score built from that law (the model-derived weight above, not an empirical fourth moment):

fit_rs <- nlmixr2(
  pk_model, admData(), est = "adgh",
  control = adghControl(studies = list(examplomycin = study), covMethod = "r,s")
)
fit_rs$covMethod   # "r,s" if the correction was applied, "r" if it degraded

Point estimates and objective are identical to the "r" fit; only the SEs move. Under correct specification J = 2H, so "r,s" returns exactly what "r" would: the two agreeing is the expected outcome for a well-specified model, not a sign that nothing happened.

It covers every residual family whose conditional law is independent across timepoints — all of them except ar(): the conditionally-normal set (add, prop, pow, combined1, combined2), the closed-form distributional ones (lnorm, pois, binom, nbinomMu, beta, t(nu > 4)) and the transform-both-sides ones (boxCox, yeoJohnson, logitNorm, probitNorm). What it cannot build — ar(), t() with nu <= 4, ordinal(), a same-subject joint study, a singular ingredient — degrades to "r" with a warning, and fit$covMethod then reads "r": it records what the covariance is, not what was asked for.

Under "r,s" the denominator of a study’s reported V stops being cosmetic — see v_denom in the aggregate data vignette.

Model comparison: AIC and BIC

Standard information criteria work directly on admFit objects. AIC and BIC compare only across estimators evaluating the same likelihood. admc, adgh and adirmc all target the exact aggregate MVN likelihood — MC, quadrature and importance-weighted estimates of one integral — so they are mutually comparable. adfo’s is linearised, and must not be set against them.

Full model (IIV on all five parameters) against a reduced one (IIV on CL and V1):

ctl <- admControl(
  studies   = list(examplomycin = study),
  n_sim     = 5000L,
  cov_n_sim   = 10000L,
  maxeval   = 300L,
  seed      = 1L
)

fit_full    <- nlmixr2(pk_model,   admData(), est = "admc", control = ctl)
fit_reduced <- nlmixr2(pk_reduced, admData(), est = "admc", control = ctl)
#> 
#> 
#> 

AIC(fit_full, fit_reduced)
#>             df       AIC
#> fit_full    11 -3668.262
#> fit_reduced  8 -3523.183
BIC(fit_full, fit_reduced)
#>             df       BIC
#> fit_full    11 -3597.732
#> fit_reduced  8 -3471.888

The model with the lower AIC/BIC is preferred; a difference > 10 is generally considered strong evidence. BIC penalises added parameters more heavily than AIC, so the two can disagree on a marginal eta.

`anova()`: a likelihood-ratio test on nested fits

The two are nested – fit_reduced is fit_full without eta.v2, eta.q and eta.ka – so a chi-squared test compares their objectives directly:

anova(fit_full, fit_reduced)
#> Likelihood-ratio test
#>        Npar  OBJF   AIC   BIC   Test  dOFV Df         p
#> admc      8 -3539 -3523 -3472     NA    NA NA        NA
#> admc_1   11 -3690 -3668 -3598 1 vs 2 151.1  3 1.541e-32

dOFV is the objective difference and Df the number of parameters the larger model adds (three etas here); p is pchisq(dOFV, Df, lower.tail = FALSE), the ordinary likelihood-ratio test.

anova() refuses fits that are not on the same footing. Different estimators score different approximations (FO-linearised, quadrature, Monte Carlo), and changing n_nodes or n_sim moves the objective with the grid or sample size rather than with the model. All three are errors rather than a number returned quietly. See ?anova.admFit for the boundary-variance case — testing whether an eta’s variance is zero — and why a negative dOFV is reported, not clamped.