Robust linear mixed effects models

Description

robustlmm provides functions for estimating linear mixed effects models in a robust way.

The main workhorse is the function rlmer; it is implemented as direct robust analogue of the popular lmer function of the lme4 package. The two functions have similar abilities and limitations. A wide range of data structures can be modeled: mixed effects models with hierarchical as well as complete or partially crossed random effects structures are possible. While the lmer function is optimized to handle large datasets efficiently, the computations employed in the rlmer function are more complex and for this reason also more expensive to compute. The two functions have the same limitations in the support of different random effect and residual error covariance structures. Both support only diagonal and unstructured random effect covariance structures.

The robustlmm package implements most of the analysis tool chain as is customary in R. The usual functions such as summary, coef, resid, etc. are provided as long as they are applicable for this type of models (see rlmerMod-class for a full list). The functions are designed to be as similar as possible to the ones in the lme4 package to make switching between the two packages easy.

Details on the implementation and example analyses are provided in the package vignette available via vignette("rlmer") (Koller 2016).

References

Manuel Koller (2016). robustlmm: An R Package for Robust Estimation of Linear Mixed-Effects Models. Journal of Statistical Software, 75(6), 1-24. doi:10.18637/jss.v075.i06

Koller M, Stahel WA (2022). "Robust Estimation of General Linear Mixed Effects Models.” In PM Yi, PK Nordhausen (eds.), Robust and Multivariate Statistical Methods, Springer Nature Switzerland AG.

Manuel Koller (2013). Robust estimation of linear mixed models. (Doctoral dissertation, Diss., Eidgenössische Technische Hochschule ETH Zürich, Nr. 20997, 2013).

Compute Asymptotic Efficiencies

Description

asymptoticEfficiency computes the theoretical asymptotic efficiency for an M-estimator for various types of equations.

Usage

asymptoticVariance(
  psi,
  equation = c("location", "scale", "eta", "tau", "mu"),
  dimension = 1
)

asymptoticEfficiency(
  psi,
  equation = c("location", "scale", "eta", "tau", "mu"),
  dimension = 1
)

findTuningParameter(
  desiredEfficiency,
  psi,
  equation = c("location", "scale", "eta", "tau", "mu"),
  dimension = 1,
  interval = c(0.15, 50),
  ...
)

Arguments

Details

The asymptotic efficiency is defined as the ratio between the asymptotic variance of the maximum likelihood estimator and the asymptotic variance of the (M-)estimator in question.

The computations are only approximate, using numerical integration in the general case. Depending on the regularity of the psi-function, these approximations can be quite crude.

References

Maronna, R. A., Martin, R. D., Yohai, V. J., & Salibián-Barrera, M. (2019). Robust statistics: theory and methods (with R). John Wiley & Sons., equation (2.25)

Rousseeuw, P. J., Hampel, F. R., Ronchetti, E. M., & Stahel, W. A. (2011). Robust statistics: the approach based on influence functions. John Wiley & Sons., Section 5.3c, Paragraph 2 (Page 286)

Bind Generated Datasets

Description

This method can be used to bind multiple datasets generated using different random genrators into one large dataset. The underlying dataset needs to be the same.

Usage

bindDatasets(..., datasetList = list(...))

Arguments

Value

merged list with generators and the contents of the prepared dataset. See 'prepareMixedEffectDataset and generateAnovaDatasets for a description of the contents.

Author(s)

Manuel Koller

See Also

splitDatasets

Examples

  datasets1 <- generateAnovaDatasets(2, 4, 4, 4)
  datasets2 <- generateAnovaDatasets(2, 4, 4, 4)
  datasets <- bindDatasets(datasets1, datasets2)
  data <- datasets$generateData(1)
  stopifnot(data$numberOfDatasets == 4,
            all.equal(datasets2$generateData(1), datasets$generateData(3),
                      check.attributes = FALSE),
            all.equal(datasets2$sphericalRandomEffects(1), datasets$sphericalRandomEffects(3)),
            all.equal(datasets2$createXMatrix(data), datasets$createXMatrix(data)),
            all.equal(datasets2$createZMatrix(data), datasets$createZMatrix(data)))

  preparedDataset <- prepareMixedEffectDataset(Reaction ~ Days + (Days|Subject), sleepstudy)
  datasets1 <- generateMixedEffectDatasets(2, preparedDataset)
  datasets2 <- generateMixedEffectDatasets(2, preparedDataset)
  datasets <- bindDatasets(datasets1, datasets2)
  data <- datasets$generateData(1)
  stopifnot(data$numberOfDatasets == 4,
            all.equal(datasets2$generateData(1), datasets$generateData(3),
                      check.attributes = FALSE),
            all.equal(datasets2$sphericalRandomEffects(1), datasets$sphericalRandomEffects(3)),
            all.equal(datasets2$createXMatrix(data), datasets$createXMatrix(data)),
            all.equal(datasets2$createZMatrix(data), datasets$createZMatrix(data)))

Change default arguments

Description

Change the default arguments for a psi_func_rcpp object

Usage

## S4 method for signature 'psi_func_rcpp'
chgDefaults(object, ...)

Arguments

Note

Note that names of named arguments are ignored. Only the order of the arguments considered when assigning new arguments.

Examples

sPsi <- chgDefaults(smoothPsi, k=2)
curve(sPsi@psi(x), 0, 3)
curve(smoothPsi@psi(x), 0, 3, col="blue", add=TRUE)

Create comparison charts for multiple fits

Description

Use compare to quickly compare the estimated parameters of the fits of multiple lmerMod or rlmerMod objects.

Usage

compare(..., digits = 3, dnames = NULL, show.rho.functions = TRUE)

## S3 method for class 'lmerMod'
getInfo(object, ...)

## S3 method for class 'rlmerMod'
getInfo(object, ...)

## S3 method for class 'comparison.table'
xtable(
  x,
  caption = NULL,
  label = NULL,
  align = NULL,
  digits = NULL,
  display = NULL,
  ...
)

