Robust estimation of linear mixed effects models, for hierarchical nested and nonnested, e.g., crossed, datasets.
The lmerNoFit
function can be used to get trivial
starting values. This is mainly used to verify the
algorithms to reproduce the fit by lmer
when starting from trivial initial values.
1 2 3 4 5 6 7 
formula 
a twosided linear formula object
describing the fixedeffects part of the model, with the
response on the left of a 
data 
an optional data frame containing the
variables named in 
... 
Additional parameters passed to lmer to find
the initial estimates. See 
method 
method to be used for estimation of theta and sigma, see Details. 
rho.e 
object of class psi_func, specifying the functions to use for the huberization of the residuals. 
rho.b 
object of class psi_func or list of such objects (see Details), specifying the functions to use for the huberization of the random effects. 
rho.sigma.e 
object of class psi_func, specifying
the weight functions to use for the huberization of the
residuals when estimating the variance components, use
the 
rho.sigma.b 
(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

rel.tol 
relative tolerance used as criteria in the fitting process. 
max.iter 
maximum number of iterations allowed. 
verbose 
verbosity of output. Ranges from 0 (none) to 3 (a lot of output) 
doFit 
logical scalar. When 
init 
optional lmerMod or rlmerModobject to use for starting values, a list with elements ‘fixef’, ‘u’, ‘sigma’, ‘theta’, or a function producing an lmerMod object. 
initTheta 
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 lme4
. The
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
vignette("rlmer")
).
Models are specified using the formula
argument,
using the same syntax as for lmer
.
Additionally, one also needs to specify what robust
scoring or weight functions are to be used (arguments
starting with 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 DASestimate in robust linear
regression.
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
(arguments rho.e
, rho.b
, rho.sigma.e
and 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
using the rho.sigma.e
and rho.sigma.b
arguments, it is up to the used to specify the
appropriate function.
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 locationscale
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
the 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.
If custom
weight functions are specified using the argument
rho.b
(rho.e
) but the argument
rho.sigma.b
(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
the psi2propII
and
chgDefaults
functions.
To specify separate weight functions rho.b
and
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.
lmer
, vignette("rlmer")
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 ~ (1Batch), Dyestuff2, method="DASvar")))
## Not run:
## Default method "DAStau"
system.time(rfm.DAStau < rlmer(Yield ~ (1Batch), Dyestuff))
summary(rfm.DAStau)
## DASvar method (faster, less accurate)
system.time(rfm.DASvar < rlmer(Yield ~ (1Batch), 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 + (1plate) + (1sample), 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 + (DaysSubject), sleepstudy,
rho.sigma.e = psi2propII(smoothPsi, k = 2.28),
rho.sigma.b = chgDefaults(smoothPsi, k = 5.11, s=10))
rlmer(Yield ~ (1Batch), Dyestuff, init = lmerNoFit)
## End(Not run)

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.