Package {robustlmm}

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

Analysis of variance for an rlmer fit.

Description

Three modes.

Single fit anova(fit)

returns a per-term robust Wald chi-square table from vcov(fit, type = vcov_type). Each term is tested marginally as T_t = \hat\beta_t^\top \hat V_t^{-1} \hat\beta_t \sim \chi^2_{k_t} under H_0: \beta_t = 0.

Nested fits differing only in fixed effects (default test = "Wald")

robust Wald restriction test on the extra coefficients. The same V used by the single-fit table is used here; vcov_type = "sandwich" carries through.

Nested fits differing in random-effects structure

the Wald and score asymptotics are invalid (boundary problem; see Self-Liang 1987 and Koller 2026 paper 2). test = "Wald" warns and switches to the parametric bootstrap automatically; test = "boot" is the default valid path. The bootstrap simulates nsim datasets from the fitted central LMM at fit0's estimates, refits both rlmer models per replicate, and uses the quasi-deviance difference D = 2 \sum [\rho_e(r_i^{(0)}/\hat\sigma_0) - \rho_e(r_i^{(1)}/\hat\sigma_1)] as the discrepancy statistic (Heritier and Ronchetti 1994; Cantoni and Ronchetti 2001; Heritier, Cantoni, Copt and Victoria-Feser 2009 sec. 5). For the special case of a single added variance component the experimental test = "score" offers a contamination-robust alternative (see Details).

Usage

## S3 method for class 'rlmerMod'
anova(
  object,
  ...,
  test = c("Wald", "boot", "score"),
  null = c("parametric", "robust"),
  vcov_type = c("default", "sandwich"),
  ddf = c("none", "satterthwaite"),
  nsim = 1000L,
  seed = NULL,
  verbose = FALSE
)

Arguments

Details

Only pairwise comparison is implemented; chains anova(fit0, fit1, fit2, ...) are not yet supported.

Small-J caveat for the sandwich. vcov_type = "sandwich" is markedly anti-conservative at small J: in a simulation study the pairwise Wald Type-I rate reached 0.15-0.20 (vs. nominal 0.05) at J = 8 and was still elevated at J = 18, returning to nominal only by J \gtrsim 50. vcov_sandwich emits a warning for J < 20; prefer vcov_type = "default" for hypothesis tests at small J.

Subject-contamination caveat for the bootstrap (experimental). The parametric-bootstrap quasi-deviance path (test = "boot") is experimental — it is currently the only exposed variance-component path. It calibrates correctly under clean Gaussian and heavy-tailed errors but is markedly anti-conservative when a subset of subjects is contaminated. In a simulation study with 10% of subjects shifted by 5 standard errors (shift_subj), Type-I climbed to 0.13-0.24 (2.5-5x nominal) across J \in \{8, 18, 50\}. The mechanism: the bootstrap simulates from the fitted central LMM at fit0's estimates, so any bias fit0's \hat\sigma_0, \hat\theta_0 picked up from the contaminated subjects propagates into a too-narrow bootstrap null. As an automatic guard, anova(test = "boot") now inspects the null fit's own random-effects robustness weights and warns when a group is heavily downweighted (smallest weight below 0.5) — a direct signal that the fit treated that group as an outlier, so the bootstrap null may be poisoned. (The robust fit absorbs a contaminated group into a downweighted random effect, which is why cooks.distance — an influence measure — does not reliably flag it here: the downweighting that removes the bias also removes the influence.) When the warning fires, treat the p-value with caution.

Experimental robust null (null = "robust"). As a mitigation for the subject-contamination anti-conservativeness above, the bootstrap path accepts null = "robust": instead of generating the parametric bootstrap from fit0's raw estimates, it generates from a contamination-cleaned null fit. Clusters the robust fit heavily downweights (smallest per-cluster random-effect weight below a tuned threshold, the same signal the automatic guard uses) are flagged, fit0 is refitted with those clusters removed, and the bootstrap is generated from the resulting de-biased (\hat\sigma, \hat\theta, \hat\beta) (with the generator's U_b = \Lambda(\theta) rebuilt at the full cleaned \hat\theta vector). The observed discrepancy D is still taken from the original untrimmed fits, so the test uses all the data; only the bootstrap null generation is cleaned. If no cluster is flagged, more than half the clusters are flagged, trimming would leave the design rank-deficient (a dropped contrast level, a collapsed grouping factor, or a constant random-slope covariate — checked by the same nonsingular-subsampling test used by the RANSAC initial estimator), or the trimmed refit fails, the call falls back to the plain parametric null (so on clean data null = "robust" closely matches null = "parametric"); the table heading states which case applied. In simulation this reduces the contamination-induced Type-I inflation at larger J while preserving power, at the cost of mild conservatism on clean data. It is experimental, validated by simulation rather than a finite-sample theorem; the table heading notes when the robust null was applied.