## S3 method for class 'xtable.comparison.table'
print(
  x,
  add.hlines = TRUE,
  latexify.namescol = TRUE,
  include.rownames = FALSE,
  ...
)

getInfo(object, ...)

Arguments

Details

The functions xtable.comparison.table and print.xtable.comparison.table are wrapper functions for the respective xtable and print.xtable functions.

The function getInfo is internally used to prepare object for producing a comparison chart in compare.

Value

getInfo returns a list with estimated coefficients, estimated variance components, sigma, deviance and parameter configuration used to fit.

See Also

xtable

print.xtable

Examples

## Not run: 
  fm1 <- lmer(Yield ~ (1|Batch), Dyestuff)
  fm2 <- rlmer(Yield ~ (1|Batch), Dyestuff)
  compare(fm1, fm2)
  require(xtable)
  xtable(compare(fm1, fm2))
  str(getInfo(fm1))

## End(Not run)

Create Dataset List From List of Data Objects

Description

Convert a list of datasets to a dataset list similar to the ones created by generateAnovaDatasets and generateMixedEffectDatasets.

Usage

createDatasetsFromList(
  datasetList,
  formula,
  trueBeta,
  trueSigma,
  trueTheta,
  ...
)

Arguments

Details

The returned list can be passed to processFit and to any of the fitDatasets functions. Splitting and binding of datasets using splitDatasets and bindDatasets is not supported.

Value

list that can be passed to processFit and to any of the fitDatasets functions. Only generateData is implemented, all the other functions return an error if called.

See Also

generateAnovaDatasets and generateMixedEffectDatasets

Examples

  data(sleepstudy)
  sleepstudy2 <- sleepstudy
  sleepstudy2[1, "Reaction"] <- sleepstudy2[1, "Reaction"] + 10
  fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
  datasets <- createDatasetsFromList(list(sleepstudy, sleepstudy2),
                                     formula = Reaction ~ Days + (Days|Subject),
                                     trueBeta = getME(fm1, "beta"),
                                     trueSigma = sigma(fm1),
                                     trueTheta = getME(fm1, "theta"))
  fitDatasets_lmer(datasets)

Create Rho-Functions With Custom Tuning Parameter

Description

Convenience function to create rho-functions with custom tuning parameter.

Usage

createRhoFunction(
  tuningParameter,
  which = c("rho.e", "rho.sigma.e", "rho.b.diagonal", "rho.sigma.b.diagonal",
    "rho.b.blockDiagonal", "rho.sigma.b.blockDiagonal"),
  rho.e = smoothPsi,
  rho.sigma.e = psi2propII(rho.e),
  rho.b.diagonal = rho.e,
  rho.sigma.b.diagonal = psi2propII(rho.b.diagonal),
  rho.b.blockDiagonal = rho.e,
  rho.sigma.b.blockDiagonal = rho.b.blockDiagonal,
  ...
)

Arguments

Details

'rho.b.diagonal' denotes the tuning parameter to be used for 'rho.b' for models with diagonal random effects covariance matrix. 'rho.b.blockDiagonal' is the tuning parameter to be used in the block diagonal case, respectively.

For arguments rho.sigma.e (and rho.sigma.b.diagonal), the Proposal 2 variant of the function specified for rho.e (and rho.b) is used.

Author(s)

Manuel Koller

Examples

  createRhoFunction(c(1.345, 2.28, 1.345, 2.28, 5.14, 5.14), "rho.sigma.e")

Extract Tuning Parameters Used In Fitting

Description

Methods to extract which tuning parameters have been used for fitting models. Use extractTuningParameter for custom configurations and extractPredefinedTuningParameter for predefined configurations provided in this package.

Usage

extractTuningParameter(
  tuningParameter,
  which = c("rho.e", "rho.sigma.e", "rho.b.diagonal", "rho.sigma.b.diagonal",
    "rho.b.blockDiagonal", "rho.sigma.b.blockDiagonal")
)

extractPredefinedTuningParameter(label, which)

Arguments

Value

scalar tuning parameter

Author(s)

Manuel Koller

Examples

  extractPredefinedTuningParameter("fitDatasets_rlmer_DAStau", "rho.e")

Fitting Functions

Description

Methods to fit various mixed effects estimators to all generated datasets.

Usage

fitDatasets_lmer(datasets, control, label, postFit, datasetIndices = "all")

fitDatasets_lmer_bobyqa(datasets, postFit, datasetIndices = "all")

fitDatasets_lmer_Nelder_Mead(datasets, postFit, datasetIndices = "all")

fitDatasets_rlmer(
  datasets,
  method,
  tuningParameter,
  label,
  postFit,
  datasetIndices = "all",
  ...,
  init
)

fitDatasets_rlmer_DAStau(datasets, postFit, datasetIndices = "all")

fitDatasets_rlmer_DAStau_lmerNoFit(datasets, postFit, datasetIndices = "all")

fitDatasets_rlmer_DASvar(datasets, postFit, datasetIndices = "all")

fitDatasets_rlmer_DAStau_noAdj(datasets, postFit, datasetIndices = "all")

fitDatasets_rlmer_DAStau_k_0_5(datasets, postFit, datasetIndices = "all")

fitDatasets_rlmer_DAStau_k_0_5_noAdj(datasets, postFit, datasetIndices = "all")

fitDatasets_rlmer_DAStau_k_2(datasets, postFit, datasetIndices = "all")

fitDatasets_rlmer_DAStau_k_2_noAdj(datasets, postFit, datasetIndices = "all")

fitDatasets_rlmer_DAStau_k_5(datasets, postFit, datasetIndices = "all")

fitDatasets_rlmer_DAStau_k_5_noAdj(datasets, postFit, datasetIndices = "all")

fitDatasets_heavyLme(datasets, postFit, datasetIndices = "all")

fitDatasets_lqmm(datasets, postFit, datasetIndices = "all")

fitDatasets_rlme(datasets, postFit, datasetIndices = "all")

fitDatasets_varComprob(
  datasets,
  control,
  label,
  postFit,
  datasetIndices = "all"
)

