Control settings for the Gauss-Hermite (GH) quadrature estimator

Creates a control object for nlmixr2(est = "adgh"). The GH estimator integrates model predictions against the random-effects prior $\eta \sim N(0, \Omega)$ using a deterministic tensor-product Gauss-Hermite quadrature grid. It is unbiased at any IIV magnitude (unlike FO), noise-free (unlike MC), and much faster than MC for models with up to ~4 etas.

Usage

adghControl(
  studies = list(),
  n_nodes = 5L,
  grad = c("analytical", "fd", "none"),
  algorithm = NULL,
  maxeval = 500L,
  ftol_rel = .Machine$double.eps^(1/2),
  print = 10L,
  seed = 12345L,
  cores = rxode2::rxCores(),
  nDisplayProgress = .Machine$integer.max,
  grad_h = 1e-04,
  grad_bounds = 5,
  cov_h = 0.001,
  cov_h_outer = .Machine$double.eps^(1/4),
  covMethod = c("r,s", "r", "none"),
  n_restarts = 1L,
  restart_sd = 0.5,
  workers = 1L,
  rxControl = NULL,
  calcTables = FALSE,
  compress = TRUE,
  ci = 0.95,
  sigdig = NULL,
  sigdigTable = NULL,
  addProp = c("combined2", "combined1"),
  optExpression = TRUE,
  sumProd = FALSE,
  literalFix = TRUE,
  returnAdmr = FALSE,
  resid_nodes = 81L,
  cov_nodes = 7L,
  cov_integration = c("on", "sparse", "off"),
  cov_sparse_level = 3L,
  xtol_rel = .Machine$double.eps^(1/2),
  ...
)

Arguments

studies

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

n_nodes

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.

grad

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.

algorithm

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.

maxeval

Maximum function evaluations (default 500).

ftol_rel

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

print

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

seed

Random seed (used for restarts).

cores

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

nDisplayProgress

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.

grad_h

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

grad_bounds

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.

cov_h

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

cov_h_outer

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

covMethod

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

n_restarts

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

restart_sd

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

workers

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

rxControl

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

calcTables, compress, ci, sigdigTable, optExpression, sumProd, literalFix

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

sigdig

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.

addProp

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

returnAdmr

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

resid_nodes

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.

cov_nodes

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.

cov_integration

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:

rule	p = 3, rho = 0.85	p = 4, rho = 0.85
sparse, level 2	6 pts, 7.8e-04 / 4.7e-02	9 pts, 9.7e-04 / 6.1e-02
product, 3 nodes	27 pts, 4.0e-05 / 1.5e-02	81 pts, 1.5e-04 / 3.2e-02
sparse, level 3	31 pts, 1.6e-06 / 5.0e-04	49 pts, 3.5e-06 / 9.5e-04

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.

cov_sparse_level

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.

xtol_rel

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

...

Unused arguments (trigger an error).

Value

An adghControl object (a named list).

Examples

ctl <- adghControl()
ctl$n_nodes
#> [1] 5
ctl$grad
#> [1] "analytical"

# More nodes for large IIV, analytical gradient
ctl2 <- adghControl(n_nodes = 7L, grad = "analytical", maxeval = 300L)

# \donttest{
library(rxode2)
library(nlmixr2)

data("examplomycin")
obs    <- examplomycin[examplomycin$EVID == 0, ]
obs    <- obs[order(obs$ID, obs$TIME), ]
times  <- sort(unique(obs$TIME))
ids    <- unique(obs$ID)
dv_mat <- do.call(rbind, lapply(ids, function(i) {
  sub <- obs[obs$ID == i, ]; sub$DV[order(sub$TIME)]
}))
E <- colMeans(dv_mat)
V <- cov.wt(dv_mat, method = "ML")$cov

pk_model <- function() {
  ini({
    tcl <- log(5); tv <- log(30)
    prop.sd <- c(0, 0.2)
    eta.cl ~ 0.09; eta.v ~ 0.04
  })
  model({
    cl <- exp(tcl + eta.cl)
    v  <- exp(tv  + eta.v)
    d/dt(central) <- -(cl/v) * central
    cp <- central / v
    cp ~ prop(prop.sd)
  })
}

fit <- nlmixr2(
  pk_model, admData(), est = "adgh",
  control = adghControl(
    studies = list(study1 = list(E = E, V = V, n = length(ids),
                                 times = times, ev = et(amt = 100)))
  )
)
#>  
#>  
#>  
#>  
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> → calculate sensitivities
#> → finding duplicate expressions in admixr2 sensitivity model...
#>  
#>  
#> === admixr2: Aggregate Data Modeling (GH) ===
#>   Obs units: 1 | Params: 5 | Nodes: 5^2=25 | Cores: 2 | Grad: Analytical | Restarts: 1
#> +----------+----------+----------+----------+----------+----------+----------+
#> |          |     -2LL |      tcl |       tv |  prop.sd |   eta.cl |    eta.v |
#> +----------+----------+----------+----------+----------+----------+----------+
#> | 0010     |  1000.18 |    6.203 |    35.45 |   0.3103 |  0.08888 |  0.05562 |
#> | 0020     |   805.78 |    6.666 |    37.33 |   0.3781 |   0.1041 |  0.05946 |
#> | 0030     |   805.77 |    6.663 |    37.35 |   0.3784 |   0.1035 |  0.05848 |
#> | 0031 ✓   |   805.77 |    6.663 |    37.35 |   0.3784 |   0.1035 |  0.05848 |
#> | 0.6 sec  |          |          |          |          |          |          |
#>   Computing covariance (R method, Analytical-Hessian, sandwich, 6 gradient evaluations)
#> → compress origData in nlmixr2 object, save 1160
#>  
#>  
# }

Control settings for the Gauss-Hermite (GH) quadrature estimator

Usage

Arguments

Value

See also

Examples