rlmer: Robust linear mixed models

Description

Robust estimation of linear mixed effects models, for hierarchical nested and non-nested, 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

rlmer(formula, data, ..., method = "DAStau",
    rho.e = smoothPsi, rho.b = smoothPsi, rho.sigma.e,
    rho.sigma.b, rel.tol = 1e-08,
    max.iter = 40 * (r + 1)^2, verbose = 0, doFit = TRUE,
    init)
  rlmerRcpp(formula, data, ..., method = "DAStau",
    rho.e = smoothPsi, rho.b = smoothPsi, 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)

Arguments

formula

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.

data

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

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 psi2propII function to specify squared weights and custom tuning parameters.

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 psi2propII to specify squared weights and custom tuning parameters or chgDefaults for regular weights for variance components including correlation parameters.

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 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.

init

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.

initTheta

parameter to initialize theta with (optional)

Value

object of class rlmerMod.

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 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.

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 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.

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 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 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 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.

Specifying (multiple) weight functions:

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.

Examples

Run this code

# NOT RUN {
## 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
  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)
# }