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. A study generated from a published model defaults to "none" and refuses an explicit covariance method because it has no sampling law. 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.

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

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. A study generated from a published model defaults to "none" and refuses an explicit covariance method because it has no sampling law. 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.

Gauss-Hermite nodes per covariate used to integrate the COVARIATE distribution when a study declares cov_dist (default 7). This is a separate dial from n_nodes, which refines the random-effect dimensions only: raising n_nodes alone leaves the covariate integration exactly where it was. Measured on a two-compartment model with an allometric weight effect and a lognormal weight distribution, 7 nodes place the marginal moments within 2e-06 (mean) and 2e-05 (covariance) of an exact reference, and the remaining error is the ODE solver's rather than the quadrature's. A wider or more skewed covariate distribution, or a more strongly non-linear covariate effect, warrants more. Measured against an exact reference on a two-compartment model with an allometric weight effect and a lognormal weight distribution: 3 nodes give 7.3e-04 / 8.2e-03 (mean / covariance), 5 give 2.8e-05 / 3.7e-04, 7 give 2.2e-06 / 2.4e-05, and 9 onwards sit at ~1.2e-06 / ~1.0e-06, which is the ODE solver's accuracy rather than the quadrature's. The default is set past that knee, and raising it further buys nothing: against a per-subject reference the accuracy is identical at 5, 9 and 15 nodes. Ignored when cov_integration = "sparse", which sets its own resolution through cov_sparse_level.

It is a nodes-per-DIRECTION budget rather than a literal node count. Where the covariates reach the model through fewer scalars than there are covariates, admixr2 integrates over those directions instead of over a product grid, and each direction is given cov_nodes * p / r nodes rounded up – MORE than cov_nodes, because a direction that absorbs several covariate axes carries their combined spread and needs proportionally more resolution to resolve it. Three covariates reaching the model as a single scalar therefore get 21 nodes on one direction at the default, not 7, and still cost 21 design points against the product grid's 343. The same budget sizes the directions of a joint random-effect/covariate design where one is used.

How a study's covariate distribution is integrated. Three states.

"on" (default) integrates on a product Gauss–Hermite grid of cov_nodes points per covariate and reduces it wherever the model permits, choosing per study without being asked. Where the random effects and the covariates span fewer directions than they have members – an allometric weight effect on the same parameter as its random effect is one direction, not two – the integral is taken over those directions instead, which for p covariates costs design points in the rank rather than cov_nodes^p. A study that does not qualify is integrated on the full grid. There is nothing to tune: a reduction is admitted only after it reproduces the design it stands in for, so it cannot trade accuracy for speed behind your back. Measured across four model shapes it is 2.5x to 17x cheaper AND 100x to 170000x more accurate than the unreduced grid.

"off" disables every reduction and integrates on the full product grid. Slower, and useful mainly as a reference when a result is in question.

"sparse" replaces the product rule with a Smolyak sparse grid of cov_sparse_level, which for p covariates costs far fewer than cov_nodes^p points and is the speed lever for models with several covariates.

Both are Gauss–Hermite rules; they differ in which product terms are kept. Measured against an exact reference (lognormal margins, allometric plus a saturable term), relative error on the mean and the covariance:

At four covariates level 3 is both cheaper than the 3-node product grid and roughly 40x more accurate, and the advantage grows with p. Level 2 is the axial rule — at one covariate it is exactly cov_nodes = 3 — and it is offered for continuity rather than recommended.

DEPENDENT covariates (cor, rho, Sigma) are handled by rotating onto the eigenvectors of the latent correlation, and correlation does not cost the sparse rule accuracy: at p = 2 its mean error is 6.6e-07 at rho = 0 and 4.8e-08 at rho = 0.85. (Level 2 behaves the other way, losing an order of magnitude to correlation, which is one reason the default is 3.) An opaque joint sampler is refused, because the rotation needs a correlation the closure does not report — declare the dependence with cor and admixr2 builds the sampler itself.

The cost of a sparse rule is SIGNED weights: they sum to 1 exactly, but the sum of their absolute values is 2.5 at level 3 for three covariates and 4.1 for four, so the answer is a difference of terms several times its own size and solver noise is amplified accordingly. A sandwich covariance whose weight matrix comes out indefinite as a result is refused rather than reported.

Smolyak level for cov_integration = "sparse" (default 3, minimum 2). Level 2 is the axial rule, level 3 adds the five-point axes and the pairwise crosses, and each further level refines again at a growing weight-magnitude cost. See cov_integration for the measured accuracy and point counts.

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

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. A study generated from a published model defaults to "none" and refuses an explicit covariance method because it has no sampling law. 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.

