Skip to contents

Control parameters for frmtmb fits

Source: R/fit.R
frmtmb_control.Rd

Control parameters for frmtmb fits

Usage

frmtmb_control(
  optimizer = "nlminb",
  optCtrl = list(iter.max = 1000, eval.max = 1000),
  restarts = 1,
  grad_tol = 0.001,
  profile = FALSE,
  sparse_x = FALSE,
  autoscale = NULL,
  check_nlev_1 = c("warning", "ignore", "stop"),
  check_olre = c("warning", "ignore", "stop"),
  check_se = c("warning", "ignore", "stop"),
  check_mi = c("warning", "ignore", "stop"),
  importance_seed = 1L,
  importance_rounds = 5L,
  importance_ess = 0.25,
  mo_search = TRUE,
  verbose = NULL,
  quad_nodes = c(11L, 15L),
  quad_max_gb = 4
)

Arguments

optimizer

"nlminb" (default), "optim" (L-BFGS-B), or a function with signature (par, fn, gr, lower, upper, control) returning a list with elements par, objective, convergence (0 = success), and optionally message - the hook for optimx, nloptr, and friends without frmtmb depending on them. Example:

nlopt <- function(par, fn, gr, lower, upper, control) {
  r <- nloptr::nloptr(par, fn, gr, lb = lower, ub = upper,
                      opts = list(algorithm = "NLOPT_LD_LBFGS",
                                  xtol_rel = 1e-10, maxeval = 2000))
  list(par = r$solution, objective = r$objective,
       convergence = as.integer(r$status < 0), message = r$message)
}
frm(..., control = frmtmb_control(optimizer = nlopt))

optCtrl

Control list passed to the built-in optimizers (stats::nlminb() / stats::optim()).

restarts

Number of times to restart the optimizer from the current optimum while the gradient remains above grad_tol.

grad_tol

How close to its optimum a fit has to stop. Read twice, in this order.

As a GRADIENT, it restarts the optimizer while the maximum absolute gradient is above it, and it is the trip-wire on the convergence warning: a fit whose largest gradient component stays under it says nothing further.

As a LOG LIKELIHOOD, it decides the warning. A fit that trips the trip-wire is measured again, this time by how much log likelihood one exact Newton step from the stopping point would still buy, over the parameters no bound holds in place, and it warns only when that is above grad_tol as well. The second reading is what makes the default mean the same thing on every design: an absolute gradient grows with the sample size, because the optimizer stops on a RELATIVE change in an objective that grows with it, while the remaining log likelihood does not. Measured over 720 correct fits on nine designs at four sample sizes (dev/gradcheck-findings.md), the gradient alone warned on 289 of them and the two readings together on none.

WHAT IT DOES NOT COVER. The first reading gates the second, so a fit whose largest gradient never reaches grad_tol is never looked at, however far short of its optimum it is. A predictor column scaled far below one does exactly that: it keeps the gradient small while the coefficient stalls near zero. Measured, with autoscale = FALSE: at a column spread of 1e-7 a gaussian fit reports convergence, a maximum absolute gradient of 1.6e-5, and a log likelihood 297 units below the standardized fit of the same data, with no warning of any kind. The default autoscale rule is what closes that gap, which is why it engages at a spread of 1e-3 (see autoscale), and diagnose() reports badly scaled columns whether or not the gradient noticed.

Lower it for a stricter fit and a stricter warning; both readings move together. diagnose() reports every number the verdict used.

profile

Experimental: move the primary (beta) coefficients into the inner (Laplace) problem, TMB's profile argument - the analog of glmmTMBControl(profile = TRUE) and glmer(nAGQ = 0). Speeds up models with many fixed effects, and like those it is an approximation: estimates differ slightly from the exact fit. Coefficient covariance comes from the joint precision. Not compatible with REML = TRUE or quadrature = TRUE; profile/uniroot confint() and hypothesis(method = "profile") need a non-profiled fit.

sparse_x

Build the parametric fixed-effect design matrices as sparse Matrix::sparse.model.matrix() objects, the analog of glmmTMB(sparseX =). Worth it when a many-level fixed factor makes the dense design dominate memory; estimates are identical either way. model.matrix() on the fit then returns a sparse matrix.

autoscale