fitDatasets_varComprob_compositeTau(datasets, postFit, datasetIndices = "all")

fitDatasets_varComprob_compositeTau_OGK(
  datasets,
  postFit,
  datasetIndices = "all"
)

fitDatasets_varComprob_compositeTau_2SGS(
  datasets,
  postFit,
  datasetIndices = "all"
)

fitDatasets_varComprob_compositeS(datasets, postFit, datasetIndices = "all")

fitDatasets_varComprob_compositeS_OGK(
  datasets,
  postFit,
  datasetIndices = "all"
)

fitDatasets_varComprob_compositeS_2SGS(
  datasets,
  postFit,
  datasetIndices = "all"
)

fitDatasets_varComprob_S(datasets, postFit, datasetIndices = "all")

fitDatasets_varComprob_S_OGK(datasets, postFit, datasetIndices = "all")

fitDatasets_varComprob_S_2SGS(datasets, postFit, datasetIndices = "all")

Arguments

Details

Existing fitting functions are:

fitDatasets_lmer: Fits datasets using lmer using its default options.

fitDatasets_lmer_bobyqa: Fits datasets using lmer using the bobyqa optimizer.

fitDatasets_lmer_Nelder_Mead: Fits datasets using lmer using the Nelder Mead optimizer.

fitDatasets_rlmer: Fits datasets using rlmer using a custom configuration. The argument 'tuningParameter' is passed to extractTuningParameter, details are documented there.

fitDatasets_rlmer_DAStau: Fits datasets using rlmer using method DAStau and smoothPsi for the rho functions. The tuning parameters are k = 1.345 for rho.e. For rho.sigma.e, the Proposal 2 variant is used using k = 2.28. The choices for rho.b and rho.sigma.b depend on whether the model uses a diagonal or a block diagonal matrix for Lambda. In the former case, the same psi functions and tuning parameters are use as for rho.e and rho.sigma.b. In the block diagonal case, rho.b and rho.sigma.b both use smoothPsi using a tuning parameter k = 5.14 (assuming blocks of dimension 2).

fitDatasets_rlmer_DAStau_lmerNoFit: Fits datasets using rlmer using the same configuration as fitDatasets_rlmer_DAStau except for that it is using lmerNoFit as initial estimator.

fitDatasets_rlmer_DASvar: Fits datasets using rlmer using method DASvar. The same rho functions and tuning parameters are used as for fitDatasets_rlmer_DAStau.

fitDatasets_rlmer_DAStau_noAdj: Fits datasets using rlmer using method DAStau. The same rho functions and tuning parameters are used as for fitDatasets_rlmer_DAStau, except for rho.sigma.e (and rho.sigma.b in the diagonal case) for which the Proposal 2 variant of smoothPsi using k = 1.345 is used.

fitDatasets_rlmer_DAStau_k_0_5: Fits datasets using rlmer using method DAStau. Use smoothPsi psi-function with tuning parameter k = 0.5 for rho.e and k = 1.47 for rho.sigma.e, the latter adjusted to reach the same asymptotic efficiency. In the diagonal case, the same are used for rho.b and rho.sigma.b as well. In the block-diagonal case, the tuning parameter k = 2.17 is used for rho.b and rho.sigma.b. The tuning parameter is chosen to reach about the same asymptotic efficiency for theta as for the fixed effects.

fitDatasets_rlmer_DAStau_k_0_5_noAdj: Fits datasets using rlmer using method DAStau. Use smoothPsi psi-function with tuning parameter k = 0.5 for rho.e and rho.sigma.e. In the diagonal case, the same are used for rho.b and rho.sigma.b as well. In the block-diagonal case, the tuning parameter k = 2.17 is used for rho.b and rho.sigma.b. The tuning parameter is chosen to reach about the same asymptotic efficiency for theta as for the fixed effects.

fitDatasets_rlmer_DAStau_k_2: Fits datasets using rlmer using method DAStau. Use smoothPsi psi-function with tuning parameter k = 2 for rho.e and k = 2.9 rho.sigma.e, the latter adjusted to reach the same asymptotic efficiency. In the diagonal case, the same are used for rho.b and rho.sigma.b as well. In the block-diagonal case, the tuning parameter k = 8.44 is used for rho.b and rho.sigma.b. The tuning parameter is chosen to reach about the same asymptotic efficiency for theta as for the fixed effects.

fitDatasets_rlmer_DAStau_k_2_noAdj: Fits datasets using rlmer using method DAStau. Use smoothPsi psi-function with tuning parameter k = 2 for rho.e and rho.sigma.e. In the diagonal case, the same are used for rho.b and rho.sigma.b as well. In the block-diagonal case, the tuning parameter k = 8.44 is used for rho.b and rho.sigma.b. The tuning parameter is chosen to reach about the same asymptotic efficiency for theta as for the fixed effects.

fitDatasets_rlmer_DAStau_k_5: Fits datasets using rlmer using method DAStau. Use smoothPsi psi-function with tuning parameter k = 5 for rho.e and k = 5.03 rho.sigma.e, the latter adjusted to reach the same asymptotic efficiency. In the diagonal case, the same are used for rho.b and rho.sigma.b as well. In the block-diagonal case, the tuning parameter k = 34.21 is used for rho.b and rho.sigma.b. The tuning parameter is chosen to reach about the same asymptotic efficiency for theta as for the fixed effects.

fitDatasets_rlmer_DAStau_k_5_noAdj: Fits datasets using rlmer using method DAStau. Use smoothPsi psi-function with tuning parameter k = 5 for rho.e and rho.sigma.e. In the diagonal case, the same are used for rho.b and rho.sigma.b as well. In the block-diagonal case, the tuning parameter k = 34.21 is used for rho.b and rho.sigma.b. The tuning parameter is chosen to reach about the same asymptotic efficiency for theta as for the fixed effects.

fitDatasets_heavyLme: Fits datasets using heavyLme from package heavy. Additional required arguments are: lmeFormula, heavyLmeRandom and heavyLmeGroups. They are passed to the formula, random and groups arguments of heavyLme.

