Package {admixr2}

Named list of study specifications (same format as admControl(): E, V, n, times, ev, optional method; or an observations list for multi-compartment fits – see admControl()).

Gradient mode. "analytical" (default) uses the closed-form FO gradient with LBFGS; "none" uses derivative-free BOBYQA; "fd" uses central finite differences of the full NLL. Forward differencing was removed in 0.4.1 – it was 10^2 to 10^4 times less accurate than a central difference at every site measured, and the one solve per parameter it saved did not pay for a gradient the optimizer struggles to descend.

The default was "none" up to 0.4.0, because the structural thetas were finite-differenced through the whole NLL and the resulting gradient was too noisy for a quasi-Newton step to pay off. They are now differentiated analytically from a second-order sensitivity model (relative error ~1e-7 against a central difference, where the finite-difference pass reached 1e-2), so LBFGS on the exact gradient is the better default. A model that cannot build that sensitivity model falls back to the finite-difference gradient automatically, and grad = "none" remains available.

nloptr algorithm, or NULL (default) to pick the default that matches grad: "NLOPT_LD_LBFGS" with a gradient, "NLOPT_LN_BOBYQA" when grad = "none". Any algorithm reported by nloptr::nloptr.print.options() is accepted. An explicit algorithm is reconciled with grad: when grad = "none" a gradient-based algorithm (⁠NLOPT_LD_*⁠ / ⁠NLOPT_GD_*⁠) falls back to "NLOPT_LN_BOBYQA"; when a gradient is requested a derivative-free algorithm (⁠NLOPT_LN_*⁠ / ⁠NLOPT_GN_*⁠) turns the gradient off. Both emit a message.

Maximum function evaluations (default 500).

Relative tolerance (default sqrt(.Machine$double.eps)).

Print-frequency for live progress (0 = silent).

Random seed (used for restarts).

OpenMP threads for rxSolve(). Defaults to rxode2::rxCores(). When workers > 1 it is a total budget, split across the workers.

Passed to rxSolve(): show the solver's text progress bar only once a single solve exceeds this many subjects. The default (.Machine$integer.max) keeps it off for clean script/vignette output; lower it (e.g. 1000L) to see progress during long fits.

Finite-difference step for unpaired struct theta gradient and FD Jacobian.

Box-constraint half-width when using gradients: the fit is confined to ⁠p0 +/- grad_bounds⁠ on the optimizer scale, which for a log-scale parameter is a factor of exp(grad_bounds) (~148 at the default 5). This bound is admixr2's, not the model's – an unbounded parameter has no other – and nloptr reports normal convergence at a box corner, so a warning is emitted if an estimate finishes on it.

Inner FD step for the gradient-based Hessian (only used when covMethod = "r" and grad != "none"). Default 1e-3.

Outer step scale for NLL-FD Hessian.

"r,s" (the DEFAULT) computes the sandwich ⁠H^-1 J H^-1⁠; "r" the numerical Hessian alone, ⁠2H^-1⁠; "none" skips the covariance. All three span the structural, residual-error and omega parameters. Omega is included because excluding it also biases the STRUCTURAL standard errors downward – a theta carrying an eta is correlated with that eta's variance. If the weakly-identified omega Cholesky makes the Hessian non-positive definite, the structural + residual sub-block is reported with a warning.

"r,s" adds a sandwich correction, ⁠H^-1 J H^-1⁠, on the same Hessian. For FO this does more than correct kurtosis: ⁠V = J Omega J' + Sigma⁠ is the covariance of an exactly normal individual law, so the reported standard errors otherwise answer to the linearisation rather than to the model. The correction scores the FO fit against the model's true nonlinear law, built post-fit on a quadrature ensemble, and so absorbs part of the linearisation error as well. Point estimates are untouched. Transform-both-sides endpoints. adfo composes the residual by a second-order expansion about the linearised moments, because FO carries no node ensemble to compose over – that is what the method is. adgh and admc compose exactly at their nodes/draws, so an adfo fit of a boxCox, yeoJohnson, logitNorm or probitNorm endpoint differs from theirs by the expansion's truncation: roughly 0.3% in V at moderate between-subject variability, rising to ~3% for a tightly-bounded logit/probit at high variability. That is a property of the estimator, not a discrepancy.

It is the default because it is the conservative choice, not the aggressive one. Under correct specification ⁠J = 2H⁠ and the sandwich returns what "r" returns, so defaulting to it costs nothing when the normal-theory assumption holds and corrects the standard errors when it does not. Anything it cannot build degrades to "r" and reports "r", so no fit loses its covariance by asking. Pass covMethod = "r" for the pre-0.4.1 behaviour.

Applies to every residual family whose conditional law is independent across timepoints, which is all of them except ar(): the conditionally-normal set (add, prop, pow, combined1, combined2), the closed-form distributional ones (lnorm, pois, binom, nbinomMu, beta, and t() with nu > 4), and the transform-both-sides ones (boxCox, yeoJohnson, logitNorm, probitNorm), whose third and fourth conditional moments come off the same quadrature that already gives their mean and variance. Refused, and degraded to "r": ar(), because it correlates the residual ACROSS timepoints and the cross terms the expansion drops are then real; t() with nu <= 4, whose kurtosis does not exist; and ordinal() and same-subject joint studies, which stack several outputs into one covariance the per-output node ensemble does not describe. These four are refusals by construction rather than failures, so the fit reports the reason as a message and falls back to "r"; a sandwich that was attempted and could not be built still warns.

"r,s" is more sensitive to an ill-conditioned Hessian than "r" is. "r" reports ⁠2H^-1⁠ and inverts H once; the sandwich reports ⁠H^-1 J H^-1⁠ and inverts it twice, so in a direction the data barely identifies any gap between J and ⁠2H⁠ is amplified quadratically. A residual SD contributing 0.01 variance against 1.7 from between-subject variability is such a direction: measured on one 1-cmt fixture at cond(H) = 3.5e5, the reported residual SE moved by a factor of 0.11 and two omega entries by 0.59 and 1.55, while the same model and design on a study the residual IS identified in (cond(H) = 247) reproduced "r" to four decimals on every parameter. Neither number is a correction there – both methods are reporting an unidentified direction, and "r,s" is louder about it. admixr2 says so: when the Hessian's reciprocal condition number falls below eps^(1/4) – the point at which squaring the conditioning reaches the bound a single inversion is already called singular at – the fit records a note naming the parameter that loads most heavily on the offending direction. It arrives on fit$runInfo and is listed by print(fit), which is where nlmixr2est routes an estimator's warnings. The sandwich is still reported, because the well-determined parameters of the same fit are unaffected; check the named parameter's relative standard error before reading its "r,s" value as a finding.