Experimental robust score test (test = "score"). For the common special case where the alternative adds exactly one independent scalar variance component relative to the null — e.g. (1|g) vs (1|g) + (0 + x|g), or diagonal structures adding one component — test = "score" runs a one-sided robust score test computed from the robust null fit only and calibrated by a score-only parametric bootstrap that refits just the null model per replicate (roughly 4x cheaper than test = "boot"'s double refits; hence the smaller default nsim = 199). Per cluster j, the whitened marginal residuals \tilde r_j = V_j^{-1/2}(y_j - X_j \hat\beta_0) and the whitened tested direction v_j = V_j^{-1/2} z_j (with V_j the null fit's marginal covariance and z_j the added term's design column) give the bounded contribution s_j = (v_j^\top \psi(\tilde r_j))^2 - \kappa_1 \|v_j\|^2 (\psi the fit's rho.e psi-function, \kappa_1 = E[\psi(Z)^2]); the statistic is the self-normalised cluster sum S = \sum_j s_j / \sqrt{\sum_j (s_j - \bar s)^2}, and the one-sided p-value is (1 + \#\{S^* \ge S\})/(n_{\mathrm{eff}} + 1). Because each s_j is psi-bounded, whole-cluster contamination shifts the statistic and its bootstrap reference law together instead of inflating the test: in simulation (Gaussian balanced designs, one scalar tested component; not a theorem) contaminated-null Type-I was 0.035 vs 0.115 for test = "boot" with 10% of clusters shifted at J = 50, clean-null Type-I 0.045, and power 0.920 vs 0.900 — with no null cleaning needed (raw vs cleaned generating parameters changed at most 1-2 decisions in 600+ paired replicates, so null = "robust" is ignored and no downweighted-group warning is issued for this path). Across an adversarial sweep the pattern held: J = 18 clean/contaminated 0.055/0.015, uncentered-x designs clean/contaminated 0.045/0.075, 5% single-observation outliers 0.060. Scope: both fits on the same data with identical fixed effects and a single shared grouping factor, and the alternative adds exactly one uncorrelated scalar component (length(theta) differs by 1); anything else — multi-component or correlated-slope alternatives, crossed or nested factors — stops with an error, use test = "boot" there. Caveat: contamination aligned with the tested direction is indistinguishable from the alternative for any test with power; inspect the attached per-cluster contributions (attr(., "boot")$s_j; the largest values identify the clusters driving the statistic) and cluster diagnostics (cooks.distance(fit, groups = ), hatvalues) when a rejection is suspect. Experimental: simulation-validated, not proven; the deviance bootstrap test = "boot" remains the default variance-component path.

Value

An "anova" data.frame; the bootstrap path attaches attr(., "boot") = list(D_boot, n_fail, nsim, D_obs). The score path attaches attr(., "boot") = list(S_boot, n_fail, nsim, S_obs, s_j) with s_j the named per-cluster contributions to the observed statistic.

References

Heritier S, Ronchetti E (1994). Robust bounded-influence tests in general parametric models. JASA 89(427): 897–904.

Cantoni E, Ronchetti E (2001). Robust inference for generalized linear models. JASA 96(455): 1022–1030.

Heritier S, Cantoni E, Copt S, Victoria-Feser MP (2009). Robust Methods in Biostatistics. Wiley.

See Also

vcov, confint

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

Default Tukey bisquare psi-function with c = 4.685.

Description

Pre-built psi_func_rcpp from makeBisquarePsi with the conventional tuning constant c = 4.685 (\approx 95\% efficiency at the normal model). The bisquare redescends comparatively fast; lqqPsi is the recommended redescender.

See Also

lqqPsi, makeRobustbasePsi.

Case-weight (Hampel) influence function of the fixed effects of a fitted rlmerMod object.

Description

Returns IF(y_i) = -\hat{J}^{-1} \psi(\hat{par}; y_i) (the Hampel object built from the score value \psi), reusing the Jacobian \hat{J} formed by implicitIF. This is the right object for robustness quantities (gross-error sensitivity, breakdown), distinct from the local-shift sensitivity returned by implicitIF for a redescending \psi.

Usage

caseweightIF(fit, idx = NULL, use.expected = FALSE)

Arguments

Details

\hat{\sigma}, \hat{\theta} are held fixed (partial IF).

Value

A list with IF_beta (p \times |idx|), IF_u (q \times |idx|), IF_b, and the per-observation score-contribution columns Gb_cw, Gu_cw.

See Also

vcov_sandwich, vcov

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)

Confidence intervals for the fixed-effect coefficients of an rlmer fit.

Description

Two routes.

method = "Wald"

the default and the only path implemented in this package. Per-coefficient closed-form CI \hat\beta_k \pm z_{1-\alpha/2} \cdot SE_k with V taken from vcov(object, type = vcov_type) – so the "sandwich" option carries through to the interval.

method = "boot" or "BCa"

a thin dispatch to confintROB (Mason, Cantoni & Ghisletta 2021, 2024) with boot.type forwarded. The wrapper subsets the returned matrix to the fixed-effect rows so the shape matches the "Wald" path; variance-component CIs from confintROB are dropped here. vcov_type is not honoured on these paths (confintROB uses its own internal covariance).

Usage

## S3 method for class 'rlmerMod'
confint(
  object,
  parm = NULL,
  level = 0.95,
  method = c("Wald", "boot", "BCa"),
  vcov_type = c("default", "sandwich"),
  df = c("none", "satterthwaite"),
  boot.type = c("wild", "parametric"),
  nsim = 1000L,
  seed = NULL,
  ...
)

Arguments

Details

Guidance (Koller 2014; Mason et al. 2024). The chi-sq-p Wald limit is adequate for J \gtrsim 20 groups; the bootstrap earns its (substantial) cost mainly at smaller J. boot.type = "wild" (the default, following confintROB) is robust to misspecification of the response covariance, while "parametric" is exact under the fitted central LMM. method = "BCa" adds the bias-correction-and-acceleration adjustment to the bootstrap percentile (preferred when the bootstrap distribution is skewed).

Small-J caveat for the sandwich. vcov_type = "sandwich" under-covers at J < 20 (in a simulation study: coverage ~0.89 at J = 8 vs. nominal 0.95); the sandwich path emits a warning. For Wald CIs at small J prefer vcov_type = "default" or use method = "boot" / "BCa".

Value

A 2-column matrix with one row per selected fixed-effect coefficient, columns "<alpha/2> %" / "<1-alpha/2> %". Attributes "method", "vcov_type" (and, for the bootstrap paths, "boot.type") record the options used.

References

Mason F, Cantoni E, Ghisletta P (2021). Parametric and bootstrap-based inference for linear mixed-effects models in the presence of outliers. Methodology 17(4): 271–293.

Mason F, Cantoni E, Ghisletta P (2024). Bootstrap confidence intervals for fixed effects in mixed-effects models with outliers. Psychological Methods.

See Also

vcov, confintROB

Cook's-distance equivalent for an rlmerMod fit (per observation or per cluster).

Description

Joint Mahalanobis influence on the fitted (\hat{\beta}, \hat{\sigma}, \hat{\theta}). With groups = NULL (default) the unit is the observation: the per-observation influence vectors are stacked into a (p + 1 + L) \times n matrix W, and the result is \sqrt{w_i^T V^{-1} w_i} with V = (1/n) W W^T. With groups set the unit is the cluster: the per-cluster influence functions \Xi = -J_{par}^{-1} S (S the per-cluster score contributions from .scoreByCluster), restricted to the (\beta, \sigma, \theta) rows, are scored the same way with V = (1/J) \Xi \Xi^T. This flags whole-cluster outliers that no single observation reveals – e.g. before relying on the bootstrap variance-component test of anova, which is anti-conservative under group contamination. If V is singular (a variance-component boundary) the Moore-Penrose pseudo-inverse is used.