fitDatasets_lqmm: Fits datasets using lqmm from package lqmm. Additional required arguments are: lmeFormula, lqmmRandom, lqmmGroup and lqmmCovariance. They are passed to the formula, random, groups and covariance arguments of lqmm. lqmmCovariance is optional, if omitted pdDiag is used.

fitDatasets_rlme: Fits datasets using rlme from package rlme.

fitDatasets_varComprob: Prototype method to fit datasets using varComprob from package robustvarComp. Additional required items in datasets are: lmeFormula, groups, varcov and lower. They are passed to the fixed, groups, varcov and lower arguments of varComprob. The running of this method produces many warnings of the form "passing a char vector to .Fortran is not portable" which are suppressed.

fitDatasets_varComprob_compositeTau: Fits datasets with the composite Tau method using varComprob from package robustvarComp. See fitDatasets_varComprob for additional details.

fitDatasets_varComprob_compositeTau_OGK: Similar to fitDatasets_varComprob_compositeTau but using covOGK as initial covariance matrix estimator.

fitDatasets_varComprob_compositeTau_2SGS: Similar to fitDatasets_varComprob_compositeTau but using 2SGS as initial covariance matrix estimator.

fitDatasets_varComprob_compositeS: Similar to fitDatasets_varComprob_compositeTau but using method composite S.

fitDatasets_varComprob_compositeS_OGK: Similar to fitDatasets_varComprob_compositeS but using covOGK as initial covariance matrix estimator.

fitDatasets_varComprob_compositeS_2SGS: Similar to fitDatasets_varComprob_compositeS but using 2SGS as initial covariance matrix estimator.

fitDatasets_varComprob_S: Similar to fitDatasets_varComprob_compositeTau but using method S and the Rocke psi-function.

fitDatasets_varComprob_S_OGK: Similar to fitDatasets_varComprob_S but using covOGK as initial covariance matrix estimator.

fitDatasets_varComprob_S_2SGS: Similar to fitDatasets_varComprob_S but using 2SGS as initial covariance matrix estimator.

Value

list of fitted models. See also lapplyDatasets which is called internally.

Author(s)

Manuel Koller

Examples

  set.seed(1)
  oneWay <- generateAnovaDatasets(1, 1, 10, 4,
                                  lmeFormula = y ~ 1,
                                  heavyLmeRandom = ~ 1,
                                  heavyLmeGroups = ~ Var2,
                                  lqmmRandom = ~ 1,
                                  lqmmGroup = "Var2",
                                  groups = cbind(rep(1:4, each = 10), rep(1:10, 4)),
                                  varcov = matrix(1, 4, 4),
                                  lower = 0)
  fitDatasets_lmer(oneWay)
  ## call rlmer with custom arguments
  fitDatasets_rlmer_custom <- function(datasets) {
    return(fitDatasets_rlmer(datasets,
                             method = "DASvar",
                             tuningParameter = c(1.345, 2.28, 1.345, 2.28, 5.14, 5.14),
                             label = "fitDatasets_rlmer_custom"))
  }
  fitDatasets_rlmer_custom(oneWay)

Generate ANOVA type datasets

Description

Generate balanced datasets with multiple factors. All combinations of all factor variables are generated, i.e., a fully crossed dataset will be generated. numberOfReplicates specifies the number of replications per unique combination.

Usage

generateAnovaDatasets(
  numberOfDatasetsToGenerate,
  numberOfLevelsInFixedFactor,
  numberOfSubjects,
  numberOfReplicates,
  errorGenerator = rnorm,
  randomEffectGenerator = rnorm,
  trueBeta = 1,
  trueSigma = 4,
  trueTheta = 1,
  ...,
  arrange = FALSE
)

Arguments

Details

numberOfLevelsInFixedFactor can either be a scalar or a vector with the number of levels for each fixed effects group. If numberOfLevelsInFixedFactor is a scalar, the value of 1 is allowed. This can be used to generate a dataset with an intercept only. If numberOfLevelsInFixedFactor is a vector with more than one entry, then all the values need to be larger than one.

numberOfSubjects can also be a scalar of a vector with the number of levels for each variance component. Each group needs to have more than one level. The vector is sorted descending before the names are assigned. This ensures that, when running lmer, the order of the random effects does not change. lmer also sorts the random effects by decending number of levels.

In order to save memory, only the generated random effects and the errors are stored. The dataset is only created on demand when the method generateData in the returned list is evaluated.

The random variables are generated in a way that one can simulate more datasets easily. When starting from the same seed, the first generated datasets will be the same as for the a previous call of generateAnovaDatasets with a smaller number of datasets to generate, see examples.

Value

list with generators and the original arguments

Author(s)

Manuel Koller

See Also

generateMixedEffectDatasets and createDatasetsFromList

Examples

  oneWay <- generateAnovaDatasets(2, 1, 5, 4)
  head(oneWay$generateData(1))
  head(oneWay$generateData(2))
  oneWay$formula
  head(oneWay$randomEffects(1))
  head(oneWay$sphericalRandomEffects(1))
  head(oneWay$errors(1))

  twoWayFixedRandom <- generateAnovaDatasets(2, 3, 5, 4)
  head(twoWayFixedRandom$generateData(1))
  twoWayFixedRandom$formula

  twoWayRandom <- generateAnovaDatasets(2, 1, c(3, 5), 4)
  head(twoWayRandom$generateData(1))
  twoWayRandom$formula

  large <- generateAnovaDatasets(2, c(10, 15), c(20, 30), 5)
  head(large$generateData(1))
  large$formula

  ## illustration how to generate more datasets
  set.seed(1)
  datasets1 <- generateAnovaDatasets(2, 1, 5, 4)
  set.seed(1)
  datasets2 <- generateAnovaDatasets(3, 1, 5, 4)
  stopifnot(all.equal(datasets1$generateData(1), datasets2$generateData(1)),
            all.equal(datasets1$generateData(2), datasets2$generateData(2)))

Generate Mixed Effects Datasets

Description

Generates mixed effects datasets using parametric bootstrap.

