Constructs a control object for est = "admc", the Monte Carlo aggregate
data modelling estimator.
Usage
admControl(
studies = list(),
n_sim = 5000L,
sampling = c("sobol", "halton", "torus", "lhs", "rnorm"),
algorithm = NULL,
maxeval = 500L,
ftol_rel = .Machine$double.eps^2,
print = 10L,
seed = 12345L,
cores = 1L,
grad = c("sens", "fd", "cfd", "none"),
grad_h = 1e-04,
cov_h = 0.001,
cov_h_outer = .Machine$double.eps^(1/5),
grad_bounds = 5,
covMethod = c("r", "none"),
cov_n_sim = 10000L,
n_restarts = 1L,
restart_sd = 0.5,
workers = 1L,
rxControl = NULL,
calcTables = FALSE,
compress = TRUE,
ci = 0.95,
sigdig = 4,
sigdigTable = NULL,
addProp = c("combined2", "combined1"),
optExpression = TRUE,
sumProd = FALSE,
literalFix = TRUE,
returnAdmr = FALSE,
...
)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)
- n_sim
Number of Monte Carlo samples per NLL evaluation.
- 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).
- seed
Random seed for reproducibility.
- cores
Number of OpenMP threads for
rxSolve().- grad
Gradient mode:
"sens"(sensitivity equations, default),"fd"(forward finite differences),"cfd"(central finite differences), or"none"(derivative-free). A warning is issued when"sens"is requested but the sensitivity model is unavailable; the estimator then falls back to forward finite differences.- grad_h
Step size for finite-difference gradient evaluation during optimization (used by
grad = "fd"or"cfd"). The default 1e-4 is near the optimal balance between truncation error (grows withh) and MC noise amplification (grows as1/h) for forward FD. Central FD ("cfd") has a slightly wider optimum around 1e-3, but 1e-4 works well for both.- 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.- grad_bounds
Box-constraint half-width when using gradients.
- covMethod
Covariance method:
"r"(numerical Hessian for structural and residual-error parameters only; omega/IIV SEs are not computed, consistent with nlmixr2 FOCEI) or"none".- 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> 1use fork workers on Unix/macOS (outside RStudio) and a PSOCK cluster on Windows or inside RStudio. RStudio on Linux:future::supportsMulticore()returnsFALSEinside RStudio even on Linux, so PSOCK is used; fork is only active when running from a terminal. Workers are stopped automatically after the restart phase so all cores are available for the Hessian step.- rxControl
rxode2::rxControl()object. Created automatically whenNULL.- calcTables, compress, ci, sigdig, sigdigTable, optExpression, sumProd, literalFix
Passed to
nlmixr2est::foceiControl()for the table/output machinery.- 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).- ...
Additional arguments (none allowed; triggers an error).
Examples
# Minimal control object -- inspect defaults
ctl <- admControl()
ctl$n_sim
#> [1] 5000
ctl$algorithm
#> [1] "NLOPT_LD_LBFGS"
# Override key settings without fitting
ctl2 <- admControl(
n_sim = 2000L,
maxeval = 300L,
grad = "fd",
seed = 42L
)
# \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); 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 = "admc",
control = admControl(
studies = list(study1 = list(E = E, V = V, n = length(ids),
times = times, ev = et(amt = 100))),
n_sim = 1000L,
maxeval = 200L
)
)
#>
#>
#>
#>
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> === admixr2: Aggregate Data Modeling (MC) ===
#> Studies: 1 | MC samples: 1000 | Params: 11 | Cores: 1 | Grad: Sens | Restarts: 1
#> +----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+
#> | | -2LL | tcl | tv1 | tv2 | tq | tka | prop.sd | eta.cl | eta.v1 | eta.v2 | eta.q | eta.ka |
#> +----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+----------+
#> | 0010 | -3644.83 | 5.145 | 11.61 | 26.16 | 10.87 | 1.243 | 0.2005 | 0.09031 | 0.09022 | 0.09008 | 0.09076 | 0.09031 |
#> | 0020 | -3690.08 | 4.991 | 10.75 | 29.19 | 9.79 | 1.082 | 0.1991 | 0.1067 | 0.1036 | 0.09614 | 0.105 | 0.1065 |
#> | 0030 | -3690.65 | 4.961 | 10.25 | 29.9 | 9.857 | 1.04 | 0.1985 | 0.1045 | 0.109 | 0.1022 | 0.1072 | 0.09699 |
#> | 0040 | -3690.72 | 4.95 | 10.16 | 30.02 | 9.814 | 1.031 | 0.1983 | 0.1044 | 0.1139 | 0.1022 | 0.1067 | 0.09433 |
#> | 0045 ✓ | -3690.73 | 4.952 | 10.12 | 30.03 | 9.816 | 1.026 | 0.1983 | 0.1045 | 0.1157 | 0.1008 | 0.1064 | 0.09239 |
#> | 8.3 sec | | | | | | | | | | | | |
#> Computing covariance (R method, Sens-Hessian, 7 gradient evaluations)
#> Note: covMethod='r' computes covariance for structural and sigma parameters only; omega (IIV) SEs are not computed (matching nlmixr2 FOCEI behavior).
#> → compress origData in nlmixr2 object, save 1160
print(fit)
#> ── nlmixr² admc ──
#>
#> OBJF AIC BIC Log-likelihood
#> admc -3690.729 -3668.729 -3598.199 1845.365
#>
#> ── Time (sec fit$time): ──
#>
#> optimize covariance elapsed
#> 1 8.299 8.851 17.15
#>
#> ── Population Parameters (fit$parFixed or fit$parFixedDf): ──
#>
#> Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tcl 1.6 0.01658 1.036 4.952 (4.794, 5.115) 33.2
#> tv1 2.314 0.08726 3.77 10.12 (8.529, 12.01) 35.0
#> tv2 3.402 0.04059 1.193 30.03 (27.74, 32.52) 32.6
#> tq 2.284 0.02142 0.938 9.816 (9.413, 10.24) 33.5
#> tka 0.02614 0.08192 313.4 1.026 (0.8742, 1.205) 31.1
#> prop.sd 0.1983 0.1983
#>
#> 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.
# }