##  "ss" "fitlist"
lme4’s algorithms scale reasonably well with
the number of observations and the number of random effect levels. The
biggest bottleneck is in the number of top-level parameters,
i.e. covariance parameters for
lmer fits or
glmer fits with
length(getME(model, "theta"))], covariance and
fixed-effect parameters for
glmer fits with
lme4 does a derivative-free (by
default) nonlinear optimization step over the top-level parameters.
For this reason, “maximal” models involving interactions of factors
with several levels each (e.g.
(stimulus*primer | subject))
will be slow (as well as hard to estimate): if the two factors have
f2 levels respectively, then the
lmer fit will need to estimate
(f1*f2)*(f1*f2+1)/2 top-level parameters.
lme4 automatically constructs the random effects model
matrix (\(Z\)) as a sparse matrix. At
present it does not allow an option for a sparse fixed-effects
model matrix (\(X\)), which is useful
if the fixed-effect model includes factors with many levels. Treating
such factors as random effects instead, and using the modular framework
?modular) to fix the variance of this random effect at a
large value, will allow it to be modeled using a sparse matrix. (The
estimates will converge to the fixed-effect case in the limit as the
variance goes to infinity.)
calc.derivs = FALSE
After finding the best-fit model parameters (in most cases using
derivative-free algorithms such as Powell’s BOBYQA or
[g]lmer does a series of finite-difference
calculations to estimate the gradient and Hessian at the MLE. These are
used to try to establish whether the model has converged reliably, and
glmer) to estimate the standard deviations of the
fixed effect parameters (a less accurate approximation is used if the
Hessian estimate is not available. As currently implemented, this
2*n^2 - n + 1 additional evaluations of
the deviance, where
n is the total number of top-level
control = [g]lmerControl(calc.derivs = FALSE) to turn off
this calculation can speed up the fit, e.g.
m0 <- lmer(y ~ service * dept + (1|s) + (1|d), InstEval, control = lmerControl(calc.derivs = FALSE))
Benchmark results for this run with and without derivatives show an
approximately 20% speedup (from 54 to 43 seconds on a Linux machine with
AMD Ryzen 9 2.2 GHz processors). This is a case with only 2 top-level
parameters, but the fit took only 31 deviance function evaluations (see
m0@optinfo$feval) to converge, so the effect of the
additional 7 (\(n^2 -n +1\)) function
evaluations is noticeable.
lmer uses the “nloptwrap” optimizer by default;
glmer uses a combination of bobyqa (
stage) and Nelder_Mead. These are reasonably good choices, although
glmer fits to
nloptwrap for both
stages may be worth a try.
allFits() gives an easy way to check the timings of a
large range of optimizers:
As expected, bobyqa - both the implementation in the
[g]lmerControl(optimizer="bobyqa")] and the one in
optimizer="nloptwrap", optCtrl = list(algorithm = "NLOPT_LN_BOBYQA"]
- are fastest.
Occasionally, the default optimizer stopping tolerances are
unnecessarily strict. These tolerances are specific to each optimizer,
and can be set via the
optCtrl argument in
[g]lmerControl. To see the defaults for
## $algorithm ##  "NLOPT_LN_BOBYQA" ## ## $xtol_abs ##  1e-08 ## ## $ftol_abs ##  1e-08 ## ## $maxeval ##  1e+05
In the particular case of the
InstEval example, this
doesn’t help much - loosening the tolerances to
xtol_abs=1e-4 only saves 2
functional evaluations and a few seconds, while loosening the tolerances
still further gives convergence warnings.
There are not many options for parallelizing
Optimized BLAS does not seem to help much.
glmmTMBmay be faster than
lme4for GLMMs with large numbers of top-level parameters, especially for negative binomial models (i.e. compared to
MixedModels.jlpackage in Julia may be much faster for some problems. You do need to install Julia.