Usage

generateMixedEffectDatasets(
  numberOfDatasetsToGenerate,
  preparedDataset,
  errorGenerator = rnorm,
  randomEffectGenerator = rnorm
)

Arguments

Value

list with generators and the contents of the prepared dataset. See prepareMixedEffectDataset and generateAnovaDatasets for a description of the contents.

Author(s)

Manuel Koller

See Also

generateAnovaDatasets, prepareMixedEffectDataset and createDatasetsFromList

Examples

  preparedDataset <- prepareMixedEffectDataset(Reaction ~ Days + (Days|Subject), sleepstudy)
  datasets <- generateMixedEffectDatasets(2, preparedDataset)
  head(datasets$generateData(1))
  head(datasets$generateData(2))
  datasets$formula
  head(datasets$randomEffects(1))
  head(datasets$sphericalRandomEffects(1))
  head(datasets$errors(1))

Generate Datasets To Create Sensitivity Curves

Description

This method creates a list of datasets that can be used to create sensitivity curves. The response of the dataset is modified according to the supplied arguments.

Usage

generateSensitivityCurveDatasets(
  data,
  observationsToChange,
  shifts,
  scales,
  center,
  formula,
  ...
)

Arguments

Details

Either shifts or scales need to be provided. Both are also possible.

The argument shifts contains all the values that shall be added to each of the observations that should be changed. One value per generated dataset.

The argument scales contains all the values that shall be used to move observations away from their center. If scales is provided, then observationsToChange needs to select more than one observation.

The returned list can be passed to processFit and to any of the fitDatasets functions. Splitting and binding of datasets using splitDatasets and bindDatasets is not supported.

Value

list that can be passed to processFit and to any of the fitDatasets functions. Only generateData is implemented, all the other functions return an error if called.

See Also

generateAnovaDatasets

Examples

  oneWay <- generateAnovaDatasets(1, 1, 10, 5)
  datasets <-
      generateSensitivityCurveDatasets(oneWay$generateData(1),
                                       observationsToChange = 1:5,
                                       shifts = -10:10,
                                       formula = oneWay$formula)
  datasets$generateData(1)

Extract or Get Generalize Components from a Fitted Mixed Effects Model

Description

Extract (or “get”) “components” – in a generalized sense – from a fitted mixed-effects model, i.e. from an object of class rlmerMod or merMod.

Usage

## S3 method for class 'rlmerMod'
getME(
  object,
  name = c("X", "Z", "Zt", "Ztlist", "mmList", "y", "mu", "u", "b.s", "b", "Gp", "Tp",
    "Lambda", "Lambdat", "Tlist", "A", "U_b", "Lind", "sigma", "flist", "fixef", "beta",
    "theta", "ST", "is_REML", "n_rtrms", "n_rfacs", "N", "n", "p", "q", "p_i", "l_i",
    "q_i", "k", "m_i", "m", "cnms", "devcomp", "offset", "lower", "rho_e", "rho_b",
    "rho_sigma_e", "rho_sigma_b", "M", "w_e", "w_b", "w_b_vector", "w_sigma_e",
    "w_sigma_b", "w_sigma_b_vector"),
  ...
)

theta(object)

Arguments

Details

The function theta is short for getME(, "theta").

The goal is to provide “everything a user may want” from a fitted rlmerMod object as far as it is not available by methods, such as fixef, ranef, vcov, etc.

Value

Unspecified, as very much depending on the name.

See Also

getCall(); more standard methods for rlmerMod objects, such as ranef, fixef, vcov, etc.: see methods(class="rlmerMod")

Examples

## shows many methods you should consider *before* using getME():
methods(class = "rlmerMod")

## doFit = FALSE to speed up example
(fm1 <- rlmer(Reaction ~ Days + (Days|Subject), sleepstudy,
              method="DASvar", doFit=FALSE))
Z <- getME(fm1, "Z")
stopifnot(is(Z, "CsparseMatrix"),
          c(180,36) == dim(Z),
	  all.equal(fixef(fm1), b1 <- getME(fm1, "beta"),
		    check.attributes=FALSE, tolerance = 0))

## A way to get *all* getME()s :
## internal consistency check ensuring that all work:
parts <- getME(fm1, "ALL")
str(parts, max=2)
stopifnot(identical(Z,  parts $ Z),
          identical(b1, parts $ beta))
stopifnot(all.equal(theta(fm1), getME(fm1, "theta")))

Lapply for generated datasets

Description

Apply function for all generated datasets.

Usage

lapplyDatasets(datasets, FUN, ..., label, POST_FUN, datasetIndices = "all")

Arguments

Value

list of results. The items in the resulting list will have two additional attributes: datasetIndex and proc.time. If FUN failed for an item, then the item will be the error as returned by try, i.e., it ill be of class try-error.

Author(s)

Manuel Koller

Examples

  oneWay <- generateAnovaDatasets(2, 1, 5, 4)
  lapplyDatasets(oneWay, function(data) sum(data$y))
  lapplyDatasets(oneWay, function(data) sum(data$y), POST_FUN = function(x) x^2)

Load And Merge Partial Results

Description

Convenience function that loads the results stored in each of the files and then calls mergeProcessedFits to merge them.

Usage

loadAndMergePartialResults(files)

Arguments

Author(s)

Manuel Koller

See Also

processDatasetsInParallel

Merge Processed Fits

Description

Combine list of processed fits into one list in matrix form.

Usage

mergeProcessedFits(processedFitList)

Arguments

Value

similar list as returned by processFit just with matrix entries instead of vectors.

Examples

  preparedDataset <-
      prepareMixedEffectDataset(Reaction ~ Days + (Days|Subject),
                                sleepstudy)
  set.seed(1)
  datasets <- generateMixedEffectDatasets(2, preparedDataset)

  fits <- fitDatasets_lmer(datasets)
  processedFits <- lapply(fits, processFit, all = TRUE)
  merged <- mergeProcessedFits(processedFits)
  str(merged)

Other methods

Description