Usage

## S3 method for class 'rlmerMod'
cooks.distance(model, groups = NULL, IF = NULL, ...)

Arguments

Details

The full IF computation is the expensive part; pre-compute it once via IF = implicitIF_full(fit) and pass it in if you need cooks.distance and influence together.

Value

Named numeric vector: one entry per observation (groups = NULL) or per cluster level.

See Also

implicitIF_full, influence, hatvalues

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",
  ...,
  formula,
  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_cs(datasets, postFit, datasetIndices = "all")

fitDatasets_rlmer_ar1(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_rlmer_DAStau_bisq(datasets, postFit, datasetIndices = "all")

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

fitDatasets_rlmer_ransac(
  datasets,
  postFit,
  datasetIndices = "all",
  K = 50L,
  sub_frac = 0.5
)

fitDatasets_rlmer_ransac_bisq(
  datasets,
  postFit,
  datasetIndices = "all",
  K = 50L,
  sub_frac = 0.5
)

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_cs: Fits datasets using rlmer with method DASvar, forcing a compound-symmetric (cs) random-effects covariance structure (the single random-effects term is rewritten to cs(...) regardless of how the data were generated). Requires lme4 >= 2.0-0. The same rho functions and tuning parameters are used as for fitDatasets_rlmer_DASvar.

fitDatasets_rlmer_ar1: Fits datasets using rlmer with method DASvar, forcing an autoregressive (ar1) random-effects covariance structure. Otherwise identical to fitDatasets_rlmer_cs. Requires lme4 >= 2.0-0.

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_rlmer_DAStau_bisq: Fits datasets using rlmer with method DAStau, replacing rho.e with the redescending bisquare psi (bisquarePsi, c = 4.685). The other rho-functions use the same smoothed Huber psi and tuning parameters as fitDatasets_rlmer_DAStau.

fitDatasets_rlmer_DAStau_sizeOBR: Fits datasets using rlmer with method DAStau and size_obr = TRUE, which replaces the finite-difference size weight in the block-diagonal V_b score equation with the Hampel-OBR form. The same rho-functions and tuning parameters are used as for fitDatasets_rlmer_DAStau. For diagonal V_b the size_obr argument is silently ignored.

fitDatasets_rlmer_ransac: Fits datasets using rlmer with method DAStau and a RANSAC-derived initial estimator (ransac_lme4). The number of random subsamples is K = 50 with subsample fraction sub_frac = 0.5. Same rho-functions and tuning parameters as fitDatasets_rlmer_DAStau.

fitDatasets_rlmer_ransac_bisq: Combines fitDatasets_rlmer_ransac (RANSAC init) with fitDatasets_rlmer_DAStau_bisq (bisquarePsi for rho.e). Designed to give redescending psi a starting value safely away from phony local minima.

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 Longitudinal Datasets

Description

Generate balanced longitudinal datasets with random intercepts and slopes. Subjects are observed at multiple time points with optional treatment groups. Treatment and its interaction with time are coded as contrasts relative to the first level.

Usage

generateLongitudinalDatasets(
  numberOfDatasetsToGenerate,
  numberOfSubjects,
  numberOfTimepoints,
  numberOfTreatmentLevels = 1L,
  timeRange = c(0, 1),
  errorGenerator = rnorm,
  randomEffectGenerator = rnorm,
  trueBeta = 0,
  trueSigma = 1,
  trueTheta = c(1, 0, 1),
  contamFun = NULL,
  ...
)

Arguments

Details

The generated data follows the model:

y_{ij} = \beta_0 + \beta_1 \cdot \text{time}_{ij} + \sum_{k=1}^{K-1} \beta_{1+k} \cdot \text{treatment}_{k,i} + \sum_{k=1}^{K-1} \beta_{K+k} \cdot \text{treatment}_{k,i} \cdot \text{time}_{ij} + b_{0i} + b_{1i} \cdot \text{time}_{ij} + \epsilon_{ij}

where K is the number of treatment groups, b_i = (b_{0i}, b_{1i})^T \sim N(0, \sigma^2 \Lambda \Lambda^T) with \Lambda being the lower-triangular Cholesky factor reconstructed from the theta vector.

The theta parameterization follows lme4 conventions:

• For a 2x2 random-effects covariance structure (intercept and slope), theta has 3 elements: \theta = (\lambda_{11}, \lambda_{21}, \lambda_{22})

• The Cholesky factor is: \Lambda = \begin{pmatrix} \lambda_{11} & 0 \\ \lambda_{21} & \lambda_{22} \end{pmatrix}

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 a previous call of generateLongitudinalDatasets with a smaller number of datasets to generate, see examples.

Treatment Assignment: Subjects are assigned to treatment groups in a balanced, deterministic manner. Subject i is assigned to treatment (i - 1) mod numberOfTreatmentLevels + 1.

Contamination Function: If contamFun is provided, it is called as contamFun(data, info) after the response y is computed. The info list contains:

• datasetIndex: the dataset index

• randomEffects: the random effects vector for this dataset

• errors: the error vector for this dataset

• trueBeta: as passed to generateLongitudinalDatasets

• trueSigma: as passed to generateLongitudinalDatasets

• trueTheta: as passed to generateLongitudinalDatasets

The function must return a data frame with the same structure (columns id, time, treatment, y). This allows arbitrary modifications including:

• Modifying the response y (e.g., adding outliers)

• Changing group assignments (e.g., moving subjects between treatments)

• Modifying time values

• Any combination of the above

Value

list with generators and the original arguments

Author(s)

Manuel Koller

See Also

generateAnovaDatasets, generateMixedEffectDatasets

Examples

  oneGroup <- generateLongitudinalDatasets(2, 10, 5)
  head(oneGroup$generateData(1))
  head(oneGroup$generateData(2))
  oneGroup$formula

  twoGroups <- generateLongitudinalDatasets(2, 20, 5, numberOfTreatmentLevels = 2)
  head(twoGroups$generateData(1))
  twoGroups$formula

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

  ## contamination example: add outliers to 10% of observations
  set.seed(42)
  contam <- generateLongitudinalDatasets(
    numberOfDatasetsToGenerate = 5,
    numberOfSubjects = 20,
    numberOfTimepoints = 5,
    contamFun = function(data, info) {
      n <- nrow(data)
      idx <- sample(n, size = ceiling(0.1 * n))
      data$y[idx] <- data$y[idx] + 10
      data
    }
  )
  head(contam$generateData(1))

  ## contamination example: reassign some subjects to different treatment
  set.seed(42)
  contamGroup <- generateLongitudinalDatasets(
    numberOfDatasetsToGenerate = 5,
    numberOfSubjects = 20,
    numberOfTimepoints = 5,
    numberOfTreatmentLevels = 2,
    contamFun = function(data, info) {
      ## move first subject from T1 to T2
      data$treatment[data$id == 1] <- "T2"
      data
    }
  )
  head(contamGroup$generateData(1), 10)

  ## medsim: simulation inspired by the medication dataset from confintROB
  ## Two subjects from treatment are mislabeled as control, and responses
  ## are truncated at a measurement floor of 100.
  contaminateMedsim <- function(data, info) {
    data$y <- pmax(data$y, 100)  # measurement floor
    data$treatment[data$id %in% c("2", "4")] <- "T1"
    data
  }
  set.seed(2000)
  medsim <- generateLongitudinalDatasets(
    numberOfDatasetsToGenerate = 100,
    numberOfSubjects = 60,
    numberOfTimepoints = 7,
    numberOfTreatmentLevels = 2,
    timeRange = c(0, 18),
    trueBeta = c(240, -3.11, -2.42, 4.00),
    trueSigma = sqrt(1229.93),
    trueTheta = c(1.310266, -0.07547461, 0.2147735),
    contamFun = contaminateMedsim
  )
  head(medsim$generateData(1))

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 Repeated-Measures Datasets With Structured Covariance

Description

Generate balanced repeated-measures datasets with a structured random-effect covariance.

Usage

generateRepeatedMeasuresDatasets(
  numberOfDatasetsToGenerate,
  numberOfSubjects,
  numberOfVisits,
  numberOfReplicates = 1L,
  structure = c("unstructured", "cs", "ar1", "diag"),
  marginalSd = 1,
  correlation = 0,
  trueBeta = 0,
  trueSigma = 1,
  errorGenerator = rnorm,
  randomEffectGenerator = rnorm
)

Arguments

Details

Each subject is observed once per level of a within-subject factor visit (with numberOfVisits levels), optionally replicated numberOfReplicates times, giving a numberOfVisits-dimensional random effect per subject through the term (0 + visit | subject). The random-effect covariance follows the requested structure:

"unstructured"

an arbitrary covariance (the (0 + visit | subject) term).

"cs"

compound symmetry: a common correlation between all visits (cs(0 + visit | subject)).

"ar1"

autoregressive: \mathrm{Cor}(i, j) = \code{correlation}^{|i - j|} (ar1(0 + visit | subject)).

"diag"

uncorrelated visits (diag(0 + visit | subject)).

The data are simulated from the chosen covariance using the same machinery as generateMixedEffectDatasets; the returned object has the identical interface (including generateData, formula and the fitDatasets_* compatibility) and stores the structured covariance in trueTheta. Structured covariances require lme4 >= 2.0-0.

Value

A list with the same structure as the return value of generateMixedEffectDatasets.

See Also

generateMixedEffectDatasets, generateLongitudinalDatasets

Examples

  if (packageVersion("lme4") >= "2.0.0") {
    datasets <- generateRepeatedMeasuresDatasets(
        1, numberOfSubjects = 30, numberOfVisits = 3, numberOfReplicates = 4,
        structure = "cs", marginalSd = c(2, 1.5, 1.2), correlation = 0.5)
    fit <- rlmer(datasets$formula, datasets$generateData(1), method = "DASvar")
    VarCorr(fit)
  }

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", "design.weights"),
  ...
)

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