Standardize badly scaled continuous predictors internally (the lme4 >= 1.1.37 feature): fit a copy of the model with each qualifying fixed-effect column centered and scaled (scaled only, in a linear predictor without an intercept), map that optimum back to the original parameterization exactly, and warm-start the ordinary fit there. Reported results are always on the original scale, so every downstream method works unchanged; the cost is a second (cheap) optimization. Columns qualify when they are parametric and take more than two distinct values; intercepts, factor contrasts, smooth bases, and mo()/mi() columns are never touched, and the whole step is a silent no-op when nothing qualifies. Compatible with profile = TRUE. Under prior or bounds, the first stage applies them to the scaled coefficients; the second stage is the fit that is reported. A random slope on a rescaled column ((1 + x | g)) is rescaled with it in us(), diag() and the Student-t blocks, its log standard deviation and its effects mapped back exactly; other covariance structures keep their slope unscaled. The pre-fit's own warnings are not shown. Under the default, a pre-fit that errors or does not converge in its coefficients never makes the reported fit less diagnostic than autoscale = FALSE: the default then reports the plain fit, unless the fit from the pre-fit has the better likelihood and either warns at least as often or is a verified optimum. Under TRUE the pre-fit's error stands and its non-convergence is a warning. TRUE always does this, FALSE never does. The default, NULL, does it when a qualifying column has a standard deviation below 1e-3 and its coefficient is optimized directly (not a mu coefficient under REML = TRUE or profile = TRUE, which the inner solver integrates), or when a column carrying such a random slope has one below 0.05. Without it such a fit can stop short and report convergence: a poisson y ~ 0 + x with x on a 1e-6 scale lost 62.6 log-likelihood units against stats::glm(), and a random slope on a column spread 0.01 lost up to 14.4. Every other fit is the fit autoscale = FALSE gives, bit for bit. A column that is far from zero against its spread, such as x + 1e6, does not need this option: whatever its value, the optimizer steps in coordinates in which each ill-conditioned fixed design has orthogonal columns.

check_nlev_1

What to do about a scalar random-effect term whose grouping factor has a single level: "warning" (default), "ignore", or "stop", following lme4's lmerControl() check vocabulary. Such a term has no variance to estimate - the single level is absorbed by the intercept - and its standard deviation collapses to zero. Structured blocks over several terms per level (ar1(), us(), the spatial covariance structures) are never flagged: one grouping level there is a single realization of a field, which is the normal way to write them.

check_olre

What to do about an observation-level random effect - one level per row - on a gaussian, student or lognormal response: "warning" (default), "ignore", or "stop". Its variance is confounded with the residual standard deviation, so only their sum is identified and the split between them is arbitrary. The check is skipped when sigma is not free to absorb it - a se() response or a constant sigma - which is the random-effects meta-analysis, and for discrete families, where an observation-level term is the usual overdispersion model. It is also skipped for the known-structure blocks gr(cov = ), gr(prec = ) and equalto(): their level covariance is not proportional to the identity, so the block and the residual are separately identified. That is the animal model with one measurement per individual.

check_se

What to do when a converged fit has parameters with no standard error: "warning" (default), "ignore", or "stop". frm() looks at the Hessian at the optimum to find out, so the cost grows with the number of parameters. Without random effects it reads the exact Hessian, about one gradient per coefficient, with no limit. With random effects it builds the finite-difference Hessian, two gradients per outer parameter, which summary() then reuses, when the model has at most 10 outer parameters or those gradients are at most a quarter of the objective and gradient evaluations the optimizer made. Otherwise the check waits for the first standard-error use (summary(), fixef(), vcov(), confint(), a prediction's standard error) and warns there; under "stop" it does not wait, so that frm() itself stops. The rule counts work and does not read the clock, so the same call warns in the same place on any machine. Measured against "ignore": 8 to 23 percent on seven small models and on mixed models with 14 to 23 outer parameters, 19 to 24 percent on a glm with 301 coefficients, and 3 percent or less on mixed models with 42 to 403 outer parameters, where the check waited. "ignore" skips it. A model without random effects then builds the exact Hessian at the first standard-error use. That is faster than the finite differences it replaces but takes more memory: on a bernoulli model with 10,000 rows the peak working set rose by 350 MB with 200 coefficients and by 870 MB with 500. The warning names each parameter and the reason: a bound holds it, the data separate the outcomes, the likelihood is flat along it (a standard deviation at 0, a mo() simplex split by a category no row is in, parameters that enter only through a combination), a step along it raises the log-likelihood, or its Hessian row is not finite. The other parameters keep their standard errors, and a prediction, emmean or hypothesis() that moves along a lost direction gets none, with one warning per call. A fit that did not converge warns about that instead. A covariance parameter at the edge of its parameter space (a group standard deviation at zero, a correlation or mixing parameter at its limit) is reported as lme4 reports a singular fit, with a message, "Boundary (singular) fit", of class frmtmb_boundary_fit; "stop" makes it an error too. When the optimizer stops with a non-zero code there (nlminb's "singular convergence") but the gradient is within grad_tol and every lost parameter is at its edge, the fit gets this message and not the convergence warning, as lme4's checkConv() does. A smooth, gp() or hsgp() term at its own limit loses those standard errors without a message, as mgcv reports a smoothing parameter at infinity.

