The lme4 allFit() idea: reuse the assembled design and refit under
each optimizer, to separate optimizer trouble from model
misspecification. Uses nlminb and optim (L-BFGS-B) always, plus
bobyqa (minqa) and NLopt L-BFGS (nloptr) when those packages are
installed.
Arguments
- fit
A
frmtmb_fit.- optimizers
Named list of optimizers (names or functions, as in
frmtmb_control()). Default: everything available.- start_from_mle
If
TRUE(default), start every optimizer at the original fit's estimates; ifFALSE, where the original fit started.- ...
Refused: an argument the method does not have is an error naming it, rather than silently changing nothing.
Value
A frmtmb_allfit object: $fits (the refits, NULL where
one errored), $errors (the error message of each refit that
failed, NA elsewhere), $times, and $original, the
log-likelihood and convergence code of the fit that was checked.
Printing it compares log-likelihoods, convergence and fixed-effect
spread.
Details
Every optimizer starts from the same point, as in lme4: by default
the original fit's estimates (lme4's start_from_mle = TRUE), or,
with start_from_mle = FALSE, the point the original fit started
from, which includes its start = and the starting values a located
prior places on a nonlinear parameter. A nonlinear model is rarely
defined at the default start of zero, so refits from there fail or
stop somewhere else, and their disagreement then says nothing about
the original fit.
The printed table compares each refit with the best log-likelihood
any fit reached, the original included. An optimizer that reports
convergence more than grad_tol (see frmtmb_control()) below that
best is marked as having stopped elsewhere, so a success code at a
worse optimum does not read as agreement. One that reaches the best
while reporting a failure code agrees, with its code shown: started
at the optimum, optim's line search can fail there (code 52).
Under a prior with a density, the table compares the value each
optimizer maximizes, logLik_penalized().
Examples
set.seed(6)
dd <- data.frame(x = rnorm(60), g = factor(rep(1:6, 10)))
dd$y <- rpois(60, exp(0.3 + 0.4 * dd$x + rnorm(6, 0, 0.4)[dd$g]))
fit <- frm(bf(y ~ x + (1 | g)) + poisson(), data = dd)
frm_allfit(fit)
#> original fit: logLik -100.9268848, code 0
#> refits started at the original fit's estimates
#>
#> optimizer logLik vs best convergence verdict seconds
#> nlminb -100.9269 -3.79e-11 0 agrees 0.01
#> optim -100.9269 -2.84e-13 0 agrees 0.01
#> bobyqa -100.9269 -4.33e-12 0 agrees 0.02
#> nloptr_lbfgs -100.9269 0.00e+00 0 agrees 0.02
#>
#> logLik spread: 3.79e-11
#> max fixed-effect spread: 4.18e-07