Robust leverage (hat values) for an rlmerMod fit.

Description

The self-leverage A_{ii} of each observation in the robust, random-effect-whitened convolution that the estimator uses internally – the robust analogue of the linear mixed-model hat value. It reduces to the classical lmer leverage at the non-robust limit (rho = cPsi); at the robust default the effective degrees of freedom \sum_i A_{ii} differ as downweighting changes each observation's pull. With groups set, the per-observation leverages are summed within each cluster, giving the cluster's leverage (its effective-df contribution).

Usage

## S3 method for class 'rlmerMod'
hatvalues(model, groups = NULL, ...)

Arguments

Value

Named numeric vector: one entry per observation (groups = NULL) or per cluster level.

See Also

cooks.distance, rlmer (the design.weights argument bounds high-leverage design points)

Implicit (numerical) influence function building blocks for an rlmer fit.

Description

Computes the local-shift sensitivity \partial \hat{\beta}/\partial y (and the analogous quantity for the spherical random effects \hat{u}) by implicit differentiation of the joint (\beta, u)-scoring equations at the fit. The returned list is a building block for caseweightIF and vcov_sandwich; users who only want robust standard errors should call vcov(object, type = "sandwich").

Usage

implicitIF(fit, use.expected = FALSE)

Arguments

Details

Prior weights and offsets are supported: the score is evaluated in the same prior-weight-whitened coordinates the estimator iterates in (U_e = \mathrm{diag}(\sqrt{w}), residuals offset-corrected).

This is the partial (\sigma, \theta held fixed) version. It is not the Hampel influence function; the Hampel object is caseweightIF, which reuses the same Jacobian but with the score value \psi on the right-hand side rather than \partial_y \psi. See Koller (2026, Paper 1) for a discussion.

Value

A list with components IF_beta, IF_u, IF_b, the Jacobian blocks Jbb, Jbu, Jub, S, and intermediate rotated-design matrices used by the cluster sandwich. IF_beta is the local-shift sensitivity matrix p \times n.

Full implicit influence function (beta, u, sigma, theta) for a fitted rlmerMod object.

Description