check_mi

What to do when the Laplace approximation of a read mi() value moves an estimate: "warning" (default), "ignore", or "stop". See bf(), section "Missing values under mi()". "ignore" skips the check.

importance_seed

Seed for the standard normal draws of frm(importance =). The draws are taken from a private random stream, so the fit neither reads nor disturbs the session's random state: the same seed gives the same answer to the last bit, whatever else has drawn random numbers.

importance_rounds

Cap on the number of times frm(importance =) refreezes its proposal. Each round rebuilds the proposal at the current estimate and optimizes again, and the loop stops as soon as either the estimates move less than 1e-3 or the freshly anchored objective is stationary at them to grad_tol. Over 648 measured fits of six families that took one round in 398, two in 237, three in 12 and four in 1, so the default leaves one in hand; an unused round costs nothing, and each round actually taken costs one tape. A fit that uses every round warns.

importance_ess

Effective sample size, as a fraction of the draw count, below which frm(importance =) warns and names the groups. The default 0.25 separates two measured regimes rather than clearing every design. A proposal that is working sits near the top of the range: the probe design behind this correction (60 groups of 8 Bernoulli rows, correlated slope) holds a worst group of 0.95 at 1000 draws, and a gaussian response, where the correction has nothing to correct, holds exactly 1 at any draw count. The designs the warning is FOR sit below that, and draws are what move them: a correlated slope in mu plus a sigma block over 12 groups of 10 rows holds a worst group of 0.14 at 500 draws and 0.35 at 2000.

After the fit, refit each mo() term whose simplex has a weight at 0 with its scale coefficient held at the other sign, and keep the better fit (see "Monotonic effects" in frm()). Each sign of the coefficient has its own local maximum, and on brms_monotonic's ls ~ mo(income) * age the search took the fits at the maximum likelihood from 174 of 200 data sets to 199. It costs up to two refits of the whole model per such term, and so grows with the model: with it the fit took 2.0 times the objective evaluations it took without it on that model, and 3.3 to 3.8 times those of frmtmb 0.68.1 on a model of 5000 rows with a random intercept over 100 groups. A weight at 0 of a simplex that a family declares (frmtmb_family(simplex_pars =), as the baseline of frmtmb.survival's cox() does) is checked the same way, at one objective value per such weight, and the model refitted from off the boundary when the likelihood rises from it. FALSE turns both off. A model without mo() terms or such a simplex is not affected.

verbose

Report fit progress through message(), one terse line per stage with its elapsed seconds, so a slow fit shows where the time went. FALSE (default) is silent and costs nothing. TRUE (or 1) reports parsing, frame assembly, the autoscale pre-fit, tape construction, each optimizer run and restart, sdreport() when se = TRUE, and a closing line with the objective, the maximum absolute gradient, and the number of convergence warnings. The fit itself opens with a line naming the family, the mode (ML or REML, plus profile, quadrature, autoscale, prior), and the optimizer. 2 adds the optimizer's own trace, by setting optCtrl$trace unless you set it yourself. That trace is printed by stats::nlminb() / stats::optim() to standard output, not through message(), and a custom optimizer function receives optCtrl unchanged, so verbose does not reach it.

Use suppressMessages() to silence a verbose fit. The refit loops in frm_bootstrap(), influence(), and frm_allfit() force verbose off, so a verbose fit does not make them print hundreds of lines; refit() and frm_multiple() report normally.

