Robust estimation of linear mixed effects models, for hierarchical nested and non-nested, e.g., crossed, datasets.
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,
size_obr = FALSE,
design.weights = NULL
)lmerNoFit(formula, data = NULL, ..., initTheta)
object of class rlmerMod.
a two-sided linear formula object describing the
fixed-effects part of the model, with the response on the left of a
~ operator and the terms, separated by + operators, on the
right. The vertical bar character "|" separates an expression for
a model matrix and a grouping factor.
an optional data frame containing the variables named in
formula. By default the variables are taken from the environment
from which lmer is called.
Additional parameters passed to lmer to find the initial
estimates. See lmer.
method to be used for estimation of theta and sigma, see Details.
a string specifying suggested choices for the arguments
rho.e, rho.sigma.e, rho.b and rho.sigma.b.
Use "RSEn" (the default) or "RSEa". Both use
smoothPsi for all the “rho” arguments. For
rho.sigma.e, squared robustness weights are used (see
psi2propII). "RSEn" uses the same tuning parameter as
for rho.e, which leads to higher robustness but lower efficiency.
"RSEa" adjusts the tuning parameter for higher asymptotic efficiency
which results in lower robustness (k = 2.28 for default rho.e).
For diagonal random effects covariance matrices, rho.sigma.b is
treated exactly as rho.sigma.e. For block diagonal random effects
covariance matrices (with correlation terms), regular robustness weights
are used for rho.sigma.b, not squared ones, as they're not needed.
But the tuning parameters are adjusted for both rho.b and
rho.sigma.b according to the dimensions of the blocks (for both
"RSEn" or "RSEa"). For a block of dimension 2 (e.g.,
correlated random intercept and slope) k = 5.14 is used.
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 psi2propII function to specify squared
weights and custom tuning parameters.
(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 psi2propII to specify squared weights and custom tuning
parameters or chgDefaults for regular weights for variance
components including correlation parameters.
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 doFit = FALSE the model is not fit
but instead a structure with the model matrices for the random-effects
terms is returned (used to speed up tests). When doFit = TRUE, the
default, the model is fit immediately.
optional lmerMod- or rlmerMod-object to use for starting values,
a list with elements ‘fixef’, ‘u’, ‘sigma’,
‘theta’, the string "ransac", or a function producing an
lmerMod object.
An lmerMod or rlmerMod passed here supplies the model, not only
the starting values: the working object is built from it, so
formula and data no longer have any effect. Pass such an
object only when it was fitted on the same data; a warning is issued if
its number of observations does not match nrow(data). To start
from values obtained elsewhere while keeping the model given by
formula and data, use the list form. Its components are
optional --- whatever is omitted keeps the value from the internal
lmer initialisation, which is the right choice for ‘u’
when the values come from a fit on different observations.
When init = "ransac", the high-breakdown RANSAC
start is obtained by calling ransac_lme4(formula, data) with
its default K = 200 subsamples and sub_frac = 0.5. Its
purpose is to raise the breakdown point of the starting value, which is
what makes redescending psi-functions usable: their estimating equations
are non-convex, so a low-breakdown start can be drawn to a bad local
solution. rlmer takes
‘fixef’, ‘sigma’ and ‘theta’ from that fit and fits
the model on the full data. The string form is random
(not reproducible without an outer set.seed); for fine control of
K, sub_frac, or seed, use rlmer_ransac.
logical scalar; if TRUE (default FALSE),
the size-controlling weight w_delta in the block-diagonal
variance-components scoring equation is replaced by the
Hampel-OBR (Stahel 1987 / Hampel et al. 1986) form
\(w_\tau(d) = \min(1, b_\tau / |d - s\kappa - a|)\), with
\(b_\tau\) taken from the tuning constant of
rho.sigma.b and \(a\) determined by Fisher consistency
under \(\chi^2_s\). The default finite-difference form
\(w_\delta(d) = (\psi(d) - \psi(d-s\kappa))/s\) is not
Hampel-OBR-optimal at the central model; this option recovers
asymptotic Hampel-OBR efficiency for the variance-component
magnitudes (a typically modest 1-2 percentage-point gain at
matched gross-error sensitivity). Has no effect for diagonal
\(V_b\) (block size 1).
Mallows-type design weights for robustness to
high-leverage design points: NULL (default, exact
current behaviour), a numeric vector of deterministic weights
\(\eta_i \in (0, 1]\) of length n, or "mcd" to
compute \(\eta_i = \min(1, \chi^2_{p^*}(0.975) / d_i^2)^{1/2}\)
from robust squared Mahalanobis distances \(d_i^2\) of the
non-constant fixed-effects covariates
(covMcd); the "mcd" tuning
(severity \(\gamma = 1\), cutoff \(c = 0.975\)) is the
simulation-backed default (see vignette("rlmer")). The
weights multiply the
e-side score contribution of each observation throughout the
estimator -- the \(\beta\), \(u\), \(\sigma\) and
\(\theta\) equations, the model vcov, the Satterthwaite
degrees of freedom and the influence
diagnostics -- bounding the influence of high-leverage design
points. For designs without continuous covariates "mcd"
reduces to \(\eta \equiv 1\) (the unmodified estimator). Active
design weights are supported for a single grouping factor
only; rlmer stops with an error if the model has more than
one. Leverage robustness costs a little efficiency at the clean model
(about 1% at the default tuning; see vignette("rlmer")).
parameter to initialize theta with (optional)
weights is passed through to
lmer and has the usual meaning,
\(\mathrm{Var}(e_i) = \sigma^2 / w_i\). Weighted
method = "DASvar" fits made with robustlmm 3.2-0 to 3.4-5
are wrong and should be refitted: the Design Adaptive Scale
counted the prior weights twice, so \(\hat\sigma\), the
variance components and every standard error derived from them
were biased by an amount that depended on the scale the weights
happened to be supplied on. The default method = "DAStau"
was affected only through its starting value, and in practice
only at rho.e = rho.b = cPsi. Unweighted fits are
unaffected in every version. See the 3.5.0 entry of
news(package = "robustlmm") for the magnitudes and the
measured exposure of each method.
Manuel Koller, with thanks to Vanda Lourenço for improvements.
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
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 DAS-estimate 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.
DAStau supports blocks of random effects of dimension at most 2;
for fits containing a larger block it falls back to DASvar with a
warning. An experimental option,
options(robustlmm.dastau.mc = TRUE), lifts this restriction via a
Monte-Carlo calibration of the consistency factors; it is
simulation-validated but not backed by a finite-sample theorem, see
robustlmm-options.
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. With the default
setting, 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
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 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
function 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
the function 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 (vignette("rlmer")).
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.e and
\(k=2.28\) for rho.sigma.e (using the squared version). For
simple variance components, the same can be used for rho.b and
rho.sigma.b. For variance components including correlation
parameters, use \(k=5.14\) for both rho.b and
rho.sigma.b. Tables of tuning parameter are given in Table 2 and 3 of
the vignette (vignette("rlmer")).
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 when
setting = "RSEn" is used. To get
higher efficiency either use setting = "RSEa" (and only set arguments
rho.e and rho.b). Or 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.
lmerNoFit: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.
lmer, vignette("rlmer")
## dropping of VC
system.time(print(rlmer(Yield ~ (1|Batch), Dyestuff2, method="DASvar")))
if (FALSE) {
## 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))
}
if (FALSE) {
## start from lmer's initial estimate, not its fit
rlmer(Yield ~ (1|Batch), Dyestuff, init = lmerNoFit)
}
Run the code above in your browser using DataLab