Extends implicitIF by adding the \sigma-DAS and \theta-DAS scoring equations to the implicit-derivative linear system. The Jacobian of the full score wrt (\beta, u, \log \sigma, \theta) is computed via numDeriv::jacobian with Richardson extrapolation, reusing the converged DAS scales on the non-\theta columns (they depend on \theta only); the Jacobian wrt the response is closed-form. Only DASvar and DAStau methods are supported (DAStau additionally requires block sizes \le 2).

Usage

implicitIF_full(fit, eps = 1e-06, use.cache = TRUE)

Arguments

Details

This function is the engine behind cooks.distance / influence; users who only want Cook's distance should call those S3 methods.

The result is cached on the fit (keyed by \theta) when use.cache = TRUE, so repeated calls – and the consumers cooks.distance, the sandwich vcov, and the Satterthwaite df of summary / anova / emmeans / confint – share a single computation.

Value

A list with components IF_beta (p \times n), IF_u (q \times n), IF_sigma (1 \times n), IF_theta (L \times n), and the model-based delta-method covariance vcov_model_delta.

See Also

implicitIF, cooks.distance

Per-observation influence on (\hat{\beta}, \hat{\sigma}, \hat{\theta}) for an rlmerMod fit.

Description

Thin wrapper that returns the stacked influence matrix used by cooks.distance. Cheaper to call this once and pass the result back to cooks.distance via IF = . than to recompute the numerical Jacobian twice.

Usage

## S3 method for class 'rlmerMod'
influence(model, do.coef = TRUE, ...)

Arguments

Value

A list with IF_beta (p \times n), IF_sigma (1 \times n), IF_theta (L \times n), and the full implicitIF_full object as full.

See Also

implicitIF_full, cooks.distance

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

Default lqq (linear-quadratic-quadratic) redescending psi-function.

Description

Pre-built psi_func_rcpp from makeRobustbasePsi("lqq") with robustbase's default lqq tuning (about 95% efficiency at the normal). The lqq psi of Koller and Stahel (2011) redescends more gradually than the bisquare and is the recommended redescender, as used by robustbase's lmrob.control(setting = "KS2014"). Use it as rho.e in rlmer; a redescender needs a good initial estimate, so pair it with init = "ransac" or rlmer_ransac.

References

Koller, M. and Stahel, W. A. (2011) Sharpening Wald-type inference in robust regression for small samples. Computational Statistics & Data Analysis 55(8), 2504–2515.

See Also

makeRobustbasePsi, bisquarePsi.

Bisquare (Tukey biweight) psi function

Description

Construct a redescending Tukey bisquare psi function.

Usage

makeBisquarePsi(c = 4.685)

Arguments

Details

Returns a psi_func_rcpp object whose psi, rho, weight and derivative slots implement Tukey's bisquare (biweight) function

\psi(x) = x (1 - (x/c)^2)^2 \quad \text{for } |x| \le c,\quad 0 \text{ otherwise}

suitable for use as rho.e (or rho.b, rho.sigma.e, rho.sigma.b) in rlmer.

Since robustlmm 3.5.0 the psi, rho, weight and derivative evaluations delegate to robustbase's compiled bisquare family (Mpsi, Mwgt, Mchi), so they match lmrob exactly; the returned values are identical (to numerical tolerance) to the previous hand-coded implementation. makeBisquarePsi is the scalar-cutoff special case of the general makeRobustbasePsi.

The Fisher consistency expectations (Erho, Epsi2, EDpsi) are computed by numerical integration against the standard normal.

The bisquare redescends comparatively fast; for a redescending fit the "lqq" psi (lqqPsi), the recommended redescender of robustbase's lmrob.control(setting = "KS2014"), is generally preferable (Koller and Stahel 2011).

Value

psi_func_rcpp object usable in rho.e, rho.b.

References

Koller, M. and Stahel, W. A. (2011) Sharpening Wald-type inference in robust regression for small samples. Computational Statistics & Data Analysis 55(8), 2504–2515.

See Also

makeRobustbasePsi, lqqPsi.

Examples

  pf <- bisquarePsi
  pf@psi(c(-6, -3, 0, 3, 6))

Redescending psi-functions from robustbase (bisquare, lqq, optimal, hampel, ggw)

Description

Construct a redescending psi-function from one of robustbase's compiled Mpsi families.

Usage

makeRobustbasePsi(family = c("bisquare", "lqq", "optimal", "hampel", "ggw"),
  cc = NULL)

Arguments

Details

Returns a psi_func_rcpp object whose psi, rho, weight and derivative slots delegate to robustbase's compiled psi-function family, so the shipped functions match lmrob exactly. The psi, its derivative and the weight are Mpsi(x, cc, family, 0), Mpsi(x, cc, family, 1) and Mwgt(x, cc, family); the rho slot is the normalised Mchi(x, cc, family) * MrhoInf(cc, family), so that rho(Inf) equals the supremum of rho (the psi_func rho convention used by rlmer). The Fisher consistency expectations (Erho, Epsi2, EDpsi) are computed by numerical integration against the standard normal.

The available families are the redescenders "bisquare" (Tukey biweight), "lqq" (linear-quadratic-quadratic), "optimal", "hampel" and "ggw" (generalised Gauss-weight). When cc is NULL the family's default tuning (\approx 95\% efficiency at the normal) is taken from robustbase::.Mpsi.tuning.default(family); for "lqq", "ggw" and "hampel" this is a short numeric specification vector rather than a single cutoff (see lmrob.control).

Of the redescenders, bisquarePsi redescends comparatively fast; the "lqq" psi of Koller and Stahel (2011), used by robustbase's lmrob.control(setting = "KS2014"), is the recommended redescender and is pre-built as lqqPsi.

Value

psi_func_rcpp object usable as rho.e (or rho.b, rho.sigma.e, rho.sigma.b) in rlmer.

References

Koller, M. and Stahel, W. A. (2011) Sharpening Wald-type inference in robust regression for small samples. Computational Statistics & Data Analysis 55(8), 2504–2515.

See Also

lqqPsi, bisquarePsi, Mpsi, lmrob.control.

Examples

  pf <- makeRobustbasePsi("lqq")
  pf@psi(c(-6, -3, 0, 3, 6))

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. https://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)

Plot longitudinal data with robustness-weight colored lines

Description