All three blocks are reported on the scale the ESTIMATES are printed on, as nlmixr2est does: structural thetas on the log/optimizer scale, residual error as an SD, and omega as the variance/covariance entries (named ⁠om.<eta>⁠ and ⁠cov.<eta_i>.<eta_j>⁠). The omega block is rotated by the full Jacobian of Omega with respect to the log-Cholesky, which is not diagonal once omega is correlated.

An adfo standard error describes scatter, not accuracy. FO linearises the model at eta = 0, and on a non-additive residual (or a saturating endpoint, or a large omega) the resulting point estimates carry a bias of several standard errors – measured 5-20 SE, giving 0% coverage for a nominal 95% interval even where the SE itself matches the sampling SD. Use adgh or admc when the uncertainty matters.

Number of optimizer restarts (1 = no multi-start).

Standard deviation for random perturbations of initial struct thetas at each restart (> 1).

Number of parallel workers (mirai daemons) for multi-restart (default 1 = sequential). Requires the mirai package.

rxode2::rxControl() object. Created automatically when NULL.

Passed to nlmixr2est::foceiControl() for the table/output machinery.

Significant digits asked of the ODE solver, or NULL (the default) to leave rxode2's own solver tolerances alone. When set, it is passed to rxode2::rxSolve()'s own sigdig argument for every solve the estimator issues – rxode2 owns the mapping to atol/rtol and has changed it between releases, which is why the digits, not the tolerances, are what travels – and to nlmixr2est::foceiControl() for the post-fit tables.

It is a speed lever, and an opt-in one because it is not free. The estimators finite-difference the solve with steps of the same order: grad_h (1e-4), cov_h (1e-3) and cov_h_outer (~2.5e-3), while sigdig = 4 maps to a relative tolerance of ~1e-4 on current rxode2. Differencing a solution whose own noise is 1e-4 with a 1e-4 step returns noise, and it surfaces as a moved objective and an indefinite covariance Hessian (every SE reported NA) rather than as an error. Most worthwhile where the gradient is fully analytic and nothing differences the solve – adfoControl(grad = "analytical") measured ~4.8x faster at sigdig = 4 with standard errors unchanged to 4 significant figures. Elsewhere, compare the objective and the standard errors against NULL before relying on it. Table formatting is unaffected either way: sigdigTable defaults to 4 regardless.

How combined additive+proportional error is parameterised in the nlmixr2 output tables: "combined2" (default, variance form) or "combined1" (SD form). Has no effect on admixr2's own estimation.

If TRUE, return a plain list instead of the full nlmixr2 fit object.

Gauss-Hermite nodes used to integrate the RESIDUAL for a transform-both-sides endpoint (boxCox, yeoJohnson, logitNorm, probitNorm), where y = g(h(f) + sigma*eps) has no closed-form mean and variance. Ignored by every other error model, which has closed forms. Default 81. Measured worst-case relative error against an independent quadrature, over all four transforms and residual SD of 0.5, 1, 2 and 3: n = 15 gives 5.7e-2, 31 gives 4.5e-3, 81 gives 5.0e-5. The error is dominated by large residual SD; at SD <= 1, n = 31 already gives 1e-7 or better.

This is an ACCURACY dial, not a speed one. The quadrature is linear in resid_nodes in isolation (~50 us at 15, 300 us at 81 for an 8-row study) but negligible beside the ODE solve: a full NLL evaluation measured 0.750 s per 60 evaluations at BOTH 31 and 81 nodes. Raise it if you have a saturating endpoint with a large residual SD; there is little to gain by lowering it.

Unused arguments (trigger an error).

Named list of study specifications (same format as admControl(): E, V, n, times, ev, optional method; or an observations list for multi-compartment fits – see admControl()).

Number of quadrature nodes per eta dimension (default 5). For a transform-both-sides endpoint (boxCox, yeoJohnson, logitNorm, probitNorm) this also controls the accuracy of the RESIDUAL composition: those endpoints have a conditional mean that is nonlinear in the structural prediction, so the residual is composed at each node and aggregated rather than expanded about the ensemble mean. Before that, n_nodes had no effect at all on a TBS fit's accuracy – the expansion's error was a floor no node count removed. Total nodes = n_nodes^n_eta. n_nodes = 5 achieves near-exact covariance moments for IIV SD up to ~0.5; n_nodes = 7 extends coverage to SD ~0.7. For models with >= 5 etas the node count grows steeply; consider reducing n_nodes or using a different estimator.

Gradient mode. "analytical" (default) uses closed-form contractions through the sensitivity equations – cheapest and exact. "fd" uses central finite differences (forward differencing was removed in 0.4.1; see adfoControl()). "none" uses derivative-free BOBYQA.

nloptr algorithm, or NULL (default) to pick the default that matches grad: "NLOPT_LD_LBFGS" with a gradient, "NLOPT_LN_BOBYQA" when grad = "none". Any algorithm reported by nloptr::nloptr.print.options() is accepted. An explicit algorithm is reconciled with grad: when grad = "none" a gradient-based algorithm (⁠NLOPT_LD_*⁠ / ⁠NLOPT_GD_*⁠) falls back to "NLOPT_LN_BOBYQA"; when a gradient is requested a derivative-free algorithm (⁠NLOPT_LN_*⁠ / ⁠NLOPT_GN_*⁠) turns the gradient off. Both emit a message.

Maximum function evaluations (default 500).

Relative tolerance (default sqrt(.Machine$double.eps)).

Print-frequency for live progress (0 = silent).

Random seed (used for restarts).

OpenMP threads for rxSolve(). Defaults to rxode2::rxCores(). When workers > 1 it is a total budget, split across the workers.

