Robust estimation of linear mixed effects models, for hierarchical nested and non-nested, e.g., crossed, datasets.
lmerNoFit function can be used to get trivial
starting values. This is mainly used to verify the
algorithms to reproduce the fit by
when starting from trivial initial values.
1 2 3 4 5 6 7
a two-sided linear formula object
describing the fixed-effects part of the model, with the
response on the left of a
an optional data frame containing the
variables named in
Additional parameters passed to lmer to find
the initial estimates. See
method to be used for estimation of theta and sigma, see Details.
object of class psi_func, specifying the functions to use for the huberization of the residuals.
object of class psi_func or list of such objects (see Details), specifying the functions to use for the huberization of the random effects.
object of class psi_func, specifying
the weight functions to use for the huberization of the
residuals when estimating the variance components, use
(optional) object of class psi_func or
list of such objects, specifying the weight functions to
use for the huberization of the random effects when
estimating the variance components (see Details). Use
relative tolerance used as criteria in the fitting process.
maximum number of iterations allowed.
verbosity of output. Ranges from 0 (none) to 3 (a lot of output)
logical scalar. When
optional lmerMod- or rlmerMod-object to use for starting values, a list with elements ‘fixef’, ‘u’, ‘sigma’, ‘theta’, or a function producing an lmerMod object.
parameter to initialize theta with (optional)
This function implements a
robust approach of fitting linear mixed effect models. It
can be used much like the function
lmer in the package
supported models are the same as for
lmer (gaussian family only). The
robust approach used is based on the robustification of
the scoring equations and an application of the Design
Adaptive Scale approach.
Example analyses and theoretical details on the method
are available in the vignette (see
Models are specified using the
using the same syntax as for
Additionally, one also needs to specify what robust
scoring or weight functions are to be used (arguments
rho.). By default a smoothed version
of the Huber function is used. Furthermore, the
method argument can be used to speed up
computations at the expense of accuracy of the results.
Currently, there are two different methods available for fitting models. They only differ in how the consistency factors for the Design Adaptive Scale estimates are computed. Available fitting methods for theta and sigma.e:
DAStau (default): For this method, the consistency
factors are computed using numerical quadrature. This is
slower but yields more accurate results. This is the
direct analogue to the DAS-estimate in robust linear
DASvar: This method computes the
consistency factors using a direct approximation which is
faster but less accurate. For complex models with
correlated random effects with more than one correlation
term, this is the only method available.
The tuning parameters of the
weight functions “rho” can be used to adjust
robustness and efficiency of the resulting estimates
rho.sigma.b). Better robustness will lead to a
decrease of the efficiency. By default, the tuning
parameters are set to yield estimates with approximately
95% efficiency for the fixed effects. The variance
components are estimated with a lower efficiency but
better robustness properties.
One has to use different weight functions and tuning
parameters for simple variance components and for such
including correlation parameters. By default, they are
chosen appropriately to the model at hand. However, when
arguments, it is up to the used to specify the
For simple variance
components and the residual error scale use the function
psi2propII to change the tuning parameters.
The is similar to Proposal II in the location-scale
problem (i.e., using the squared robustness weights of
the location estimate for the scale estimate; otherwise
the scale estimate is not robust).
For random effects modeled with correlation
parameters (referred to as nondiagonal case below), use
chgDefaults function to change the
tuning parameters. The parameter estimation problem is
multivariate, unlike the case without correlation where
the problem was univariate. For the employed estimator,
this amounts to switching from simple scale estimates to
estimating correlation matrices. Therefore different
weight functions have to be used. Squaring of the weights
(using the function
psi2propII) is no
longer necessary. To yield estimates with the same
efficiency, the tuning parameters for the nondiagonal are
generally larger than for the simple case. As a rule of
thumb, one may use the squared tuning parameters of the
simple case for the nondiagonal case.
Tables of tuning factors are given in the vignette
vignette("rlmer")). For the smoothed Huber
function the tuning parameters to get approximately 95%
efficiency are k=2.28 for simple variance
components and k=5.11 for variance
components including correlation parameters.
weight functions are specified using the argument
rho.e) but the argument
rho.sigma.e) is missing, then
the squared weights are used for simple variance
components and the regular weights are used for variance
components including correlation parameters. The same
tuning parameters will be used, to get higher efficiency
one has to specify the tuning parameters by hand using
To specify separate weight functions
rho.sigma.b for different variance components, it
is possible to pass a list instead of a psi_func object.
The list entries correspond to the groups as shown by
VarCorr(.) when applied to the model fitted with
lmer. A set of correlated random effects count as
just one group.
object of class rlmerMod.
Manuel Koller, with thanks to Vanda Lourenço for improvements.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
## dropping of VC system.time(print(rlmer(Yield ~ (1|Batch), Dyestuff2, method="DASvar"))) ## Not run: ## Default method "DAStau" system.time(rfm.DAStau <- rlmer(Yield ~ (1|Batch), Dyestuff)) summary(rfm.DAStau) ## DASvar method (faster, less accurate) system.time(rfm.DASvar <- rlmer(Yield ~ (1|Batch), Dyestuff, method="DASvar")) ## compare the two compare(rfm.DAStau, rfm.DASvar) ## Fit variance components with higher efficiency ## psi2propII yields squared weights to get robust estimates rlmer(diameter ~ 1 + (1|plate) + (1|sample), Penicillin, rho.sigma.e = psi2propII(smoothPsi, k = 2.28), rho.sigma.b = psi2propII(smoothPsi, k = 2.28)) ## use chgDefaults for variance components including ## correlation terms (regular, non squared weights suffice) rlmer(Reaction ~ Days + (Days|Subject), sleepstudy, rho.sigma.e = psi2propII(smoothPsi, k = 2.28), rho.sigma.b = chgDefaults(smoothPsi, k = 5.11, s=10)) rlmer(Yield ~ (1|Batch), Dyestuff, init = lmerNoFit) ## End(Not run)