Creates a visualization of longitudinal data with one facet per treatment group. Subject trajectories are colored by their robustness weight from a robust mixed-effects model fit, with darker lines indicating lower weights (potential outliers). Fixed-effect predictions are overlaid as reference lines.

Usage

plotLongitudinalBySubject(
  data,
  formula = NULL,
  idVar = "id",
  timeVar = "time",
  treatmentVar = "treatment",
  responseVar = "y",
  rlmerArgs = list(),
  lineAlpha = 0.6,
  fixedLineWidth = 1.2,
  lowColor = "black",
  highColor = "lightgray",
  fixedLineColor = "firebrick",
  fixedLinetype = "solid",
  title = NULL,
  xlab = NULL,
  ylab = NULL
)

Arguments

Details

The function fits a robust linear mixed-effects model using rlmer and extracts the robustness weights for the random effects. Subjects with low weights (shown in darker colors) are those whose random effects deviate substantially from the assumed distribution.

The fixed-effect prediction lines show the population-average trajectory for each treatment group, ignoring random effects.

Value

A ggplot object.

See Also

rlmer, generateLongitudinalDatasets

Examples

## Not run: 
  ## Using the medication dataset from confintROB
  library(confintROB)
  plotLongitudinalBySubject(
    medication,
    idVar = "id",
    treatmentVar = "treat",
    responseVar = "pos"
  )

  ## Using simulated data
  set.seed(123)
  simdat <- generateLongitudinalDatasets(
    numberOfDatasetsToGenerate = 1,
    numberOfSubjects = 40,
    numberOfTimepoints = 7,
    numberOfTreatmentLevels = 2,
    timeRange = c(0, 18),
    trueBeta = c(200, -2, -5, 3),
    trueSigma = 30
  )
  plotLongitudinalBySubject(simdat$generateData(1))

## End(Not run)

Predictions and confidence/prediction intervals for an rlmerMod fit.

Description

By default (interval = "none") returns the same numeric vector as the lme4-style point prediction (preserved byte-for-byte from the pre-existing predict method). With interval = "confidence" or "prediction", returns a data frame with columns fit, lwr, upr, se.

Usage

## S3 method for class 'rlmerMod'
predict(
  object,
  newdata = NULL,
  re.form = NULL,
  ReForm,
  REForm,
  REform,
  terms = NULL,
  type = c("link", "response"),
  allow.new.levels = FALSE,
  na.action = na.pass,
  interval = c("none", "confidence", "prediction"),
  level = 0.95,
  ...
)

Arguments

Details

The fixed-effect contribution to the SE uses the robust cluster sandwich vcov(object, type = "sandwich"); the random-effects contribution (when REs are part of the prediction) is computed from the partial influence function of \hat{u} (the local-shift sensitivity). For interval = "prediction" the additional residual variance \hat{\sigma}^2 is added. Intervals are fit +/- z * se with z = \Phi^{-1}((1 + level)/2); for small numbers of clusters the normal approximation may under-cover, and a bootstrap CI (via confint(..., method = "boot")) is preferable.

Value

A numeric vector (when interval = "none") or a data frame with columns fit, lwr, upr, se.

See Also

vcov

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

Redescender basin (support-preservation) radius

Description

Support-preservation (basin) radius for a fixed-effects start.

Usage

ransac_basin_radius(object, cc = NULL, rho = NULL)

Arguments

Details

Computes the radius r^\star(c) = c\,\sigma / (2 \max_j \|x_j\|) of the ball around the initial fixed-effects estimate within which every observation keeps its redescender-support membership, so the population Hessian stays positive definite (Koller and Stahel; the RANSAC-RSE basin theorem). Here c is the rejection point of the redescending \psi — the smallest x > 0 with \psi(x) = 0. For the bisquare this is the tuning cutoff (default 4.685); the geometry generalises to any finite-rejection-point redescender (e.g. lqqPsi) by finding its rejection point numerically. A redescending \psi is safe to engage from a start only if the eventual estimate stays within this radius of the (high-breakdown) start — otherwise the iteration may leave the basin and converge to a phony solution.

Value

the basin radius r^\star(c) in the units of \beta (a scalar), with attribute "max_xnorm".

RANSAC initial estimator for LMM

Description

RANSAC-style random subsample initial estimator for linear mixed-effects models.

Usage

ransac_lme4(
  formula,
  data,
  K = 200L,
  sub_frac = 0.5,
  scale_fn = robustbase::Qn,
  adaptive = TRUE,
  patience = 50L,
  K_min = 50L,
  tol = 0.001,
  stratify = TRUE,
  n_keep = 1L,
  seed = NULL,
  verbose = FALSE
)

Arguments

Details

For K random subsamples of the data, fit a classical lmer on each subsample, score by a robust scale of residuals computed on the full data, and return the lmer fit minimising that score.

The motivation is that for redescending psi-functions (e.g. lqqPsi, the recommended redescender, or the faster-redescending bisquarePsi) the rlmer optimiser benefits from a starting value close to the true parameters. A bad initial estimate can produce phony local minima (e.g. random-effects correlation pinned at +/- 1; see Koller and Stahel 2022, Section 4.4). RANSAC is a classical way of generating a high-breakdown-point initial estimate by subsampling.

Value

list with fit (best lmerMod), scale (its score), subset (its row indices), scales (length-K vector of scores; NA for draws not run under an adaptive early stop), K (the requested cap), K_used (draws actually run), n_sub, n_singular (total number of collinear candidate observations skipped by the nonsingular-subsampling draw across all subsamples; see below).

Nonsingular subsampling

