Robust linear mixed models
Description
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.
Usage
1 2 3 4 5 6 7 
Arguments
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) 
Details
 Overview:
This function implements a robust approach of fitting linear mixed effect models. It can be used much like the function
lmer
in the packagelme4
. The supported models are the same as forlmer
(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 forlmer
. Additionally, one also needs to specify what robust scoring or weight functions are to be used (arguments starting withrho.
). By default a smoothed version of the Huber function is used. Furthermore, themethod
argument can be used to speed up computations at the expense of accuracy of the results. Computation methods:
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.

 Weight functions:
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
andrho.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
andrho.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 functionpsi2propII
) 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. Specifying (multiple) weight functions:
If custom weight functions are specified using the argument
rho.b
(rho.e
) but the argumentrho.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 thepsi2propII
andchgDefaults
functions.To specify separate weight functions
rho.b
andrho.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 byVarCorr(.)
when applied to the model fitted withlmer
. A set of correlated random effects count as just one group.
Value
object of class rlmerMod.
Author(s)
Manuel Koller, with thanks to Vanda Lourenço for improvements.
See Also
lmer
, vignette("rlmer")
Examples
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)
