Control fitting and fitted-object storage

drm_control() collects optimizer settings and storage choices for drmTMB(). Select nlminb() (the default) or a stats::optim() method with optimizer_method. Use optimizer for that method's control settings, or optimizer_preset for named nlminb() budgets that keep ordinary defaults fast while making complex refits easier to write. Use the storage flags when a fitted object should keep less R-side state, for example during large-data experiments where the original data frame and TMB automatic-differentiation object are expensive to retain.

Usage

drm_control(
  optimizer = list(),
  se = TRUE,
  se_report_covariance = TRUE,
  se_skip_delta_method = FALSE,
  se_group_sd = FALSE,
  keep_data = TRUE,
  keep_model_frame = TRUE,
  keep_tmb_object = TRUE,
  sparse_fixed = FALSE,
  aggregate_gaussian = FALSE,
  logsigma_clamp = c(-12, 12),
  logsigma_clamp_margin = 3,
  optimizer_preset = c("default", "careful", "robust"),
  newton_polish = TRUE,
  multi_start = 1L,
  fallback_optimizer = NULL,
  start = NULL,
  parallel = c("none", "multicore", "cluster"),
  workers = NULL,
  optimizer_method = "nlminb"
)

Arguments

optimizer: Named list passed to the control argument of stats::nlminb() or stats::optim(), according to optimizer_method. Direct optim() methods use maxit = 1000L unless overridden, with their native tolerance defaults. "L-BFGS-B" uses factr and pgtol; the other methods use reltol and abstol. Unknown controls and controls for other methods are rejected. fnscale must be finite and positive: fitting minimizes the negative log-likelihood. lower and upper are not control settings and are rejected; model parameters are transformed internally, and no user-supplied box bounds are provided.
se: Logical; compute standard errors and fixed-effect covariance with TMB::sdreport() after optimization. Set to FALSE to keep fitted coefficients, fitted values, residuals, predictions, simulations, and profile-likelihood paths while skipping Wald standard errors, stats::vcov(), and Wald confidence intervals. Profile-likelihood intervals still require keep_tmb_object = TRUE.
se_report_covariance: Logical; passed to the getReportCovariance argument of TMB::sdreport(). The default TRUE builds the full covariance matrix of every ADREPORTed quantity. Models with a direct-SD surface (sd_phylo(...) ~ .) report one standard deviation per group, so this matrix is n_group x n_group and its memory cost grows with the square of the number of groups; at ten thousand tips it dominates the fit. Set to FALSE to keep per-quantity standard errors while skipping their joint covariance.
se_skip_delta_method: Logical; passed to the skip.delta.method argument of TMB::sdreport(). Set to TRUE to skip standard errors for ADREPORTed quantities entirely while retaining fixed-effect standard errors and stats::vcov(). This is the cheapest route to Wald inference on the fixed effects for large structured models.
se_group_sd: Logical; report delta-method standard errors for the per-group direct-SD surface (sd_phylo(...) ~ .). Defaults to FALSE. The surface has one standard deviation per group, so ADREPORTing it makes the joint ADREPORT covariance n_group x n_group. Under REML = TRUE the fixed effects are integrated into the Laplace random block and stats::vcov() reads exactly that joint covariance, so at ten thousand tips a bivariate fit needs tens of gigabytes for it. The fitted per-group standard deviations themselves are always available (they are recomputed from the parameters); only their standard errors are opt-in. Set to TRUE to restore the pre-0.3.0 behaviour.
keep_data: Logical; keep the complete-case model data in the fitted object. Set to FALSE to drop fit$data and fit$model$data after fitting. Prediction, fitted values, residuals, simulation, and basic summaries still use the stored model matrices and response vectors.
keep_model_frame: Logical; keep model frames in the fitted object. Set to FALSE to drop fit$model$model_frame and random-effect scale model frames after fitting. Prediction, fitted values, residuals, simulation, sigma(), rho12(), corpairs(), and check_drm() use stored model matrices, terms, response vectors, offsets, and response-name metadata.
keep_tmb_object: Logical; keep the TMB automatic-differentiation object in fit$obj. Set to FALSE to reduce fitted-object size after optimization. check_drm() will then report the fixed-gradient check as a note because it cannot re-evaluate the gradient without fit$obj, and profile-likelihood confidence intervals will be unavailable.
sparse_fixed: Logical; opt-in control for sparse fixed-effect design matrices. This option is limited to univariate Gaussian mu fixed effects with no random effects and intercept-only sigma.
aggregate_gaussian: Logical; opt-in control for sufficient- statistic row aggregation in univariate Gaussian fixed-effect models. This option rejects random effects, structured effects, known sampling covariance, bivariate models, non-Gaussian families, non-unit likelihood weights, and combined sparse fixed-effect matrices.
logsigma_clamp: Numeric c(lo, hi) band, or NULL. Bounds the per-observation Gaussian log(sigma) with an identity-in-band soft-clamp, a numerical guard against scale overflow on near-degenerate per-group scale models. The default c(-12, 12) is identity for any standardized response; widen it for legitimately huge-variance unstandardized data, or set NULL to disable the guard entirely. It is a numerical guard only and does not change identifiability.
logsigma_clamp_margin: Positive number; the soft-clamp saturation margin beyond logsigma_clamp (default 3, saturating to [-15, 15] for the default band). Ignored when logsigma_clamp = NULL.
optimizer_preset: Optimizer-budget preset. "default" adds no optimizer controls, "careful" sets iter.max = 1000 and eval.max = 1000, and "robust" sets iter.max = 5000 and eval.max = 5000. The optimizer escalates this ladder automatically when a preset does not converge. Presets apply only to optimizer_method = "nlminb".
newton_polish: Logical (default TRUE). stats::nlminb()'s own rel.tol/x.tol stopping rules trigger on a relative change in the objective or the step, not on the gradient norm, so on a flat or weakly-curved surface (for example a variance component identified by few groups) nlminb() can report convergence while the true gradient at that point is still far from zero. Rather than tightening rel.tol/x.tol – which in the tested flat-surface case instead produced a spurious "singular convergence" code without reliably landing closer to the true optimum – the default fixed rel.tol/x.tol (nlminb()'s own) are left alone, and a Newton polish step is applied to whatever the selected optimizer returns: the exact TMB gradient plus a finite-differenced Hessian of that gradient (stats::optimHess()) drive a few Newton steps until the gradient's max absolute component is at or below a small threshold, accepting a step only if it does not increase the objective. Cost is a handful of extra gradient evaluations at the reported optimum (skipped entirely once the gradient is already small); set to FALSE to restore the pre-0.7.1 behaviour of reporting nlminb()'s raw optimum unpolished. The same polish applies to direct optim() methods; set newton_polish = FALSE to inspect their raw results.
multi_start: Whole number >= 1 (default 1). With multi_start > 1, each nlminb() preset or the directly selected optim() method is run from multi_start starting points – the principled start plus reproducibly perturbed starts – and the lowest- objective result is kept. Opt-in robustness for weakly identified models; multi_start = 1 is the single-start fit and is unchanged. Joint models with both scale and zero-inflation predictors (for example sigma ~ ... together with zi ~ ... on Poisson or nbinom2() families) can admit multiple local optima; when a fit looks unstable or convergence checks flag a boundary, try multi_start > 1 before changing the model. The starts run one after another unless parallel is set.
fallback_optimizer: NULL (default) or one stats::optim() method ("BFGS", "L-BFGS-B", "Nelder-Mead", "CG"). When set, and no nlminb() preset converges, drmTMB() tries this optimizer as a final attempt; a different algorithm can succeed on a numerically awkward but identified problem. Opt-in; NULL keeps the nlminb()-only ladder. Available only with optimizer_method = "nlminb".
start: NULL (default) or a named list of starting values, keyed by public start labels. Three label families are supported: "fixef:<dpar>:<column>" for a fixed-effect coefficient (already on the model's internal/link scale, e.g. log(sigma) for the default sigma link); "sd:<dpar>:<term>" for a random-effect standard deviation, given on the natural (positive) scale and transformed internally via log(); and "cor:<dpar>:<term>" for a random-effect correlation, given on the natural (-1, 1) scale and transformed internally via atanh(). Labels are validated against the parsed formula and family; unknown labels error before optimization starts. A start only changes where the optimizer begins, never what is reported: unnamed components keep the ordinary family-builder default. Latent random-effect values (u) are not addressable through this contract.
parallel, workers: Backend and worker count for the multi_start starts only: "none" (default) runs them one after another, "multicore" forks Unix workers, and "cluster" uses socket (PSOCK) workers on any platform. A single fit, and multi_start = 1, always run on one core. Parallel starts each begin from the same initial TMB state, so with random effects they can differ from the serial run at the inner solver's tolerance; fixed-effect fits are identical. workers = NULL uses about half the detected cores, capped at multi_start and at 10. See drmTMB-parallel.
optimizer_method: One exact method name: "nlminb" (default), "BFGS", "L-BFGS-B", "Nelder-Mead", or "CG". A direct optim() method runs immediately, without the nlminb() preset ladder or fallback. These methods require optimizer_preset = "default" and fallback_optimizer = NULL. The selected method is recorded in fit$optimizer_used$optimizer and fit$optimizer_attempts$optimizer; its attempt has optimizer_preset = NA because no preset applies. optim() evaluation counts are stored in fit$opt$evaluations; it does not report an iteration count.

Value

A drm_control object.

Details

For optimizer-only settings, control = list(eval.max = 1000) remains valid. When using drm_control(), put optimizer arguments inside optimizer = list(...); do not pass eval.max directly to drm_control(). Under engine = "julia" the forwarded set is a closed whitelist: optimizer$g_tol and optimizer$algorithm on the base bridge, plus optimizer$q4_vcov on the bivariate q = 4 phylogenetic route. On that q4 route optimizer$g_tol becomes DRM.jl's outer-gradient tolerance (q4_g_tol) and optimizer$algorithm is refused, because that route's optimiser has no solver-selection setting. Every other drm_control() field is rejected with an error naming it – never silently ignored – because it describes an nlminb()/TMB program that DRM.jl does not run; use engine = "tmb" for those. The rejected set is derived from drm_control() itself, so a control added here is rejected on the Julia path without a matching edit to the bridge. The structured, bivariate q2 structured, and cross-family Julia routes are narrower still: they accept only a default control, so even optimizer$g_tol is refused there. Presets "careful" and "robust" expand to explicit iter.max and eval.max controls for nlminb(). Values in optimizer override values from the selected preset.

When the selected preset uses the standard nlminb() budget and no explicit optimizer controls, drmTMB() escalates the preset ladder ("default" -> "careful" -> "robust") when an attempt either raises an error (such as a non-finite gradient) or does not converge cleanly (a nonzero convergence code or a non-finite objective). The first cleanly-converged attempt is returned; if no preset converges, the best (lowest-objective) attempt is returned and the fit-time convergence warning flags it. Every attempted preset is recorded in fit$optimizer_attempts and the selected one in fit$optimizer_used.

Examples

dat <- data.frame(y = rnorm(20), x = rnorm(20))
fit <- drmTMB(
  bf(y ~ x, sigma ~ 1),
  data = dat,
  control = drm_control(
    optimizer_preset = "careful",
    se = FALSE,
    keep_data = FALSE,
    keep_model_frame = FALSE,
    keep_tmb_object = FALSE
  )
)