Relative parameter tolerance for convergence (default sqrt(.Machine$double.eps)).

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. A study generated from a published model defaults to "none" and refuses an explicit covariance method because it has no sampling law. 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.

Relative parameter tolerance for convergence (default sqrt(.Machine$double.eps)).

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 object.

Simulation size and seed for the predicted moments. Defaults to the fit's own.

"source" to collapse a conditional source back together, "stratum" to keep its strata apart.

Named covariates. A continuous one takes a named vector, e.g. WT = c(mean = 75, sd = 16) or CRCL = c(median = 92, iqr = c(62, 118)). A binary one takes a single named proportion, e.g. SEX = c(male = 0.55), which becomes levels 0/1 with that probability on 1.

Correlations between covariate PAIRS, named A.B, e.g. cor = c(WT.CRCL = 0.45). Pairs not named are independent, so a partial table needs no identity padding. A full matrix is accepted too.

"lnorm" (default) or "normal", for the continuous margins. Lognormal is the usual choice for a positive covariate — a normal margin wide enough to matter puts mass at or below zero, which is NaN inside any power or log term.

A data frame of individual covariates to derive the table FROM, instead of typing it out — a digitised baseline listing, or the cohort itself in a simulation study. Each numeric column becomes a margin (a 0/1 column becomes a proportion, everything else a continuous margin matching that column's mean and SD), and every continuous PAIR gets its correlation — taken on the LATENT scale, so on the logs for a lognormal margin, which is the conversion easiest to get wrong by hand. Anything named in ... or cor overrides what the data would have given, so a column you would rather state yourself simply gets stated.

admStudy() objects, optionally named.

The published model, as an nlmixr2-style function (or a parsed rxUi). Omit for a digitised study.

Named vector of the paper's parameter estimates, on the scale the model's ini() is written on — so tcl = log(5.2) for a clearance reported as 5.2. Omitted parameters keep the model's own ini() value.

Digitised aggregate data: the reported mean profile E, with its spread as sd (per timepoint), sem (converted using n), or a full covariance V.

Number of subjects the study reports on. Taken from population when that is the cohort's data frame.

Observation times.

Dose amount, as shorthand for a single-dose ev.

A dosing event table from rxode2::et(), for anything dose cannot express.

The population the study enrolled — see admPopulation(). Covariates the model reads are integrated over it. A data frame of enrolled subjects is accepted directly and read into a table, taking n from its row count. EVERY column that is not a reserved data column (ID, TIME, DV, AMT, ...) becomes a covariate, including ones this model does not read — a covariate another source fitted is still evidence here, and each one adds a quadrature axis, so drop the columns you do not want integrated over.

Named list pinning a covariate at a value, for a study reported in one subgroup, e.g. at = list(SEX = 1). A pinned covariate must be omitted from population.

Covariate the paper reports results SEPARATELY by, e.g. by = "SEX". Expands a published model into one study per level; digitised subgroup profiles must be supplied as separate studies.

Quadrature resolution for a conditional CONTINUOUS covariate, i.e. how many nodes it is cut into. Default 9. Raise it for precision, at a multiplicative cost: a source conditional on two continuous covariates becomes strata_nodes^2 studies, and each is solved. Discrete covariates are unaffected — their levels are exact and cost one stratum each. Note that the objective VALUE depends on this, so two fits are only comparable at the same resolution; admixr2 records it per study and anova() checks it.

Where the nodes are cut on the span of the source and analysis models rather than the product grid, two error sources remain and neither is reduced by raising strata_nodes. A range puts the nodes on a fine cloud, which carries a quadrature error of its own, and the recombination that thins that cloud preserves a set of polynomial moments rather than the integrand the objective contains. strata_nodes refines the moment matching against the same cloud. admixr2 sizes and certifies the cloud at admission — against a coarser cloud, and against nonlinear probes the recombination did not match by construction — and keeps the product grid when it cannot. Untruncated margins reproduce the product objective to ~1e-08; a truncated one to ~3e-05.

Optional named list giving the covariate span the source ENROLLED, e.g. range = list(WT = c(52, 118)). The declared distribution is truncated to it, whether that covariate ends up conditional (strata are then cut from the truncated margin) or marginal (the truncated margin is integrated over). Needed when that distribution is wider than the enrolment — a ⁠mean +/- SD⁠ transcribed from a baseline table has tails where nobody was — and not wanted when population is the patients themselves, since the margins are then fitted to exactly who was enrolled and truncating would score the paper against a sub-population of its own. An unnamed value is allowed only when the source itself is conditional on exactly one covariate — which covariate an unnamed range belongs to is a question about the source, so the answer does not change with the model being fitted to it.

Optional display name; otherwise taken from the argument name in admStudies().

Which denominator the supplied spread uses: "unbiased" (n - 1) or "ml" (n). Usually leave it unset — the currency you wrote the study in already says which it is, and admStudy() resolves it and shows the answer when the study is printed.

sd/sem default to "unbiased", because a published spread is the n - 1 one; V defaults to "ml", because handing over a covariance matrix is the deliberate act of someone who computed it, and the likelihood wants the ML denominator. Set it explicitly when your source breaks that pattern — a sd you computed yourself with the ML denominator, say. A published model is always "ml": its moments are generated at that denominator, so "unbiased" is refused rather than applied to a V that is already right.

An admFit.

Further admFits to compare it with.

One argument per covariate, named as the model reads it. Each is one of:

c(mean = , sd = )

Mean and SD on the covariate's own scale — what a baseline table reports. Normal by default; dist = "lnorm" moment-matches a lognormal margin instead, which is what you want for a covariate that must stay positive (see dist).

c(mu = , sd = ) / c(meanlog = , sdlog = )

A normal or lognormal margin given directly on its own scale.

c(label = prob, ...)

A categorical covariate: names are the level labels, values their proportions. Levels are coded ⁠0, 1, ...⁠ in the order given.

a list(...)

The canonical form, for anything else — including list(quantile = f) for an arbitrary margin.

Alternatively a single data.frame with covariate, mean and sd columns (and optionally dist): a baseline-characteristics table transcribed as-is.

Correlation between the covariates: a scalar for two of them, or a correlation matrix (named, in any order). Realised through a Gaussian copula on the declared margins.

Not accepted yet, and an error if given. It is the place an arbitrary sampler will attach — a function receiving the matrix of uniforms admixr2 supplies and returning one named column per covariate, on each covariate's own scale, which is the shape a vine copula produces — but the contract around it is not settled, so the argument refuses rather than half-works. Use cor for dependence: it builds a Gaussian copula over the declared margins.

Default margin for the c(mean = , sd = ) form: "normal" (default) or "lnorm". A per-covariate dist wins over it.

Choose "lnorm" when the model needs the covariate to stay positive. A normal margin is unbounded below and the quadrature reaches 3.75 standard deviations, so any covariate with a coefficient of variation above about 0.27 (the guard is mu - 3.75 * sd <= 0) gets a node at or below zero, where an allometric or log term is NaN. covDist() warns when that would happen.

A covariate specification, as given to a study. Each element names a covariate and describes its distribution, either as mean and sd on the covariate's own scale — NORMAL by default, so pass dist = "lnorm" for the positive margin an allometric term needs (only covDist() defaults to lognormal) — or as a quantile function, or as values (with optional probs) for a discrete covariate. A cor entry — a scalar for two covariates, or a correlation matrix — links them through a Gaussian copula. A joint function takes the matrix of uniforms and returns one column per covariate, which is how an arbitrary copula, an R-vine included, is supplied.

Number of subjects to draw.

Number of random effects in the model the specification belongs to. The draws are deterministic and come from Sobol dimensions after the random effects', so passing the model's own n_eta reproduces exactly the values a fit would use. The default of 0 is right for inspecting a specification on its own.

A covariate specification, as given to a study — see covDraw() for the full grammar.

Character vector naming the covariates to stratify on. The rest are left in each stratum's own cov_dist, to be marginalised over their distribution conditional on that stratum.

Strata per stratified covariate (default 9). A discrete covariate ignores it and is cut at its levels, exactly, with the level probabilities as weights.

For a continuous one it sets how finely the covariate is resolved, and the rule depends on what is known about the joint distribution. Where the latent structure is known — an explicit cor, or no declared dependence, which IS independence — each stratum is a Gauss-Hermite node, held at a POINT, and the remaining covariates follow their exact conditional law within it. An opaque joint sampler cannot be conditioned, and there the covariate is cut into that many equiprobable BINS instead, each described by the pool members that landed in it.

This is a convergence parameter, not a modelling choice. The well-defined object is the limit as the count grows; any finite value approximates it, and under misspecification the answer can jump between basins rather than drift. Raise it until the estimates stop moving rather than picking a value you like. Cost is n_nodes^p studies for p stratified covariates, so this is cheap in one covariate and expensive in three.

How much the objective depends on it turns on which rule applies. Where the latent structure is known — an explicit cor, or no declared dependence, which is independence — conditioning uses a Gauss-Hermite grid and the objective is stable: 0.03 units across counts of 5 to 100, so objective, AIC, BIC and likelihood ratios are comparable across resolutions.

A declared discrete covariate no longer forces the slow rule: as long as it is latently independent of the others it enumerates exactly at its levels, which took a study declaring one sex covariate from 515 units of movement, still rising at a count of 50, to 0.009.

An opaque joint sampler cannot be conditioned and falls back to equiprobable bins, where it is not stable: 51 units at a count of 100 and still moving, with the estimate itself wandering (0.681 at 5, 9, 15 and 100; 0.697 at 50). There, raise the count until the estimates settle, and note that anova() refuses to compare two fits built at different ones.

Total sample size to divide among the strata. The default of 1 returns the raw stratum weights.

Minimum size of the deterministic pool the strata are cut from. Each stratum's covariate distribution is that pool restricted to its own bin, so this sets how finely the within-stratum distribution is resolved — and, because no stratum can produce a covariate value beyond the pool's extremes, how far into the tails it reaches. The pool is grown automatically with the number of strata, and a stratum left with too few draws is an error rather than a silently coarse answer; n_pool_cell in the result reports what each one got.

It applies only where a pool is used at all: an opaque joint sampler, or a stratified covariate latently correlated with a marginalised one. Where the conditional law is closed form the stratum carries the DECLARED specs themselves, no sample stands in for them, and n_pool_cell is NA.

Optional named list giving the range each stratified covariate was actually ENROLLED over, e.g. list(WT = c(52, 118)) — publications routinely report a min-max or a median with an IQR.

Without it the strata are cut over the whole declared distribution, and a source is credited with evidence at covariate values it never sampled. The overstatement is exactly var(declared) / var(enrolled): 3.4x for a source spanning ±1 SD, 12.4x at ±0.5 SD. Estimates stay correct — what inflates is confidence, and it only shows once a second source disagrees. Omitting it warns for that reason.

More strata resolve the covariate range more finely but do not buy accuracy, and each one costs a solve: on a matched one-covariate fit the coefficient came back at 0.7000 / 0.7002 / 0.7005 / 0.7005 for 3 / 4 / 10 / 16 strata against a true 0.700, while the fit took 7.2 / 4.1 / 8.0 / 11.5 seconds. Raise it when the covariate effect is strongly nonlinear over the range; the only hard limit is that a stratum must keep enough pooled draws to describe itself, which admixr2 checks.

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().

cov_dist

(Optional) the covariate distribution this study's subjects span — see covDraw() for the grammar. The generated E/V are MARGINAL over it, which is what a publication reports. Needs datagenControl(method = "mc") or "gh"; "fo" integrates over the random effects only and is refused. Prefer "gh", which adds no Monte Carlo noise to data that is meant to BE the reference.

stratify

(Optional) TRUE to stratify on every covariate this study's OWN data-generating model ESTIMATED a coefficient for, marginalising the rest — the split is read from the model, so it cannot disagree with it. Two sources sharing one cov_dist therefore stratify differently, each according to what it fitted. A covariate the model merely READS — weight at a fixed allometric exponent — is refused: there is no fitted effect to recover there. A character vector names the covariates explicitly instead, and is taken at its word. NOTHING IS DERIVED WITHOUT THIS: a study that asks for no stratification is generated MARGINAL over whatever it declared, which is what a publication reports. The study is expanded into one ordinary study per covariate stratum, named ⁠<study>_s1⁠, ⁠<study>_s2⁠, ..., each pinned at its own covariate value, carrying its own effective size n_k (the quadrature weight times n, summing to n), and marginalising the remaining covariates over their distribution CONDITIONAL on that stratum. Use it when the source reports — or, being a published model, can report — summaries by covariate subgroup. Stratify only on what the source actually fitted: a source with no term in a covariate has no contrast in it to give, and nodes that vary it manufacture a null one. See covStrata().

strata_nodes

(Optional) strata per stratified covariate (default 9); a discrete covariate is cut at its levels instead. Each stratum costs a solve and more of them do not buy accuracy: a matched one-covariate fit recovers 0.7000 / 0.7002 / 0.7005 at 3 / 4 / 10 strata against a true 0.700. See covStrata().

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.

Gauss-Hermite nodes per covariate when method = "gh" integrates a study's cov_dist (default 7). Total covariate points are cov_nodes^p for p covariates. Ignored by "mc", which draws a covariate value per simulated subject instead, and by "fo", which cannot integrate a covariate at all.

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", "covariate", "nll", "par"). Defaults to all five. Note "cov" (predicted vs observed covariance) and "covariate" (covariate effect) are different panels.

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

Random seed for reproducibility.

Unused.

An admFit object.

Passed to print.nlmixr2FitCore.