Passed to rxSolve(): show the solver's text progress bar only once a single solve exceeds this many subjects. The default (.Machine$integer.max) keeps it off for clean script/vignette output; lower it (e.g. 1000L) to see progress during long fits.

Finite-difference step for unpaired struct theta gradient and FD Jacobian fallback.

Box-constraint half-width when using gradients: the fit is confined to ⁠p0 +/- grad_bounds⁠ on the optimizer scale, which for a log-scale parameter is a factor of exp(grad_bounds) (~148 at the default 5). This bound is admixr2's, not the model's – an unbounded parameter has no other – and nloptr reports normal convergence at a box corner, so a warning is emitted if an estimate finishes on it.

Inner FD step for the gradient-based Hessian (only used when covMethod = "r" and grad != "none").

Outer step scale for numerical Hessian. Default eps^(1/4) (tighter than admc's eps^(1/5) because the GH surface is noise-free).

"r,s" (the DEFAULT) computes the sandwich ⁠H^-1 J H^-1⁠; "r" the numerical Hessian alone, ⁠2H^-1⁠; "none" skips the covariance. All three span the structural, residual-error and omega parameters. Omega is included because excluding it also biases the STRUCTURAL standard errors downward – a theta carrying an eta is correlated with that eta's variance. If the weakly-identified omega Cholesky makes the Hessian non-positive definite, the structural + residual sub-block is reported with a warning.

"r,s" adds a sandwich correction, ⁠H^-1 J H^-1⁠, on the same Hessian. The aggregate objective scores the reported mean and covariance as though the subjects behind them were multivariate normal; they are not, because the model is nonlinear in the random effects, so the sampling law of ⁠(E, V)⁠ is not the one the objective assumes. "r,s" scores that law from the model instead. Point estimates are untouched – only the reported uncertainty changes – and under correct specification it reduces to "r" exactly. It is the default because it is the conservative choice, not the aggressive one. Under correct specification ⁠J = 2H⁠ and the sandwich returns what "r" returns, so defaulting to it costs nothing when the normal-theory assumption holds and corrects the standard errors when it does not. Anything it cannot build degrades to "r" and reports "r", so no fit loses its covariance by asking. Pass covMethod = "r" for the pre-0.4.1 behaviour.

Applies to every residual family whose conditional law is independent across timepoints, which is all of them except ar(): the conditionally-normal set (add, prop, pow, combined1, combined2), the closed-form distributional ones (lnorm, pois, binom, nbinomMu, beta, and t() with nu > 4), and the transform-both-sides ones (boxCox, yeoJohnson, logitNorm, probitNorm), whose third and fourth conditional moments come off the same quadrature that already gives their mean and variance. Refused, and degraded to "r": ar(), because it correlates the residual ACROSS timepoints and the cross terms the expansion drops are then real; t() with nu <= 4, whose kurtosis does not exist; and ordinal() and same-subject joint studies, which stack several outputs into one covariance the per-output node ensemble does not describe. These four are refusals by construction rather than failures, so the fit reports the reason as a message and falls back to "r"; a sandwich that was attempted and could not be built still warns.

"r,s" is more sensitive to an ill-conditioned Hessian than "r" is. "r" reports ⁠2H^-1⁠ and inverts H once; the sandwich reports ⁠H^-1 J H^-1⁠ and inverts it twice, so in a direction the data barely identifies any gap between J and ⁠2H⁠ is amplified quadratically. A residual SD contributing 0.01 variance against 1.7 from between-subject variability is such a direction: measured on one 1-cmt fixture at cond(H) = 3.5e5, the reported residual SE moved by a factor of 0.11 and two omega entries by 0.59 and 1.55, while the same model and design on a study the residual IS identified in (cond(H) = 247) reproduced "r" to four decimals on every parameter. Neither number is a correction there – both methods are reporting an unidentified direction, and "r,s" is louder about it. admixr2 says so: when the Hessian's reciprocal condition number falls below eps^(1/4) – the point at which squaring the conditioning reaches the bound a single inversion is already called singular at – the fit records a note naming the parameter that loads most heavily on the offending direction. It arrives on fit$runInfo and is listed by print(fit), which is where nlmixr2est routes an estimator's warnings. The sandwich is still reported, because the well-determined parameters of the same fit are unaffected; check the named parameter's relative standard error before reading its "r,s" value as a finding.

All three blocks are reported on the scale the ESTIMATES are printed on, as nlmixr2est does: structural thetas on the log/optimizer scale, residual error as an SD, and omega as the variance/covariance entries (named ⁠om.<eta>⁠ and ⁠cov.<eta_i>.<eta_j>⁠). The omega block is rotated by the full Jacobian of Omega with respect to the log-Cholesky, which is not diagonal once omega is correlated.

Number of optimizer restarts (1 = no multi-start).

SD of random perturbations of initial struct thetas at each restart.

Number of parallel workers (mirai daemons) for multi-restart (default 1 = sequential). Requires the mirai package.

rxode2::rxControl() object. Created automatically when NULL.

Passed to nlmixr2est::foceiControl() for the table/output machinery.

Significant digits asked of the ODE solver, or NULL (the default) to leave rxode2's own solver tolerances alone. When set, it is passed to rxode2::rxSolve()'s own sigdig argument for every solve the estimator issues – rxode2 owns the mapping to atol/rtol and has changed it between releases, which is why the digits, not the tolerances, are what travels – and to nlmixr2est::foceiControl() for the post-fit tables.

It is a speed lever, and an opt-in one because it is not free. The estimators finite-difference the solve with steps of the same order: grad_h (1e-4), cov_h (1e-3) and cov_h_outer (~2.5e-3), while sigdig = 4 maps to a relative tolerance of ~1e-4 on current rxode2. Differencing a solution whose own noise is 1e-4 with a 1e-4 step returns noise, and it surfaces as a moved objective and an indefinite covariance Hessian (every SE reported NA) rather than as an error. Most worthwhile where the gradient is fully analytic and nothing differences the solve – adfoControl(grad = "analytical") measured ~4.8x faster at sigdig = 4 with standard errors unchanged to 4 significant figures. Elsewhere, compare the objective and the standard errors against NULL before relying on it. Table formatting is unaffected either way: sigdigTable defaults to 4 regardless.

How combined additive+proportional error is parameterised in the nlmixr2 output tables: "combined2" (default) or "combined1".

If TRUE, return a plain list instead of the full nlmixr2 fit object.

Gauss-Hermite nodes used to integrate the RESIDUAL for a transform-both-sides endpoint (boxCox, yeoJohnson, logitNorm, probitNorm), where y = g(h(f) + sigma*eps) has no closed-form mean and variance. Ignored by every other error model, which has closed forms. Default 81. Measured worst-case relative error against an independent quadrature, over all four transforms and residual SD of 0.5, 1, 2 and 3: n = 15 gives 5.7e-2, 31 gives 4.5e-3, 81 gives 5.0e-5. The error is dominated by large residual SD; at SD <= 1, n = 31 already gives 1e-7 or better.

This is an ACCURACY dial, not a speed one. The quadrature is linear in resid_nodes in isolation (~50 us at 15, 300 us at 81 for an 8-row study) but negligible beside the ODE solve: a full NLL evaluation measured 0.750 s per 60 evaluations at BOTH 31 and 81 nodes. Raise it if you have a saturating endpoint with a large residual SD; there is little to gain by lowering it.

Unused arguments (trigger an error).

Named list of study specifications. Each element is a list with:

• E – observed mean vector

• V – observed covariance matrix or variance vector (auto-detected)

• n – sample size

• times – numeric vector of observation times

• ev – rxode2::et() dosing event table

• method – "cov" or "var" (optional; auto-detected from V)

• v_denom – "ml" (default) or "unbiased", declaring which denominator the supplied V uses. The likelihood is the exact one for n iid draws only under the ML (n) covariance, which is what cov.wt(method = "ML") and datagen() produce. A published SD is the unbiased (n - 1) SD, so a digitised figure gives V = SD^2 on the n - 1 scale: declare v_denom = "unbiased" and admixr2 converts it. Declared per study, since a meta-analysis routinely mixes a digitised source with a model-derived one and the two need not share a denominator. At n = 60 the factor is 1.7%; it matters more the smaller n is, and more again for any method that scores the reported covariance against its own sampling law.

Multi-compartment (multiple observed outputs). To fit several observed compartments simultaneously (e.g. plasma and brain/CSF), give the study an observations list instead of top-level E/V/times. Each entry is one observed output with its own output (the model prediction variable, e.g. "cp" or "cCSF"), times, E, V and – for independent fits – ev and n. Pass the endpoint names to admData(), e.g. admData(c("cp", "cCSF")), so nlmixr2 recognises every endpoint. There are two modes:

• Independent – each observed output has its own n/ev (separate experiments / subjects, e.g. a plasma study and a brain study combined for meta-analysis). The outputs are independent likelihood blocks and the aggregate ⁠-2LL⁠ is their sum.

• Joint (same subjects) – the outputs are measured on the SAME subjects. Give the study a shared n and ev, and a joint covariance either as a study-level full matrix V (blocks in observations order) or as per-output marginal V plus a cross list of cross-covariance blocks keyed "outA:outB" (each length(times_A) x length(times_B); omitted pairs are zero). The compartments are then scored by a single MVN over the stacked vector with shared random effects. est = "adirmc" does not support multiple observed outputs; use "admc", "adfo" or "adgh".

Long format (one row per endpoint/time). As an alternative to the observations list, a study may carry a data frame that keys each observed summary by endpoint, the way nlmixr2 keys observations by DVID/CMT. The frame needs an endpoint column (DVID, CMT or output), a time column (TIME), a mean column (E) and – unless a joint V is given – a variance column (V) or an SD column (SD). It is normalised into exactly the same units as the observations form, so the two are interchangeable:

# independent blocks: per-row variances; optional per-endpoint `n` column
# and per-endpoint `ev` (a list of event tables keyed by endpoint)
list(n = 60L, ev = ev,
     data = data.frame(DVID = c("cp", "cp", "cCSF"), TIME = c(1, 2, 2),
                       E = c(9.1, 7.4, 2.2), V = c(1.2, 0.9, 0.1)))

# joint (same subjects): ONE stacked covariance whose rows/cols align with
# the rows of `data` -- no `cross` blocks to assemble by hand
list(n = 60L, ev = ev, data = data.frame(DVID = ..., TIME = ..., E = ...),
     V = V_joint)

A study-level V (or an explicit joint = TRUE) marks the endpoints as same-subject; without one, each endpoint is an independent likelihood block. Endpoints are stacked in the order they first appear in data.

Number of Monte Carlo samples per NLL evaluation.

Maximum inner optimiser iterations per phase.

Sampling method for eta draws: "sobol" (Sobol, default), "halton" (Halton), "torus" (Kronecker/torus), "lhs" (Latin hypercube), or "rnorm" (iid normal).

nloptr algorithm string, or NULL (default) to pick the default that matches grad: "NLOPT_LD_LBFGS" with a gradient, "NLOPT_LN_BOBYQA" when grad = "none". Any algorithm reported by nloptr::nloptr.print.options() is accepted (e.g. "NLOPT_LD_MMA", "NLOPT_LN_NELDERMEAD"). An explicit algorithm is reconciled with grad: when grad = "none" a gradient-based algorithm (⁠NLOPT_LD_*⁠ / ⁠NLOPT_GD_*⁠) falls back to "NLOPT_LN_BOBYQA"; when a gradient is requested a derivative-free algorithm (⁠NLOPT_LN_*⁠ / ⁠NLOPT_GN_*⁠) turns the gradient off. Both emit a message.

Maximum number of optimizer function evaluations.

Relative function-value tolerance for convergence.

Print progress every this many evaluations (0 = silent).

Inflate proposal Omega by this factor (>= 1).

Random seed for reproducibility.

Number of OpenMP threads for rxSolve(). Defaults to rxode2::rxCores(). rxSolve() parallelises over subjects, so this is the main speed lever for the MC estimators; when workers > 1 it is a total budget, split across the workers.

Passed to rxSolve(): the solver shows its text progress bar only once a single solve exceeds this many subjects. The default (.Machine$integer.max) keeps the bar off, which is what you want for scripts, vignettes and logs; lower it (e.g. 1000L) to see solver progress during long interactive fits.

Gradient mode for the inner optimiser: "analytical" (default, closed-form weight-path gradient), "none" (derivative-free BOBYQA), or "fd" (central finite differences). Note: "sens" is not available for the IRMC estimator.

Kappa correction method for models with non-mu-referenced struct thetas: "exact" (default, re-evaluates population prediction f(theta, 0) via rxSolve at each inner step), "linearized" (precomputes J = df/d(theta) once per outer iteration using f(theta, 0) as baseline — zero rxSolve per inner step), or "linearized_gh" (same linear approximation but baseline and Jacobian use Gauss-Hermite quadrature E_GH[f(theta, eta)] instead of f(theta, 0) — more accurate baseline at any IIV magnitude, still zero rxSolve per inner step).

Number of GH nodes per eta dimension for kappa_method = "linearized_gh" (default 5). Total quadrature points = kappa_n_nodes^n_eta. Ignored for other kappa methods.

Step size for the inner optimiser's finite-difference gradient (grad = "fd"). Defaults to 1e-6, not the 1e-4 the other three controls use: the IRMC inner NLL is deterministic given fixed proposals, so there is no Monte Carlo noise to step over and the truncation-versus-noise balance that sets 1e-4 elsewhere does not apply. The inner step was a hard-coded 1e-6 until it was made to honour grad_h; inheriting the coarser default would have changed the gradient the loop was tuned for.

Inner FD step for the gradient-based Hessian (only used when covMethod = "r" and grad != "none"). Each gradient evaluation has MC noise of order sigma / cov_h; the Hessian divides that noise by the outer step, giving total noise ⁠sigma / (cov_h * cov_h_outer * |p|)⁠. cov_h = 1e-3 balances truncation error and noise amplification. Increase to 1e-2 if the Hessian is non-positive definite.

Outer step scale for the numerical Hessian. The actual step for parameter p is ⁠max(|p|, 0.1) * cov_h_outer⁠. Applied to both the gradient-FD Hessian (grad != "none") and the NLL-FD Hessian (grad = "none"). Default eps^(1/5) (~2.5e-3) is larger than the textbook eps^(1/4) to account for MC noise in NLL and gradient evaluations; empirically it matches the analytical (sensitivity-equation) Hessian ground truth. Increase (e.g. to 5e-3 or 1e-2) if the Hessian is non-positive definite.

Numeric vector of box-constraint half-widths, one per phase. Phases progressively tighten the search region.

Convergence criterion: phase ends when ⁠|approx - exact| < convcrit⁠.

Stop a phase after this many consecutive worsening iterations.

"r,s" (the DEFAULT) computes the sandwich ⁠H^-1 J H^-1⁠; "r" the numerical Hessian alone, ⁠2H^-1⁠; "none" skips the covariance. All three span the structural, residual-error and omega parameters, and are reported on the scale the estimates are printed on. admControl() documents what the sandwich is, why it is the conservative default, and why it is more sensitive than "r" to an ill-conditioned Hessian; the same implementation runs here. What does NOT carry over is the family coverage, because adirmc itself is narrower: the estimator accepts only add, prop, pow, combined1, combined2 and lnorm residuals, and "r,s" applies to all six. The count, beta, transform-both-sides, ordinal and ar() endpoints admControl() lists are refused by est = "adirmc" itself, before any covariance is reached – they are not models whose sandwich degrades to "r" here, they are models this estimator does not fit.

Number of MC samples for the covariance (Hessian) step. More samples reduce MC noise in NLL evaluations. The NLL-based Hessian (grad = "none") uses a central second difference of the NLL with the same Sobol sequence (CRN) at every perturbed point, so noise largely cancels and cov_n_sim = 10000 (default) is sufficient for most models.

Number of optimization restarts. Runs in parallel when workers > 1.

Standard deviation of structural theta perturbations for restart initialisation.

Number of parallel workers for multi-restart. 1 (default) runs restarts sequentially. Values ⁠> 1⁠ run the restarts on a pool of background R processes (mirai daemons), which behaves the same way on every platform. Requires the mirai package. Workers are stopped automatically after the restart phase so all cores are available for the Hessian step; if a fit is interrupted, admStopWorkers() cleans up any survivors.

rxode2::rxControl() object. Created automatically when NULL.

Passed to nlmixr2est::foceiControl() for the table/output machinery.

Significant digits asked of the ODE solver, or NULL (the default) to leave rxode2's own solver tolerances alone. When set, it is passed to rxode2::rxSolve()'s own sigdig argument for every solve the estimator issues – rxode2 owns the mapping to atol/rtol and has changed it between releases, which is why the digits, not the tolerances, are what travels – and to nlmixr2est::foceiControl() for the post-fit tables.

It is a speed lever, and an opt-in one because it is not free. The estimators finite-difference the solve with steps of the same order: grad_h (1e-4), cov_h (1e-3) and cov_h_outer (~2.5e-3), while sigdig = 4 maps to a relative tolerance of ~1e-4 on current rxode2. Differencing a solution whose own noise is 1e-4 with a 1e-4 step returns noise, and it surfaces as a moved objective and an indefinite covariance Hessian (every SE reported NA) rather than as an error. Most worthwhile where the gradient is fully analytic and nothing differences the solve – adfoControl(grad = "analytical") measured ~4.8x faster at sigdig = 4 with standard errors unchanged to 4 significant figures. Elsewhere, compare the objective and the standard errors against NULL before relying on it. Table formatting is unaffected either way: sigdigTable defaults to 4 regardless.

How combined additive+proportional error is parameterised in the nlmixr2 output tables: "combined2" (default, variance form) or "combined1" (SD form). Has no effect on admixr2's own estimation; passed to nlmixr2est::foceiControl() for the table/output machinery only.

If TRUE, return a plain list instead of a full nlmixr2 fit object (useful for debugging).

Gauss-Hermite nodes used to integrate the RESIDUAL for a transform-both-sides endpoint (boxCox, yeoJohnson, logitNorm, probitNorm), where y = g(h(f) + sigma*eps) has no closed-form mean and variance. Ignored by every other error model, which has closed forms. Default 81. Measured worst-case relative error against an independent quadrature, over all four transforms and residual SD of 0.5, 1, 2 and 3: n = 15 gives 5.7e-2, 31 gives 4.5e-3, 81 gives 5.0e-5. The error is dominated by large residual SD; at SD <= 1, n = 31 already gives 1e-7 or better.

This is an ACCURACY dial, not a speed one. The quadrature is linear in resid_nodes in isolation (~50 us at 15, 300 us at 81 for an 8-row study) but negligible beside the ODE solve: a full NLL evaluation measured 0.750 s per 60 evaluations at BOTH 31 and 81 nodes. Raise it if you have a saturating endpoint with a large residual SD; there is little to gain by lowering it.

Additional arguments (none allowed; triggers an error).

Named list of study specifications. Each element is a list with:

• E – observed mean vector

• V – observed covariance matrix or variance vector (auto-detected)

• n – sample size

• times – numeric vector of observation times

• ev – rxode2::et() dosing event table

• method – "cov" or "var" (optional; auto-detected from V)

• v_denom – "ml" (default) or "unbiased", declaring which denominator the supplied V uses. The likelihood is the exact one for n iid draws only under the ML (n) covariance, which is what cov.wt(method = "ML") and datagen() produce. A published SD is the unbiased (n - 1) SD, so a digitised figure gives V = SD^2 on the n - 1 scale: declare v_denom = "unbiased" and admixr2 converts it. Declared per study, since a meta-analysis routinely mixes a digitised source with a model-derived one and the two need not share a denominator. At n = 60 the factor is 1.7%; it matters more the smaller n is, and more again for any method that scores the reported covariance against its own sampling law.

Multi-compartment (multiple observed outputs). To fit several observed compartments simultaneously (e.g. plasma and brain/CSF), give the study an observations list instead of top-level E/V/times. Each entry is one observed output with its own output (the model prediction variable, e.g. "cp" or "cCSF"), times, E, V and – for independent fits – ev and n. Pass the endpoint names to admData(), e.g. admData(c("cp", "cCSF")), so nlmixr2 recognises every endpoint. There are two modes:

• Independent – each observed output has its own n/ev (separate experiments / subjects, e.g. a plasma study and a brain study combined for meta-analysis). The outputs are independent likelihood blocks and the aggregate ⁠-2LL⁠ is their sum.

• Joint (same subjects) – the outputs are measured on the SAME subjects. Give the study a shared n and ev, and a joint covariance either as a study-level full matrix V (blocks in observations order) or as per-output marginal V plus a cross list of cross-covariance blocks keyed "outA:outB" (each length(times_A) x length(times_B); omitted pairs are zero). The compartments are then scored by a single MVN over the stacked vector with shared random effects. est = "adirmc" does not support multiple observed outputs; use "admc", "adfo" or "adgh".

Long format (one row per endpoint/time). As an alternative to the observations list, a study may carry a data frame that keys each observed summary by endpoint, the way nlmixr2 keys observations by DVID/CMT. The frame needs an endpoint column (DVID, CMT or output), a time column (TIME), a mean column (E) and – unless a joint V is given – a variance column (V) or an SD column (SD). It is normalised into exactly the same units as the observations form, so the two are interchangeable:

# independent blocks: per-row variances; optional per-endpoint `n` column
# and per-endpoint `ev` (a list of event tables keyed by endpoint)
list(n = 60L, ev = ev,
     data = data.frame(DVID = c("cp", "cp", "cCSF"), TIME = c(1, 2, 2),
                       E = c(9.1, 7.4, 2.2), V = c(1.2, 0.9, 0.1)))

# joint (same subjects): ONE stacked covariance whose rows/cols align with
# the rows of `data` -- no `cross` blocks to assemble by hand
list(n = 60L, ev = ev, data = data.frame(DVID = ..., TIME = ..., E = ...),
     V = V_joint)

A study-level V (or an explicit joint = TRUE) marks the endpoints as same-subject; without one, each endpoint is an independent likelihood block. Endpoints are stacked in the order they first appear in data.

Number of Monte Carlo samples per NLL evaluation.

Sampling method for eta draws: "sobol" (Sobol, default), "halton" (Halton), "torus" (Kronecker/torus), "lhs" (Latin hypercube), or "rnorm" (iid normal).

nloptr algorithm string, or NULL (default) to pick the default that matches grad: "NLOPT_LD_LBFGS" with a gradient, "NLOPT_LN_BOBYQA" when grad = "none". Any algorithm reported by nloptr::nloptr.print.options() is accepted (e.g. "NLOPT_LD_MMA", "NLOPT_LN_NELDERMEAD"). An explicit algorithm is reconciled with grad: when grad = "none" a gradient-based algorithm (⁠NLOPT_LD_*⁠ / ⁠NLOPT_GD_*⁠) falls back to "NLOPT_LN_BOBYQA"; when a gradient is requested a derivative-free algorithm (⁠NLOPT_LN_*⁠ / ⁠NLOPT_GN_*⁠) turns the gradient off. Both emit a message.

Maximum number of optimizer function evaluations.

Relative function-value tolerance for convergence.

Print progress every this many evaluations (0 = silent).

Random seed for reproducibility.

Number of OpenMP threads for rxSolve(). Defaults to rxode2::rxCores(). rxSolve() parallelises over subjects, so this is the main speed lever for the MC estimators; when workers > 1 it is a total budget, split across the workers.

Passed to rxSolve(): the solver shows its text progress bar only once a single solve exceeds this many subjects. The default (.Machine$integer.max) keeps the bar off, which is what you want for scripts, vignettes and logs; lower it (e.g. 1000L) to see solver progress during long interactive fits.

Gradient mode: "sens" (sensitivity equations, default), "fd" (central finite differences; forward was removed in 0.4.1), or "none" (derivative-free). A warning is issued when "sens" is requested but the sensitivity model is unavailable; the estimator then falls back to central finite differences.

Step size for finite-difference gradient evaluation during optimization (used by grad = "fd"). This is the FALLBACK step: the step is normally measured per parameter by the Shi (2021) procedure, and grad_h is what a parameter falls back to when that measurement cannot be made (a direction the objective is flat in, or a failed noise estimate).

Inner FD step for the gradient-based Hessian (only used when covMethod = "r" and grad != "none"). Each gradient evaluation has MC noise of order sigma / cov_h; the Hessian divides that noise by the outer step, giving total noise ⁠sigma / (cov_h * cov_h_outer * |p|)⁠. cov_h = 1e-3 balances truncation error and noise amplification. Increase to 1e-2 if the Hessian is non-positive definite.

Outer step scale for the numerical Hessian. The actual step for parameter p is ⁠max(|p|, 0.1) * cov_h_outer⁠. Applied to both the gradient-FD Hessian (grad != "none") and the NLL-FD Hessian (grad = "none"). Default eps^(1/5) (~2.5e-3) is larger than the textbook eps^(1/4) to account for MC noise in NLL and gradient evaluations; empirically it matches the analytical (sensitivity-equation) Hessian ground truth. Increase (e.g. to 5e-3 or 1e-2) if the Hessian is non-positive definite.

Box-constraint half-width when using gradients: the fit is confined to ⁠p0 +/- grad_bounds⁠ on the optimizer scale, which for a log-scale parameter is a factor of exp(grad_bounds) (~148 at the default 5). This bound is admixr2's, not the model's – an unbounded parameter has no other – and nloptr reports normal convergence at a box corner, so a warning is emitted if an estimate finishes on it.

"r,s" (the DEFAULT) computes the sandwich ⁠H^-1 J H^-1⁠; "r" the numerical Hessian alone, ⁠2H^-1⁠; "none" skips the covariance. All three span the structural, residual-error and omega parameters. Omega is included because excluding it also biases the STRUCTURAL standard errors downward – a theta carrying an eta is correlated with that eta's variance. If the weakly-identified omega Cholesky makes the Hessian non-positive definite, the structural + residual sub-block is reported with a warning.

"r,s" adds a sandwich correction, ⁠H^-1 J H^-1⁠, on the same Hessian. The aggregate objective scores the reported mean and covariance as though the subjects behind them were multivariate normal; they are not, because the model is nonlinear in the random effects. "r,s" scores that law from the model instead, on a quadrature ensemble rather than on this fit's own MC draws, so the reported uncertainty carries no sampling noise of its own. Point estimates are untouched, and under correct specification it reduces to "r" exactly. It is the default because it is the conservative choice, not the aggressive one. Under correct specification ⁠J = 2H⁠ and the sandwich returns what "r" returns, so defaulting to it costs nothing when the normal-theory assumption holds and corrects the standard errors when it does not. Anything it cannot build degrades to "r" and reports "r", so no fit loses its covariance by asking. Pass covMethod = "r" for the pre-0.4.1 behaviour.

Applies to every residual family whose conditional law is independent across timepoints, which is all of them except ar(): the conditionally-normal set (add, prop, pow, combined1, combined2), the closed-form distributional ones (lnorm, pois, binom, nbinomMu, beta, and t() with nu > 4), and the transform-both-sides ones (boxCox, yeoJohnson, logitNorm, probitNorm), whose third and fourth conditional moments come off the same quadrature that already gives their mean and variance. Refused, and degraded to "r": ar(), because it correlates the residual ACROSS timepoints and the cross terms the expansion drops are then real; t() with nu <= 4, whose kurtosis does not exist; and ordinal() and same-subject joint studies, which stack several outputs into one covariance the per-output node ensemble does not describe. These four are refusals by construction rather than failures, so the fit reports the reason as a message and falls back to "r"; a sandwich that was attempted and could not be built still warns.

"r,s" is more sensitive to an ill-conditioned Hessian than "r" is. "r" reports ⁠2H^-1⁠ and inverts H once; the sandwich reports ⁠H^-1 J H^-1⁠ and inverts it twice, so in a direction the data barely identifies any gap between J and ⁠2H⁠ is amplified quadratically. A residual SD contributing 0.01 variance against 1.7 from between-subject variability is such a direction: measured on one 1-cmt fixture at cond(H) = 3.5e5, the reported residual SE moved by a factor of 0.11 and two omega entries by 0.59 and 1.55, while the same model and design on a study the residual IS identified in (cond(H) = 247) reproduced "r" to four decimals on every parameter. Neither number is a correction there – both methods are reporting an unidentified direction, and "r,s" is louder about it. admixr2 says so: when the Hessian's reciprocal condition number falls below eps^(1/4) – the point at which squaring the conditioning reaches the bound a single inversion is already called singular at – the fit records a note naming the parameter that loads most heavily on the offending direction. It arrives on fit$runInfo and is listed by print(fit), which is where nlmixr2est routes an estimator's warnings. The sandwich is still reported, because the well-determined parameters of the same fit are unaffected; check the named parameter's relative standard error before reading its "r,s" value as a finding.

All three blocks are reported on the scale the ESTIMATES are printed on, as nlmixr2est does: structural thetas on the log/optimizer scale, residual error as an SD, and omega as the variance/covariance entries (named ⁠om.<eta>⁠ and ⁠cov.<eta_i>.<eta_j>⁠). The omega block is rotated by the full Jacobian of Omega with respect to the log-Cholesky, which is not diagonal once omega is correlated.

Number of MC samples for the covariance (Hessian) step. More samples reduce MC noise in NLL evaluations. The NLL-based Hessian (grad = "none") uses a central second difference of the NLL with the same Sobol sequence (CRN) at every perturbed point, so noise largely cancels and cov_n_sim = 10000 (default) is sufficient for most models.

Number of optimization restarts. Runs in parallel when workers > 1.

Standard deviation of structural theta perturbations for restart initialisation.

Number of parallel workers for multi-restart. 1 (default) runs restarts sequentially. Values ⁠> 1⁠ run the restarts on a pool of background R processes (mirai daemons), which behaves the same way on every platform. Requires the mirai package. Workers are stopped automatically after the restart phase so all cores are available for the Hessian step; if a fit is interrupted, admStopWorkers() cleans up any survivors.

rxode2::rxControl() object. Created automatically when NULL.

Passed to nlmixr2est::foceiControl() for the table/output machinery.

Significant digits asked of the ODE solver, or NULL (the default) to leave rxode2's own solver tolerances alone. When set, it is passed to rxode2::rxSolve()'s own sigdig argument for every solve the estimator issues – rxode2 owns the mapping to atol/rtol and has changed it between releases, which is why the digits, not the tolerances, are what travels – and to nlmixr2est::foceiControl() for the post-fit tables.

It is a speed lever, and an opt-in one because it is not free. The estimators finite-difference the solve with steps of the same order: grad_h (1e-4), cov_h (1e-3) and cov_h_outer (~2.5e-3), while sigdig = 4 maps to a relative tolerance of ~1e-4 on current rxode2. Differencing a solution whose own noise is 1e-4 with a 1e-4 step returns noise, and it surfaces as a moved objective and an indefinite covariance Hessian (every SE reported NA) rather than as an error. Most worthwhile where the gradient is fully analytic and nothing differences the solve – adfoControl(grad = "analytical") measured ~4.8x faster at sigdig = 4 with standard errors unchanged to 4 significant figures. Elsewhere, compare the objective and the standard errors against NULL before relying on it. Table formatting is unaffected either way: sigdigTable defaults to 4 regardless.

How combined additive+proportional error is parameterised in the nlmixr2 output tables: "combined2" (default, variance form) or "combined1" (SD form). Has no effect on admixr2's own estimation; passed to nlmixr2est::foceiControl() for the table/output machinery only.

If TRUE, return a plain list instead of a full nlmixr2 fit object (useful for debugging).

Gauss-Hermite nodes used to integrate the RESIDUAL for a transform-both-sides endpoint (boxCox, yeoJohnson, logitNorm, probitNorm), where y = g(h(f) + sigma*eps) has no closed-form mean and variance. Ignored by every other error model, which has closed forms. Default 81. Measured worst-case relative error against an independent quadrature, over all four transforms and residual SD of 0.5, 1, 2 and 3: n = 15 gives 5.7e-2, 31 gives 4.5e-3, 81 gives 5.0e-5. The error is dominated by large residual SD; at SD <= 1, n = 31 already gives 1e-7 or better.

This is an ACCURACY dial, not a speed one. The quadrature is linear in resid_nodes in isolation (~50 us at 15, 300 us at 81 for an 8-row study) but negligible beside the ODE solve: a full NLL evaluation measured 0.750 s per 60 evaluations at BOTH 31 and 81 nodes. Raise it if you have a saturating endpoint with a large residual SD; there is little to gain by lowering it.

Additional arguments (none allowed; triggers an error).

Optional character vector of observed output (endpoint) names for a multi-compartment model with several prediction lines (e.g. c("cp", "cCSF")). One observation row is emitted per endpoint, keyed by name in the DVID column (nlmixr2's endpoint identifier; CMT is NA on those rows), so nlmixr2's data translation recognises every endpoint. These rows carry a non-NA placeholder DV (1) because nlmixr2's multi-endpoint translator rejects an all-NA-DV dataset; the placeholder is purely for dispatch and never enters the reported objective (each estimator overwrites it with its own aggregate -2LL). When NULL (default) the single-endpoint dummy frame is returned unchanged.

An admFit.

Further admFits to compare it with.

A named list of study specifications. Each element is a list with:

model

An nlmixr2-style model function with ini() and model() blocks. Serves as the data-generating model for this study. May differ between studies. Can be omitted if a top-level default is supplied via the model argument.

times

Numeric vector of observation times.

ev

A dosing event table created with rxode2::et().

n

(Optional) integer sample size; stored as metadata and used when supplying the result to admControl().

observations

(Optional) a named list to generate data for several observed outputs (multi-compartment). Each entry gives one output's output (model prediction variable, e.g. "cp"), times, and optionally ev/n (inherited from the study otherwise). When present, the study result carries a matching observations list of per-output E/V, ready to pass straight to admControl(studies = ...).

Optional default model function used for any study that does not supply its own model element. At least one of model or each study's model must be non-NULL.

A datagenControl() object.

Moment approximation used to generate E and V: "mc" (default) draws Monte Carlo samples over the IIV distribution, as in est = "admc"; "fo" uses the deterministic First-Order expansion (mu = f(theta, 0), ⁠V = J Omega J' + Sigma⁠), matching est = "adfo"; "gh" uses deterministic Gauss-Hermite quadrature over the random-effects prior, matching est = "adgh" – unbiased at any IIV magnitude and noise-free. Use "fo" or "gh" for design evaluation where the data-generating and data-analytic models must coincide.

Number of Monte Carlo samples used to approximate population moments. Ignored when method = "fo" or "gh".

Number of Gauss-Hermite nodes per eta dimension for method = "gh" (default 5). Total nodes = n_nodes^n_eta. Ignored for "mc" and "fo".

Quasi-random sampling method: "sobol" (default), "halton", "torus", "lhs", or "rnorm". Ignored when method = "fo" or "gh".

Integer seed. Applied before stochastic methods ("rnorm", "lhs"). Ignored when method = "fo" or "gh".

Number of rxSolve threads.

Include the raw ⁠n_sim x length(times)⁠ prediction matrix as ⁠$samples⁠ in each study's output. No effect when method = "fo" or "gh" (those methods draw no samples).

Gauss-Hermite nodes used to integrate the RESIDUAL for a transform-both-sides endpoint (boxCox, yeoJohnson, logitNorm, probitNorm), where y = g(h(f) + sigma*eps) has no closed-form mean and variance. Ignored by every other error model. Default 81 – the same default the four estimator controls use, so datagen() and the fit it feeds agree unless you deliberately change one of them. See admControl() for the measured convergence.

nlmixr2 environment containing ui and control.

Unused.

nlmixr2 environment containing ui and control.

Unused.

nlmixr2 environment containing ui and control.

Unused.

nlmixr2 environment containing ui and control.

Unused.

An admFit object returned by nlmixr2() with est = "adfo", est = "admc", est = "adgh", or est = "adirmc".

Character vector selecting which panel types to produce. Any subset of c("mean", "cov", "nll", "par"). Defaults to all four.

Number of MC samples for the final prediction. Defaults to the value used during fitting. Only used when "mean" or "cov" is in which.

Random seed for reproducibility.

Unused.

An admFit object.

Passed to print.nlmixr2FitCore.