With categorical predictors a rank-deficient subsample often does not make lmer error — it silently drops the aliased (all-zero) fixed-effect columns of a dropped factor level — so a degenerate candidate would otherwise enter the scale competition with the wrong number of parameters. Rather than draw and then repair, each subsample is drawn nonsingular by construction using the nonsingular-subsampling algorithm of Koller and Stahel (2017, Algorithm 1; the method behind robustbase::lmrob.control(setting = "KS2014")). The draw units (clusters when a grouping factor is available, else single observations) are randomly permuted; a Gaxpy-variant LU decomposition with partial pivoting and column skipping (implemented in C++) then walks the permuted observations and greedily selects the first p that are linearly independent, skipping any observation collinear with those already chosen. This yields a full-rank fixed-effects core whenever the full design has full column rank. Whole clusters carrying that core are the mandatory seed; further clusters are then added in the permuted order until the retained data are identifiable — (i) a full-rank fixed-effects design [guaranteed by the core], (ii) at least two observed levels of every random-effects grouping factor, and (iii) non-constant numeric random-slope variables — and finally up to the target subsample size. Because adding clusters never lowers the rank or removes a level this is monotone and terminates. n_singular counts the collinear candidate observations skipped by the LU across all draws (a measure of the collinearity encountered); it is 0 on a purely continuous, full-rank design, where the algorithm selects exactly the first p permuted observations, so the draw is a uniform random cluster subsample — statistically identical to plain random subsampling. Note that a nonsingular start does not guarantee the fit stays nonsingular through the rlmer refinement: a redescending psi can zero-weight observations back into a rank-deficient design (Koller and Stahel 2017, Remark 2), which rlmer_ransac checks for and warns about.

References

Koller, M. and Stahel, W. A. (2017) Nonsingular subsampling for regression S estimators with categorical predictors. Computational Statistics 32(2), 631–646.

Examples

  set.seed(1)
  res <- ransac_lme4(Reaction ~ Days + (Days | Subject),
                      data = sleepstudy, K = 30)
  res$scale

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)

Resolve a cluster specification for a fitted rlmerMod object.

Description

Used by vcov_sandwich. For a single grouping factor (nested design) the cluster-robust sandwich is exact; for crossed factors no single clustering nests every random-effect block, so a warning is issued and the sandwich is approximate.

Usage

resolveCluster(fit, cluster, n)

Arguments

Value

A length-n factor of cluster memberships.

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,
  size_obr = FALSE,
  design.weights = NULL
)

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.

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.

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.

Coefficient-table degrees of freedom

By default (df = "auto") summary(object) reports a Satterthwaite-type df and a Pr(>|t|) column for the fixed effects whenever computing it is cheap – either the underlying influence function is already cached on the fit (from an earlier cooks.distance, vcov sandwich, confint or summary call), or its deterministic size workload is within the cutoff getOption("robustlmm.summary.df.max", 5000). On larger fits "auto" falls back to the historic Estimate / Std. Error / t value table (no p-values, as for lmer) and prints a one-line note on how to request the df anyway. Use df = "satterthwaite" to always compute it (which may be slow on large data) or df = "none" to never compute it. The cutoff is a dimensionless function of the problem size (number of observations, parameters and the df method), so the same fit behaves identically on every machine; raise or lower robustlmm.summary.df.max to show the df on larger or only on smaller fits. See robustlmm-options for this option. Unlike lmerTest, the degrees of freedom are derived from the robust, influence-function-based covariance of the variance parameters, so they stay honest under contamination. The Satterthwaite construction follows Giesbrecht and Burns (1985) and Fai and Cornelius (1996); inserting a robust sandwich covariance into the moment-matching ratio parallels the cluster-robust approach of Bell and McCaffrey (2002) and Pustejovsky and Tipton (2018). The approximation is reliable only for a moderate number of grouping levels; summary prints a note recommending confint(object, method = "boot") when the smallest grouping factor has fewer than 20 levels. The feature supports single-factor, nested (e.g. (1 | school/class)) and crossed (e.g. (1 | subject) + (1 | item)) designs with diagonal random effects, including Mallows-weighted fits, with the default vcov. Nested designs use a one-way cluster sandwich over the coarsest grouping factor; crossed designs use a Cameron-Gelbach-Miller multiway cluster-robust covariance, projected to the nearest positive-semidefinite matrix. Designs not covered (e.g. crossed random slopes) fall back to t-values with an explanatory note. The same Satterthwaite df is used by emmeans results (emmeans, emtrends, contrast) for this class when the default vcov is in effect; a user-supplied vcov. keeps the asymptotic (Inf, z-based) degrees of freedom. When a variance component is estimated at 0 (a boundary fit), the df is computed conditional on that component being held at the boundary – it equals the df of the model with that component dropped – and summary notes this; only a genuinely non-identifiable (singular) fit suppresses the df.

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)

rlmer with RANSAC initial estimator

Description

Fit rlmer with a RANSAC-derived initial estimator.

Usage

rlmer_ransac(
  formula,
  data,
  K = 200L,
  sub_frac = 0.5,
  n_starts = 1L,
  phony_threshold = 0.99,
  seed = NULL,
  max_tries = 5L,
  ...
)

Arguments

Details

Convenience wrapper that calls ransac_lme4 to obtain a starting value and then passes it to rlmer's init argument.

With n_starts > 1 it runs a multi-start consensus: rlmer is fitted from each of the n_starts best distinct RANSAC candidate starts, and the returned fit is the lowest-residual-scale one whose random-effects covariance is interior (|\hat\rho| \le phony_threshold). This samples several basins of a redescending \psi and so recovers the interior solution when the single best start happens to fall into the phony |\hat\rho| \to 1 attractor. If every start lands phony, the best-scoring fit is returned with a warning. The per-start summary is attached as attr(fit, "consensus").

The nonsingular subsample guarantees a full-rank fixed-effects design at the start only; a redescending \psi can zero-weight observations during the refinement and collapse the positive-weight design to rank-deficiency (Koller and Stahel 2017, Remark 2). After the fit converges its e-side robustness weights are inspected and a warning is issued if the design restricted to the positively-weighted observations is rank-deficient, suggesting a positive-weight \psi (e.g. rho.e = smoothPsi) or a different start.

Value

rlmerMod object.

Examples

  
  fit <- rlmer_ransac(Reaction ~ Days + (Days | Subject),
                       data = sleepstudy, K = 30)
  

Global options consulted by robustlmm

Description

robustlmm reads a small number of global options, set with options() and queried with getOption(). With one experimental exception (the Monte-Carlo DAS-tau calibration below), none of them change the fitted estimates; they only tune optional diagnostics and the default behaviour of summary (see rlmerMod-class).

Degrees-of-freedom options

These control whether summary(object) computes the robust Satterthwaite degrees of freedom and Pr(>|t|) column by default (df = "auto"); see the “Coefficient-table degrees of freedom” section of rlmerMod-class.