quad_nodes

Gauss-Hermite nodes per level for frm(quadrature = TRUE) over two nested blocks: one count for both levels, or the outer level's count and then the inner level's, from 1 to 50. The default c(11, 15) gives more nodes to the inner level, whose small groups are the least Gaussian. On 12 groups of 6 weighted Bernoulli rows in 24 inner groups, the log-likelihood error was 7e-6 at the default, 2.5e-7 at 21 nodes per level and 6e-10 at 31, against 0.85 for the Laplace approximation. Memory and time grow with the rows times the product of the two counts. A model with one block per row does not read this setting: its rule is adaptive Gauss-Kronrod.

quad_max_gb

Limit, in GB, on the memory that frm(quadrature = TRUE) over two nested blocks may add. The fit predicts that memory before the Laplace warm start, from three small probe tapes and a cost per family measured on 20 families, and stops with a message that gives the prediction and the remedies when it is above the limit. On those families the prediction was 1.22 to 1.61 times the measured peak. A family not measured is predicted with the largest measured cost, so its prediction can be several times its peak. At the default nodes, 3000 Bernoulli rows took about 0.5 GB, 3000 Student rows 5.3 GB, and 40 censored Beta rows 3.7 GB.

Value

A list of control settings.

Examples

set.seed(1)
n <- 200
dd <- data.frame(x = rnorm(n), g = factor(rep(1:10, 20)))
dd$y <- rnorm(n, 1 + 0.5 * dd$x + rnorm(10, 0, 0.8)[dd$g], 1)

# another optimizer, with its own control list
fit <- frm(bf(y ~ x + (1 | g)) + gaussian(), data = dd,
           control = frmtmb_control(optimizer = "optim",
                                    optCtrl = list(maxit = 500)))
fit$opt$convergence
#> [1] 0

# a tighter gradient criterion, with restarts from the current
# optimum until it is met
frm(bf(y ~ x + (1 | g)) + gaussian(), data = dd,
    control = frmtmb_control(grad_tol = 1e-4, restarts = 3))
#>  Family: gaussian 
#>  Links: mu = identity; sigma = log
#> 
#> Formula: y ~ x + (1 | g) 
#>    Data: dd (Number of observations: 200) 
#>  Method: ML   logLik: -296.63   AIC: 601.26   BIC: 614.454 
#> 
#> Multilevel Hyperparameters:
#> ~g (Number of levels: 10) 
#>               Estimate Est.Error l-95% CI u-95% CI
#> sd(Intercept)        1      0.23     0.63     1.58
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> Intercept     1.44      0.32     0.81     2.07    4.46  8.1e-06
#> x             0.52      0.08     0.37     0.67    6.74  1.6e-11
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI
#> sigma     0.99      0.05     0.89     1.09

# badly scaled predictors: fit an internally standardized copy first,
# then warm-start the reported fit from it
dd$xbig <- dd$x * 1e5
frm(bf(y ~ xbig + (1 | g)) + gaussian(), data = dd,
    control = frmtmb_control(autoscale = TRUE))
#>  Family: gaussian 
#>  Links: mu = identity; sigma = log
#> 
#> Formula: y ~ xbig + (1 | g) 
#>    Data: dd (Number of observations: 200) 
#>  Method: ML   logLik: -296.63   AIC: 601.26   BIC: 614.454 
#> 
#> Multilevel Hyperparameters:
#> ~g (Number of levels: 10) 
#>               Estimate Est.Error l-95% CI u-95% CI
#> sd(Intercept)        1      0.23     0.63     1.58
#> 
#> Regression Coefficients:
#>           Estimate Est.Error l-95% CI u-95% CI z value Pr(>|z|)
#> Intercept     1.44      0.32     0.81     2.07    4.46  8.1e-06
#> xbig          0.00      0.00     0.00     0.00    6.74  1.6e-11
#> 
#> Further Distributional Parameters:
#>       Estimate Est.Error l-95% CI u-95% CI
#> sigma     0.99      0.05     0.89     1.09

# the object is a plain list, so it can be built once and reused
ctrl <- frmtmb_control(check_nlev_1 = "ignore")
ctrl$optimizer
#> [1] "nlminb"