Other miscellaneous utilities for instances of the PsiFunction class.

Usage

## S4 method for signature 'Rcpp_SmoothPsi'
show(object)
## S4 method for signature 'Rcpp_HuberPsi'
show(object)
## S4 method for signature 'Rcpp_PsiFunction'
show(object)
## S4 method for signature 'Rcpp_PsiFunctionToPropIIPsiFunctionWrapper'
show(object)

Arguments

Examples

show(smoothPsi)

Compute Partial Moments

Description

Computes a partial moment for the standard normal distribution. This is the expectation taken not from -Infinity to Infinity but just to z.

Usage

partialMoment_standardNormal(z, n)

Arguments

References

Winkler, R. L., Roodman, G. M., & Britney, R. R. (1972). The Determination of Partial Moments. Management Science, 19(3), 290–296. http://www.jstor.org/stable/2629511, equation (2.5)

Examples

  partialMoment_standardNormal(0, 2)

Plot an Object of the "Psi Function" Class

Description

The plot method objects of class PsiFunction simply visualizes the \rho(), \psi(), and weight functions and their derivatives.

Usage

## S4 method for signature 'Rcpp_SmoothPsi'
plot(x, y,
     which = c("rho", "psi", "Dpsi", "wgt", "Dwgt"),
     main = "full", 
     col = c("black", "red3", "blue3", "dark green", "light green"),
     leg.loc = "right", ...)
## S4 method for signature 'Rcpp_HuberPsi'
plot(x, y,
     which = c("rho", "psi", "Dpsi", "wgt", "Dwgt"),
     main = "full", 
     col = c("black", "red3", "blue3", "dark green", "light green"),
     leg.loc = "right", ...)
## S4 method for signature 'Rcpp_PsiFunction'
plot(x, y,
     which = c("rho", "psi", "Dpsi", "wgt", "Dwgt"),
     main = "full", 
     col = c("black", "red3", "blue3", "dark green", "light green"),
     leg.loc = "right", ...)
## S4 method for signature 'Rcpp_PsiFunctionToPropIIPsiFunctionWrapper'
plot(x, y,
     which = c("rho", "psi", "Dpsi", "wgt", "Dwgt"),
     main = "full", 
     col = c("black", "red3", "blue3", "dark green", "light green"),
     leg.loc = "right", ...)

Arguments

Note

If you want to specify your own title, use main=FALSE, and a subsequent title(...) call.

See Also

psi-functions.

Examples

plot(huberPsiRcpp)
plot(huberPsiRcpp, which=c("psi", "Dpsi", "wgt"),
     main="short", leg = "topleft")

plot(smoothPsi)
## Plotting aspect ratio = 1:1 :
plot(smoothPsi, asp=1, main="short",
     which = c("psi", "Dpsi", "wgt", "Dwgt"))

Plot Method for "rlmerMod" objects.

Description

Diagnostic plots for objects of class rlmerMod and lmerMod.

Usage

## S3 method for class 'rlmerMod'
plot(
  x,
  y = NULL,
  which = 1:4,
  title = c("Fitted Values vs. Residuals", "Normal Q-Q vs. Residuals",
    "Normal Q-Q vs. Random Effects", "Scatterplot of Random Effects for Group \"%s\""),
  multiply.weights = FALSE,
  add.line = c("above", "below", "none"),
  ...
)

## S3 method for class 'rlmerMod_plots'
print(x, ask = interactive() & length(x) > 1, ...)

Arguments

Details

The robustness weights for estimating the fixed and random effects are used in the plots, e.g., the ones returned by getME(object, "w_e") and getME(object, "w_b").

Value

a list of plots of class ggplot that can be used for further modification before plotting (using print).

See Also

getME, ggplot

Examples

## Not run: 
  rfm <- rlmer(Yield ~ (1|Batch), Dyestuff)
  plot(rfm)
  fm <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
  plot.rlmerMod(fm)

## End(Not run)

Prepare Dataset for Parametric Bootstrap

Description

This function runs lmer and extracts all information needed to generate new datasets using parametric bootstrap later.

Usage

prepareMixedEffectDataset(
  formula,
  data,
  REML = TRUE,
  overrideBeta,
  overrideSigma,
  overrideTheta,
  ...
)

Arguments

Value

List that can be passed to generateMixedEffectDatasets.

Author(s)

Manuel Koller

Examples

  preparedDataset <- prepareMixedEffectDataset(Reaction ~ Days + (Days|Subject), sleepstudy)
  str(preparedDataset)

Process Datasets in Parallel

Description

Convenience function to run simulation study in parallel on a single machine.

Usage

processDatasetsInParallel(
  datasets,
  path,
  baseFilename,
  fittingFunctions,
  chunkSize,
  saveFitted = FALSE,
  checkProcessed = FALSE,
  createMinimalSaveFile = FALSE,
  ncores = 1,
  clusterType = "PSOCK",
  ...
)

Arguments

Details

The merged results are saved in a file taking the name <path>/<baseFilename>-processed.Rdata. You can delete the intermediate result files with the numbers (the chunk index) in the name.

To run on multiple machines, use saveDatasets to save datasets into multiple files. Then call processFile on each of them on the designated machine. Finally, load and merge the results together using loadAndMergePartialResults.

Value

The list of all processed results merged together.

To help reproduciblility, the output of toLatex(sessionInfo(), locale = FALSE) is stored in the sessionInfo attribute.

Author(s)

Manuel Koller

See Also

saveDatasets, processFile

Process File of Stored Datasets

Description

Call this function for each file stored using saveDatasets. If a file hasn't been processed yet, then it is processed and a new file with the postfix “processed” is created containing the results.

Usage

processFile(
  file,
  fittingFunctions,
  saveFitted = FALSE,
  checkProcessed = FALSE,
  createMinimalSaveFile = FALSE,
  datasets,
  ...
)

Arguments

Details

In case the raw fits may have to be inspected or processFit may be called with another set of arguments, then set saveFitted to TRUE. In that case, another file with the postfix “fitted” is created. Remove the files with postfix “processed” and run processFile again. The fits will not be re-done but instead loaded from the file with postfix “fitted”.

