Constructs a control object for est = "adirmc", the Iterative Reweighting
Monte Carlo estimator.
Usage
adirmcControl(
studies = list(),
n_sim = 2500L,
outer_iter = 50L,
sampling = c("sobol", "halton", "torus", "lhs", "rnorm"),
algorithm = NULL,
maxeval = 5000L,
ftol_rel = .Machine$double.eps,
print = 1L,
omega_expansion = 1,
seed = 12345L,
cores = rxode2::rxCores(),
nDisplayProgress = .Machine$integer.max,
grad = c("analytical", "none", "fd"),
kappa_method = c("exact", "linearized", "linearized_gh"),
kappa_n_nodes = 5L,
grad_h = 1e-06,
cov_h = 0.001,
cov_h_outer = .Machine$double.eps^(1/5),
phases = c(2, 1, 0.5, 0.01),
convcrit = 1e-05,
max_worse = 5L,
covMethod = c("r", "none"),
cov_n_sim = 10000L,
n_restarts = 1L,
restart_sd = 0.2,
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,
...
)Arguments
- studies
Named list of study specifications. Each element is a list with:
E– observed mean vectorV– observed covariance matrix or variance vector (auto-detected)n– sample sizetimes– numeric vector of observation timesev–rxode2::et()dosing event tablemethod–"cov"or"var"(optional; auto-detected fromV)
Multi-compartment (multiple observed outputs). To fit several observed compartments simultaneously (e.g. plasma and brain/CSF), give the study an
observationslist instead of top-levelE/V/times. Each entry is one observed output with its ownoutput(the model prediction variable, e.g."cp"or"cCSF"),times,E,Vand – for independent fits –evandn. Pass the endpoint names toadmData(), 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-2LLis their sum.Joint (same subjects) – the outputs are measured on the SAME subjects. Give the study a shared
nandev, and a joint covariance either as a study-level full matrixV(blocks inobservationsorder) or as per-output marginalVplus acrosslist of cross-covariance blocks keyed"outA:outB"(eachlength(times_A)xlength(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
observationslist, a study may carry adataframe that keys each observed summary by endpoint, the way nlmixr2 keys observations byDVID/CMT. The frame needs an endpoint column (DVID,CMToroutput), a time column (TIME), a mean column (E) and – unless a jointVis given – a variance column (V) or an SD column (SD). It is normalised into exactly the same units as theobservationsform, 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 explicitjoint = 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 indata.- n_sim
Number of Monte Carlo samples per NLL evaluation.
- outer_iter
Maximum inner optimiser iterations per phase.
- sampling
Sampling method for eta draws:
"sobol"(Sobol, default),"halton"(Halton),"torus"(Kronecker/torus),"lhs"(Latin hypercube), or"rnorm"(iid normal).- algorithm
nloptr algorithm string, or
NULL(default) to pick the default that matchesgrad:"NLOPT_LD_LBFGS"with a gradient,"NLOPT_LN_BOBYQA"whengrad = "none". Any algorithm reported bynloptr::nloptr.print.options()is accepted (e.g."NLOPT_LD_MMA","NLOPT_LN_NELDERMEAD"). An explicit algorithm is reconciled withgrad: whengrad = "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 number of optimizer function evaluations.
- ftol_rel
Relative function-value tolerance for convergence.
Print progress every this many evaluations (0 = silent).
- omega_expansion
Inflate proposal Omega by this factor (>= 1).
- seed
Random seed for reproducibility.
- cores
Number of OpenMP threads for
rxSolve(). Defaults torxode2::rxCores().rxSolve()parallelises over subjects, so this is the main speed lever for the MC estimators; whenworkers > 1it is a total budget, split across the workers.- nDisplayProgress
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.- grad
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_method
Kappa correction method for models with non-mu-referenced struct thetas:
"exact"(default, re-evaluates population predictionf(theta, 0)via rxSolve at each inner step),"linearized"(precomputesJ = df/d(theta)once per outer iteration usingf(theta, 0)as baseline — zero rxSolve per inner step), or"linearized_gh"(same linear approximation but baseline and Jacobian use Gauss-Hermite quadratureE_GH[f(theta, eta)]instead off(theta, 0)— more accurate baseline at any IIV magnitude, still zero rxSolve per inner step).- kappa_n_nodes
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.- grad_h
Step size for the inner optimiser's finite-difference gradient (
grad = "fd"). Defaults to1e-6, not the1e-4the 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 sets1e-4elsewhere does not apply. The inner step was a hard-coded1e-6until it was made to honourgrad_h; inheriting the coarser default would have changed the gradient the loop was tuned for.- cov_h
Inner FD step for the gradient-based Hessian (only used when
covMethod = "r"andgrad != "none"). Each gradient evaluation has MC noise of ordersigma / cov_h; the Hessian divides that noise by the outer step, giving total noisesigma / (cov_h * cov_h_outer * |p|).cov_h = 1e-3balances truncation error and noise amplification. Increase to1e-2if the Hessian is non-positive definite.- cov_h_outer
Outer step scale for the numerical Hessian. The actual step for parameter
pismax(|p|, 0.1) * cov_h_outer. Applied to both the gradient-FD Hessian (grad != "none") and the NLL-FD Hessian (grad = "none"). Defaulteps^(1/5)(~2.5e-3) is larger than the textbookeps^(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. to5e-3or1e-2) if the Hessian is non-positive definite.- phases
Numeric vector of box-constraint half-widths, one per phase. Phases progressively tighten the search region.
- convcrit
Convergence criterion: phase ends when
|approx - exact| < convcrit.- max_worse
Stop a phase after this many consecutive worsening iterations.
- covMethod
Covariance method:
"r"(numerical Hessian over the structural, residual-error and omega parameters) or"none". 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.All three blocks are reported on the scale the ESTIMATES are printed on, as
nlmixr2estdoes: structural thetas on the log/optimizer scale, residual error as an SD, and omega as the variance/covariance entries (namedom.<eta>andcov.<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.- cov_n_sim
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 andcov_n_sim = 10000(default) is sufficient for most models.- n_restarts
Number of optimization restarts. Runs in parallel when
workers > 1.- restart_sd
Standard deviation of structural theta perturbations for restart initialisation.
- workers
Number of parallel workers for multi-restart.
1(default) runs restarts sequentially. Values> 1run the restarts on a pool of background R processes (mirai daemons), which behaves the same way on every platform. Requires themiraipackage. 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.- rxControl
rxode2::rxControl()object. Created automatically whenNULL.- 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 torxode2::rxSolve()'s ownsigdigargument for every solve the estimator issues – rxode2 owns the mapping toatol/rtoland has changed it between releases, which is why the digits, not the tolerances, are what travels – and tonlmixr2est::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) andcov_h_outer(~2.5e-3), whilesigdig = 4maps 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 (everySEreportedNA) 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 atsigdig = 4with standard errors unchanged to 4 significant figures. Elsewhere, compare the objective and the standard errors againstNULLbefore relying on it. Table formatting is unaffected either way:sigdigTabledefaults to 4 regardless.- addProp
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 tonlmixr2est::foceiControl()for the table/output machinery only.- returnAdmr
If
TRUE, return a plain list instead of a full nlmixr2 fit object (useful for debugging).- resid_nodes
Gauss-Hermite nodes used to integrate the RESIDUAL for a transform-both-sides endpoint (
boxCox,yeoJohnson,logitNorm,probitNorm), wherey = 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_nodesin isolation (~50 us at 15, 300 us at 81 for an 8-row study) but negligible beside the ODE solve: a full NLL evaluation measured 0.750 s per 60 evaluations at BOTH 31 and 81 nodes. Raise it if you have a saturating endpoint with a large residual SD; there is little to gain by lowering it.- ...
Additional arguments (none allowed; triggers an error).
Details
Multi-compartment fits (a study observations list with several observed
outputs) are not supported by adirmc; use est = "admc", "adfo", or
"adgh" for those. Single-output studies are fit as usual.
Examples
# Inspect defaults
ctl <- adirmcControl()
ctl$phases
#> [1] 2.00 1.00 0.50 0.01
ctl$omega_expansion
#> [1] 1
# Tighter phases, more restarts
ctl2 <- adirmcControl(
n_sim = 1000L,
omega_expansion = 1.5,
phases = c(2, 1, 0.5, 0.01),
n_restarts = 3L
)
# \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 <- diag(diag(cov.wt(dv_mat, method = "ML")$cov))
pk_model <- function() {
ini({
tcl <- log(5); tv1 <- log(12); tv2 <- log(25)
tq <- log(12); tka <- log(1.2)
prop.sd <- c(0, 0.2)
eta.cl ~ 0.09; eta.v1 ~ 0.09; eta.v2 ~ 0.09
eta.q ~ 0.09; eta.ka ~ 0.09
})
model({
cl <- exp(tcl + eta.cl); v1 <- exp(tv1 + eta.v1)
v2 <- exp(tv2 + eta.v2); q <- exp(tq + eta.q)
ka <- exp(tka + eta.ka)
d/dt(depot) <- -ka * depot
d/dt(central) <- ka * depot - (cl/v1 + q/v1) * central + (q/v2) * peripheral
d/dt(peripheral) <- (q/v1) * central - (q/v2) * peripheral
cp <- central / v1
cp ~ prop(prop.sd)
})
}
fit <- nlmixr2(
pk_model, admData(), est = "adirmc",
control = adirmcControl(
studies = list(study1 = list(E = E, V = V, n = length(ids),
times = times, ev = et(amt = 100))),
n_sim = 500L
)
)
#>
#>
#>
#>
#> ℹ 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...
#> → optimizing duplicate expressions in admixr2 sensitivity model...
#>
#>
#>
#>
#>
#>
#> === admixr2: Aggregate Data Modeling (IR-MC) ===
#> Studies: 1 | MC samples: 500 | Phases: 4 | Iters/phase: 50 | Expansion: 1.00 | Grad: analytic+Sens-Hessian | Restarts: 1
#> +----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+
#> | | -2LL | tcl | tv1 | tv2 | tq | tka | prop.sd | eta.cl | eta.v1 | eta.v2 | eta.q | eta.ka |
#> +-- Phase 1: Wide (+/-2.00) -------------------------------------------------------------------------------------------------------------------+
#> | 0001 | 728.22 | 3.552 | 6.49 | 36.82 | 8.796 | 0.5216 | 0.1867 | 0.4043 | 0.1828 | 0.665 | 0.02393 | 0.2747 |
#> | 0002 | -1247.61 | 4.966 | 6.662 | 33.45 | 7.927 | 0.6548 | 0.1217 | 0.1298 | 0.02474 | 0.1283 | 0.02614 | 0.2091 |
#> | 0003 | -1194.49 | 5 | 8.026 | 33.06 | 8.569 | 0.7865 | 0.1911 | 0.09852 | 0.1406 | 0.07392 | 0.006252 | 0.1188 |
#> | 0004 | -1238.75 | 4.872 | 7.888 | 31.64 | 8.864 | 0.7894 | 0.1762 | 0.1499 | 0.07673 | 0.1169 | 0.006838 | 0.1372 |
#> | 0005 | -1264.67 | 4.955 | 7.829 | 31.56 | 8.832 | 0.79 | 0.1859 | 0.1141 | 0.07267 | 0.08511 | 0.006873 | 0.1292 |
#> | 0006 | -1266.43 | 4.959 | 8.001 | 31.16 | 8.888 | 0.8065 | 0.1863 | 0.119 | 0.07429 | 0.08569 | 0.006851 | 0.1323 |
#> | 0007 | -1266.44 | 4.959 | 8.004 | 31.15 | 8.885 | 0.8061 | 0.1863 | 0.119 | 0.07429 | 0.0857 | 0.006851 | 0.1322 |
#> | 0008 ✓ | -1266.44 | 4.959 | 8.002 | 31.15 | 8.887 | 0.8061 | 0.1863 | 0.119 | 0.07429 | 0.08571 | 0.00685 | 0.1322 |
#> +-- Phase 2: Focused (+/-1.00) ----------------------------------------------------------------------------------------------------------------+
#> | 0009 | -1266.44 | 4.959 | 8.002 | 31.15 | 8.887 | 0.8061 | 0.1863 | 0.119 | 0.07429 | 0.08571 | 0.00685 | 0.1322 |
#> | 0010 | -1266.44 | 4.959 | 7.991 | 31.15 | 8.886 | 0.8054 | 0.1862 | 0.119 | 0.07433 | 0.08583 | 0.006848 | 0.1319 |
#> | 0011 | -1266.44 | 4.959 | 7.991 | 31.16 | 8.886 | 0.8054 | 0.1862 | 0.119 | 0.07433 | 0.08583 | 0.006848 | 0.1319 |
#> | 0012 | -1266.44 | 4.956 | 7.984 | 31.2 | 8.882 | 0.8044 | 0.1863 | 0.1191 | 0.07447 | 0.08599 | 0.006845 | 0.132 |
#> | 0013 | -1266.44 | 4.957 | 7.983 | 31.2 | 8.88 | 0.8043 | 0.1862 | 0.1191 | 0.07448 | 0.086 | 0.006845 | 0.1319 |
#> | 0014 | -1266.44 | 4.956 | 7.982 | 31.2 | 8.881 | 0.8043 | 0.1862 | 0.1191 | 0.07448 | 0.086 | 0.006845 | 0.1319 |
#> | 0015 | -1266.45 | 4.957 | 7.974 | 31.19 | 8.881 | 0.8038 | 0.1861 | 0.1191 | 0.07453 | 0.0861 | 0.006843 | 0.1317 |
#> | 0016 | -1266.45 | 4.957 | 7.974 | 31.19 | 8.881 | 0.8038 | 0.1861 | 0.1191 | 0.07453 | 0.0861 | 0.006843 | 0.1317 |
#> | 0017 | -1266.22 | 4.923 | 7.869 | 31.59 | 8.843 | 0.7939 | 0.185 | 0.1201 | 0.08131 | 0.09585 | 0.006669 | 0.1294 |
#> | 0018 | -1266.41 | 4.939 | 7.867 | 31.58 | 8.843 | 0.7933 | 0.1854 | 0.1188 | 0.07741 | 0.09027 | 0.006751 | 0.1287 |
#> | 0019 | -1266.47 | 4.937 | 7.862 | 31.54 | 8.859 | 0.7941 | 0.1856 | 0.119 | 0.07746 | 0.09029 | 0.00675 | 0.1289 |
#> | 0020 | -1266.48 | 4.939 | 7.872 | 31.51 | 8.851 | 0.7944 | 0.1858 | 0.1191 | 0.07757 | 0.09032 | 0.00675 | 0.1293 |
#> | 0021 | -1266.48 | 4.938 | 7.874 | 31.52 | 8.848 | 0.7945 | 0.1858 | 0.1191 | 0.07757 | 0.09032 | 0.00675 | 0.1293 |
#> | 0022 ✓ | -1266.48 | 4.938 | 7.874 | 31.52 | 8.849 | 0.7944 | 0.1858 | 0.1191 | 0.07758 | 0.09032 | 0.00675 | 0.1293 |
#> +-- Phase 3: Fine-tuning (+/-0.50) ------------------------------------------------------------------------------------------------------------+
#> | 0023 | -1266.48 | 4.938 | 7.874 | 31.52 | 8.85 | 0.7945 | 0.1858 | 0.1191 | 0.07758 | 0.09032 | 0.00675 | 0.1293 |
#> | 0024 ✓ | -1266.48 | 4.938 | 7.874 | 31.52 | 8.85 | 0.7945 | 0.1858 | 0.1191 | 0.07758 | 0.09032 | 0.00675 | 0.1293 |
#> +-- Phase 4: Precision (+/-0.01) --------------------------------------------------------------------------------------------------------------+
#> | 0025 | -1266.48 | 4.938 | 7.883 | 31.51 | 8.85 | 0.7952 | 0.1856 | 0.119 | 0.07768 | 0.09035 | 0.006749 | 0.1295 |
#> | 0026 | -1266.48 | 4.938 | 7.883 | 31.5 | 8.851 | 0.7953 | 0.1856 | 0.1191 | 0.07768 | 0.09036 | 0.006749 | 0.1295 |
#> | 0027 | -1266.48 | 4.939 | 7.884 | 31.5 | 8.85 | 0.7953 | 0.1856 | 0.1191 | 0.07769 | 0.09036 | 0.006749 | 0.1295 |
#> | 0028 | -1266.48 | 4.939 | 7.884 | 31.5 | 8.851 | 0.7953 | 0.1856 | 0.1191 | 0.07769 | 0.09036 | 0.006749 | 0.1295 |
#> | 0029 ✓ | -1266.48 | 4.939 | 7.885 | 31.5 | 8.851 | 0.7954 | 0.1856 | 0.1191 | 0.07769 | 0.09037 | 0.006749 | 0.1295 |
#> | 1.4 sec | | | | | | | | | | | | |
#> Computing covariance (R method, MC NLL, Sens-Hessian, 12 gradient evaluations)
#> → compress origData in nlmixr2 object, save 1160
#>
#>
print(fit)
#> ── nlmixr² adirmc ──
#>
#> OBJF AIC BIC Log-likelihood
#> adirmc -1266.481 -1244.481 -1173.951 633.2406
#>
#> ── Time (sec fit$time): ──
#>
#> optimize covariance other elapsed
#> 1 1.424 13.914 0 15.338
#>
#> ── Population Parameters (fit$parFixed or fit$parFixedDf): ──
#>
#> Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tcl 1.597 0.04576 2.865 4.939 (4.515, 5.402) 35.57 NaN
#> tv1 2.065 0.3649 17.67 7.885 (3.857, 16.12) 28.42 NaN
#> tv2 3.450 0.1515 4.392 31.50 (23.41, 42.39) 30.75 NaN
#> tq 2.180 0.1191 5.463 8.851 (7.008, 11.18) 8.229 NaN
#> tka -0.2290 0.3296 143.9 0.7954 (0.4169, 1.517) 37.18 NaN
#> prop.sd 0.1856 0.01563 8.418 0.1856 (0.1550, 0.2163)
#>
#> Covariance Type (fit$covMethod): r
#> No correlations in between subject variability (BSV) matrix
#> Full BSV covariance (fit$omega) or correlation (fit$omegaR; diagonals=SDs)
#> Distribution stats (mean/skewness/kurtosis/p-value) available in fit$shrink
#> Censoring (fit$censInformation): No censoring
#> Minimization message (fit$message):
#> NLOPT_XTOL_REACHED: Optimization stopped because xtol_rel or xtol_abs (above) was reached.
# }