Robust estimation of linear mixed effects models, for hierarchical nested and non-nested, e.g., crossed, datasets.
rlmerRcpp( formula, data, ..., method = c("DAStau", "DASvar"), setting, rho.e, rho.b, rho.sigma.e, rho.sigma.b, rel.tol = 1e-08, max.iter = 40 * (r + 1)^2, verbose = 0, doFit = TRUE, init ) rlmer( formula, data, ..., method = c("DAStau", "DASvar"), setting, rho.e, rho.b, rho.sigma.e, rho.sigma.b, rel.tol = 1e-08, max.iter = 40 * (r + 1)^2, verbose = 0, doFit = TRUE, init ) lmerNoFit(formula, data = NULL, ..., initTheta)
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
method to be used for estimation of theta and sigma, see Details.
a string specifying suggested choices for the arguments
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 the
(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).
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 the Robust Scoring Equations estimator for 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
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
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 (arguments
Better robustness will lead to a decrease of the efficiency. With the default
setting = "RSEn", 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
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
rho.sigma.b arguments, it is up to
the user to specify the appropriate function. See
asymptoticEfficiency for methods to find tuning parameters
that yield a given asymptotic efficiency.
For simple variance components and the residual error scale use the
psi2propII to change the tuning parameters. This is
similar to Proposal 2 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 multi-dimensional blocks of random effects modeled, e.g.,
a model with correlated random intercept and slope, (referred to as
block diagonal 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
psi2propII) is no longer necessary. To yield
estimates with the same efficiency, the tuning parameters for the
block diagonal are larger than for the simple case. Tables of tuning parameters
are given in Table 2 and 3 of the vignette (
For a more robust estimate, use
setting = "RSEn" (the default). For
higher efficiency, use
setting = "RSEa". The settings described in
the following paragraph are used when
setting = "RSEa" is specified.
For the smoothed Huber function the tuning parameters to get approximately
95% efficiency are k=1.345 for
rho.sigma.e (using the squared version). For
simple variance components, the same can be used for
rho.sigma.b. For variance components including correlation
parameters, use k=5.14 for both
rho.sigma.b. Tables of tuning parameter are given in Table 2 and 3 of
the vignette (
If custom weight functions are specified using the argument
rho.e) but the argument
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 when
setting = "RSEn" is used. To get
higher efficiency either use
setting = "RSEa" (and only set arguments
rho.b). Or specify the tuning parameters by hand
To specify separate weight functions
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.
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.
object of class rlmerMod.
Manuel Koller, with thanks to Vanda Lourenço for improvements.
## dropping of VC system.time(print(rlmer(Yield ~ (1|Batch), Dyestuff2, method="DASvar"))) ## new Rcpp implementation system.time(print(rlmerRcpp(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 ## this is the same as using rlmer's argument `setting = "RSEa"` 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) ## this is the same as using rlmer's argument `setting = "RSEa"` rlmer(Reaction ~ Days + (Days|Subject), sleepstudy, rho.sigma.e = psi2propII(smoothPsi, k = 2.28), rho.b = chgDefaults(smoothPsi, k = 5.14, s=10), rho.sigma.b = chgDefaults(smoothPsi, k = 5.14, s=10)) ## End(Not run) ## Not run: ## start from lmer's initial estimate, not its fit rlmer(Yield ~ (1|Batch), Dyestuff, init = lmerNoFit) ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.