Value

The list of all processed results merged together.

To help reproduciblility, the output of toLatex(sessionInfo(), locale = FALSE) is stored in the sessionInfo attribute.

Author(s)

Manuel Koller

Process Fitted Objects

Description

Methods to process fitted objects and convert into a data structure that is useful in post-processing.

Usage

processFit(
  obj,
  all = FALSE,
  coefs = TRUE,
  stdErrors = all,
  tValues = all,
  sigma = TRUE,
  thetas = TRUE,
  b = all,
  meanB = all,
  meanAbsB = all,
  residuals = all,
  converged = TRUE,
  numWarnings = all,
  procTime = all,
  ...
)

## S3 method for class 'lmerMod'
processFit(
  obj,
  all = FALSE,
  coefs = TRUE,
  stdErrors = all,
  tValues = all,
  sigma = TRUE,
  thetas = TRUE,
  b = all,
  meanB = all,
  meanAbsB = all,
  residuals = all,
  converged = TRUE,
  numWarnings = all,
  procTime = all,
  ...
)

## S3 method for class 'rlmerMod'
processFit(
  obj,
  all = FALSE,
  coefs = TRUE,
  stdErrors = all,
  tValues = all,
  sigma = TRUE,
  thetas = TRUE,
  b = all,
  meanB = all,
  meanAbsB = all,
  residuals = all,
  converged = TRUE,
  numWarnings = all,
  procTime = all,
  ...
)

## S3 method for class 'heavyLme'
processFit(
  obj,
  all = FALSE,
  coefs = TRUE,
  stdErrors = all,
  tValues = all,
  sigma = TRUE,
  thetas = TRUE,
  b = all,
  meanB = all,
  meanAbsB = all,
  residuals = all,
  converged = TRUE,
  numWarnings = all,
  procTime = all,
  ...
)

## S3 method for class 'lqmm'
processFit(
  obj,
  all = FALSE,
  coefs = TRUE,
  stdErrors = all,
  tValues = all,
  sigma = TRUE,
  thetas = TRUE,
  b = all,
  meanB = all,
  meanAbsB = all,
  residuals = all,
  converged = TRUE,
  numWarnings = all,
  procTime = all,
  ...
)

## S3 method for class 'rlme'
processFit(
  obj,
  all = FALSE,
  coefs = TRUE,
  stdErrors = all,
  tValues = all,
  sigma = TRUE,
  thetas = TRUE,
  b = all,
  meanB = all,
  meanAbsB = all,
  residuals = all,
  converged = TRUE,
  numWarnings = all,
  procTime = all,
  ...
)

## S3 method for class 'varComprob'
processFit(
  obj,
  all = FALSE,
  coefs = TRUE,
  stdErrors = all,
  tValues = all,
  sigma = TRUE,
  thetas = TRUE,
  b = all,
  meanB = all,
  meanAbsB = all,
  residuals = all,
  converged = TRUE,
  numWarnings = all,
  procTime = all,
  isInterceptCorrelationSlopeModel,
  ...
)

Arguments

Details

Warning. processFit.varComprob uses simplistic logic to convert from the parameterisation used in the robustvarComp package to theta as used in lmer and rlmer. If there are three variance components, the code assumes that they are intercept, correlation and slope. Otherwise the code assumes that the variance components are independent. Exports b and residuals are not supported.

Value

List with extracted values, most items can be suppressed to save disk space.

Examples

  set.seed(1)
  oneWay <- generateAnovaDatasets(1, 1, 10, 4,
                                  lmeFormula = y ~ 1,
                                  heavyLmeRandom = ~ 1,
                                  heavyLmeGroups = ~ Var2,
                                  lqmmRandom = ~ 1,
                                  lqmmGroup = "Var2",
                                  groups = cbind(rep(1:4, each = 10), rep(1:10, 4)),
                                  varcov = matrix(1, 4, 4),
                                  lower = 0)
  processFit(fitDatasets_lmer(oneWay)[[1]], all = TRUE)
  processFit(fitDatasets_rlmer_DASvar(oneWay)[[1]], all = TRUE)
  ## Not run: 
    processFit(fitDatasets_heavyLme(oneWay)[[1]], all = TRUE)
  
## End(Not run)
  if (require(lqmm)) {
    processFit(fitDatasets_lqmm(oneWay)[[1]], all = TRUE)
  }
  ## Not run: 
    processFit(fitDatasets_varComprob_compositeTau(oneWay)[[1]], all = TRUE)
  
## End(Not run)

Classical, Huber and smoothed Huber psi- or rho-functions

Description

\psi-functions are used by rlmer in the estimating equations and to compute robustness weights. Change tuning parameters using chgDefaults and convert to squared robustness weights using the psi2propII function.

Usage

## see examples

Details

The “classical” \psi-function cPsi can be used to get a non-robust, i.e., classical, fit. The psi slot equals the identity function, and the rho slot equals quadratic function. Accordingly, the robustness weights will always be 1 when using cPsi.

The Huber \psi-function huberPsi is identical to the one in the package robustbase. The psi slot equals the identity function within \pm k (where k is the tuning parameter). Outside this interval it is equal to \pm k. The rho slot equals the quadratic function within \pm k and a linear function outside.

The smoothed Huber \psi-function is very similar to the regular Huber \psi-function. Instead of a sharp bend like the Huber function, the smoothed Huber function bends smoothly. The first tuning contant, k, can be compared to the tuning constant of the original Huber function. The second tuning constant, s, determines the smoothness of the bend.

See Also

chgDefaults and psi2propII for changing tuning parameters; psi_func-class for a more detailed description of the slots;

Examples

plot(cPsi)
plot(huberPsiRcpp)
plot(smoothPsi)
curve(cPsi@psi(x), 0, 3, col="blue")
curve(smoothPsi@psi(x), 0, 3, add=TRUE)
curve(huberPsiRcpp@psi(x), 0, 3, add=TRUE, col="green")

Convert to Proposal 2 weight function