robustlmm.summary.df.max

Numeric, default 5000. The size cutoff for computing the Satterthwaite df under df = "auto". The cost of the df is a deterministic dimension-only “workload” W – one O(n) score evaluation per parameter-Jacobian column: W = (p + q + 1 + L) n for method "DASvar" and W = 40\,L\,n for "DAStau", where n is the number of observations, p the number of fixed effects, q the number of random effects and L the number of variance parameters. If W exceeds this cutoff and no influence function is cached on the fit, summary falls back to the plain Estimate / Std. Error / t value table and prints a note. The rule is dimensionless, so the same fit behaves identically on every machine. The default 5000 computes the df by default up to about n = 170 for a single random intercept fit with "DASvar" (n = 125 for "DAStau"). Set it higher to show the df on larger fits, or to 0 to always skip the automatic computation (you can still request it with summary(object, df = "satterthwaite")).

Monte-Carlo DAS-tau calibration (EXPERIMENTAL)

These options enable and tune an experimental Monte-Carlo calibration of the DAS-tau fixed point in rlmer. Without it, method = "DAStau" computes the consistency factors for non-diagonal random-effect blocks by Gauss-Hermite quadrature, which is limited to blocks of dimension \le 2; fits containing a larger block fall back to method = "DASvar" with a warning. The Monte-Carlo path lifts this restriction: it computes the same self-consistent fixed point by plain Monte-Carlo integration, which works for any block dimension, including structured (cs/ar1) and unstructured blocks of dimension > 2. It is simulation-validated – in the companion ar1 simulation study it removes the small calibration residual that the DASvar approximation leaves in the fitted correlation – but it is not backed by a finite-sample theorem, and it has been validated on clean Gaussian data only. For blocks of dimension \le 2 the classical quadrature path remains the default and the Monte-Carlo path offers no improvement there.

robustlmm.dastau.mc

Logical, default FALSE. Master switch. When TRUE, method = "DAStau" uses the Monte-Carlo calibration for all non-diagonal random-effect blocks of dimension > 2 (instead of falling back to "DASvar" for the whole fit). The Monte-Carlo sample is drawn once per fit (common random numbers), deterministically seeded and moment-matched to exact zero mean and identity second moment, so repeated fits are identical and the caller's .Random.seed is left untouched. rlmer emits a message when the experimental path is active.

robustlmm.dastau.mc.all

Logical, default FALSE. Research switch: also use the Monte-Carlo path for blocks of dimension 2, replacing the Gauss-Hermite quadrature. Intended only for comparing the two calibration paths; it offers no improvement over the quadrature.

robustlmm.dasmc.nsim

Integer, default 1e5. Number of Monte-Carlo draws. Larger values reduce the (deterministic, seed-dependent) residual calibration error at linear cost in time and memory.

robustlmm.dasmc.seed

Integer, default 20260703. Seed for the common-random-numbers draw. Fits are deterministic given this option; change it (e.g. per replicate in a simulation) to decorrelate the residual Monte-Carlo calibration error across fits. The global RNG state is saved and restored around the draw.

Developer options

robustlmm.check_rhs_optimisation

Logical, default FALSE. When TRUE, rlmer cross-checks the vectorised right-hand-side computation in the block-diagonal \theta update against an explicit per-block loop and stops on any discrepancy. Intended for development and debugging only; it adds redundant work and is not needed in normal use.

See Also

rlmerMod-class, rlmer

Examples

## show the df on larger fits (raise the size cutoff)
## Not run: 
options(robustlmm.summary.df.max = 20000)

## End(Not run)

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

Variance-covariance matrix of the fixed effects of an rlmerMod fit.

Description

By default returns the same object as lme4's vcov.merMod (the linearised model-based covariance). With type = "sandwich", returns the robust cluster-sandwich vcov_sandwich.

Usage

## S3 method for class 'rlmerMod'
vcov(
  object,
  type = c("default", "sandwich"),
  cluster = NULL,
  correction = c("G1", "none"),
  ...
)

Arguments

Details

The cluster sandwich is exact for a single (nested) grouping factor; for crossed factors it is approximate (a warning is issued, see vcov_sandwich).

Small-J caveat for the sandwich. The G1 correction is necessary but not sufficient at very small J: in a simulation study sandwich CI coverage drops to ~0.89 at J = 8 (vs. nominal 0.95). vcov_sandwich emits a warning for J < 20; for inference at small J prefer type = "default" or pair the sandwich CI with a bootstrap calibration (e.g. confintROB).

Value

A p \times p covariance matrix for \hat{\beta}.

See Also

vcov_sandwich

Robust cluster-sandwich covariance of the fixed effects of a fitted rlmerMod object.

Description

Computes the robust score sandwich \hat{V}_{IF} = \hat{A}^{-1} \hat{B} \hat{A}^{-T}, where \hat{A} is the Schur-complement (marginal) Jacobian of the profiled \beta-score and \hat{B} = \sum_j s_j s_j^T sums the per-cluster \beta-score contributions s_j = \sum_{i \in j} x_i \psi_e(\hat{r}_i). Equal to the user-facing vcov(object, type = "sandwich").

Usage

vcov_sandwich(fit, cluster = NULL, correction = c("G1", "none"))

Arguments

Details

Exact for a single (nested) grouping factor; approximate for crossed factors (a warning is issued via resolveCluster). With few clusters, set correction = "G1" (default) for the J/(J-1) small-sample scaling.

Small-J caveat. The G1 correction is necessary but not sufficient at very small J: in a simulation study CI coverage drops to ~0.89 at J = 8 (vs. nominal 0.95), and Wald-style hypothesis tests using the sandwich are anti-conservative (Type-I ~3-4x nominal). The function emits a warning for J < 20. For inference at small J prefer vcov_type = "default" or pair the sandwich CI with a bootstrap calibration (e.g. confintROB); the sandwich is most useful at J \gtrsim 50.

\hat{\sigma}, \hat{\theta} are held fixed (partial sandwich); the returned variance is the leading-order fixed-effects covariance.

Value

A p \times p covariance matrix for \hat{\beta}, with dimnames from the fixed-effect coefficient names and attribute "n.clusters".

See Also

vcov, caseweightIF

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