Description

Converts the psi_func object into a function that corresponds to Proposal 2, i.e., a function of the squared weights. The other elements of the psi_func object are adapted accordingly.

Usage

psi2propII(object, ..., adjust = FALSE)

## S4 method for signature 'psi_func_rcpp'
psi2propII(object, ..., adjust = FALSE)

Arguments

Examples

par(mfrow=c(2,1))
plot(smoothPsi)
plot(psi2propII(smoothPsi))

Get residuals

Description

The per-observation residuals are returned, i.e., the difference of the observation and the fitted value including random effects. With type one can specify whether the weights should be used or not.

Usage

## S3 method for class 'rlmerMod'
residuals(object, type = c("response", "weighted"), scaled = FALSE, ...)

Arguments

Examples

## Not run: 
  fm <- rlmer(Yield ~ (1|Batch), Dyestuff)
  stopifnot(all.equal(resid(fm, type="weighted"),
                      resid(fm) * getME(fm, "w_e")))

## End(Not run)

Robust Scoring Equations Estimator for Linear Mixed Models

Description

Robust estimation of linear mixed effects models, for hierarchical nested and non-nested, e.g., crossed, datasets.

Usage

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)

Arguments

Details

Overview:

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.

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

Recommended tuning parameters:

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

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

Value

object of class rlmerMod.

Author(s)

Manuel Koller, with thanks to Vanda Lourenço for improvements.

See Also

lmer, vignette("rlmer")

Examples

## dropping of VC
system.time(print(rlmer(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)

rlmerMod Class

Description

Class "rlmerMod" of Robustly Fitted Mixed-Effect Models

Details

A robust mixed-effects model as returned by rlmer.

Objects from the Class

Objects are created by calls to rlmer.

Methods

Almost all methods available from objects returned from lmer are also available for objects returned by rlmer. They usage is the same.

It follows a list of some the methods that are exported by this package:

• coef

• deviance (disabled, see below)

• extractAIC (disabled, see below)

• family

• fitted

• fixef

• formula

• getInfo

• isGLMM

• isLMM

• isNLMM

• isREML

• logLik (disabled, see below)

• model.frame

• model.matrix

• nobs

• plot

• predict

• ranef (only partially implemented)

• residuals

• sigma

• summary

• terms

• update

• VarCorr

• vcov

• weights

Disabled methods

A log likelihood or even a pseudo log likelihood is not defined for the robust estimates returned by rlmer. Methods that depend on the log likelihood are therefore not available. For this reason the methods deviance, extractAIC and logLik stop with an error if they are called.

See Also

rlmer; corresponding class in package lme4: merMod

Examples

showClass("rlmerMod")

## convert an object of type 'lmerMod' to 'rlmerMod'
## to use the methods provided by robustlmm
fm <- lmer(Yield ~ (1|Batch), Dyestuff)
rfm <- as(fm, "rlmerMod")
compare(fm, rfm)

Save datasets

Description

Saves dataset to one or more files.

Usage

saveDatasets(datasets, path = getwd(), file, chunkSize)

Arguments

Details

The file will be saved to path/filename.Rdata.

If chunkSize is not missing, the filename is interpreted as format specifier and passed onto sprintf. One argument is given, the index of the chunk.

Value

filename or vector of filenames.

Author(s)

Manuel Koller

Shorten Labels

Description

Shorten labels created by the various fitDatasets functions, for use in plotting, etc.

Usage

shortenLabelsKS2022(labels)

Arguments

Details

The labels are shortened as they are in the simulation study published in Koller and Stahel (2022).

Value

Vector of shortened labels

Author(s)

Manuel Koller

References

Koller M, Stahel WA (2022). "Robust Estimation of General Linear Mixed Effects Models.” In PM Yi, PK Nordhausen (eds.), Robust and Multivariate Statistical Methods, Springer Nature Switzerland AG.

Examples

  labels <- c("fitDatasets_lmer", "fitDatasets_rlmer_DAStau",
              "fitDatasets_rlmer_DAStau_noAdj",
              "fitDatasets_varComprob_compositeTau_OGK",
              "fitDatasets_varComprob_S_OGK",
              "fitDatasets_heavyLme",
              "fitDatasets_lqmm")
  shortenLabelsKS2022(labels)

Split Datasets Into Chunks

Description

Method that splits up dataset objects into smaller chunks, so that they can be processed separately.

Usage

splitDatasets(datasets, chunkSize = 50)

Arguments

Value

list of dataset lists with generators and the contents of the original dataset. See prepareMixedEffectDataset and generateAnovaDatasets for a description of the contents. There is one additional entry in the list:

Author(s)

Manuel Koller

See Also

bindDatasets

Examples

  oneWay <- generateAnovaDatasets(18, 1, 5, 4)
  datasetList <- splitDatasets(oneWay, 5)
  data <- datasetList[[4]]$generateData(1)
  stopifnot(all.equal(oneWay$generateData(16), datasetList[[4]]$generateData(1),
                      check.attributes = TRUE),
            all.equal(oneWay$sphericalRandomEffects(16),
                      datasetList[[4]]$sphericalRandomEffects(1)),
            all.equal(oneWay$createXMatrix(data), datasetList[[4]]$createXMatrix(data)),
            all.equal(oneWay$createZMatrix(data), datasetList[[4]]$createZMatrix(data)))

Access Simulation Study Code

Description

This is a convenience function to make it simple to access the simulation study script files that are shipped with robustlmm.

Usage

viewCopyOfSimulationStudy(
  study = c("sensitivityCurves.R", "consistencyAndEfficiencyDiagonal.R",
    "consistencyAndEfficiencyBlockDiagonal.R", "breakdown.R", "convergence.R",
    "robustnessDiagonal.R", "robustnessBlockDiagonal.R"),
  destinationPath = getwd(),
  overwrite = FALSE
)

Arguments

Details

The function creates a copy of the script file that can be safely edited without changing the original file.

Examples

## Not run: 
  viewCopyOfSimulationStudy("sensitivityCurves")

## End(Not run)