| Title: | Robust Linear Mixed Effects Models |
|---|---|
| Description: | Implements the Robust Scoring Equations estimator to fit linear mixed effects models robustly. Robustness is achieved by modification of the scoring equations combined with the Design Adaptive Scale approach. |
| Authors: | Manuel Koller [aut, cre] |
| Maintainer: | Manuel Koller <[email protected]> |
| License: | GPL-2 |
| Version: | 3.5.1-1 |
| Built: | 2026-10-26 01:28:06 UTC |
| Source: | https://github.com/cran/robustlmm |
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).
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).
Three modes.
anova(fit)
a per-term robust Wald table:
each term is tested with , from vcov(fit, type = vcov_type),
against (default covariance) or a small-sample
reference (sandwich; see ddf).
a robust Wald test of the extra coefficients with the same covariance.
Wald asymptotics
fail on the boundary of the parameter space (Self and Liang 1987),
so test = "Wald" switches to test = "boot": a
parametric bootstrap from fit0's estimates of the
quasi-deviance difference
(Heritier and Ronchetti 1994; Cantoni and Ronchetti 2001; Heritier
et al. 2009).
test = "score" is an experimental robust score test for
added independent variance components.
Only pairwise comparisons are implemented.
## S3 method for class 'rlmerMod' anova( object, ..., test = c("Wald", "boot", "score"), null = c("parametric", "robust"), nuisance = c("threshold", "raw"), vcov_type = c("default", "sandwich"), ddf = c("auto", "none", "satterthwaite", "cluster"), nsim = 1000L, seed = NULL, verbose = FALSE, cluster = NULL )## S3 method for class 'rlmerMod' anova( object, ..., test = c("Wald", "boot", "score"), null = c("parametric", "robust"), nuisance = c("threshold", "raw"), vcov_type = c("default", "sandwich"), ddf = c("auto", "none", "satterthwaite", "cluster"), nsim = 1000L, seed = NULL, verbose = FALSE, cluster = NULL )
object |
An |
... |
A second |
test |
|
null |
Bootstrap generator for |
nuisance |
How the |
vcov_type |
Covariance for the Wald paths, passed to
|
ddf |
Reference distribution for the Wald paths. |
nsim |
Bootstrap replicates; default 1000 for |
seed |
Optional RNG seed for the bootstrap. |
verbose |
Bootstrap progress messages. |
cluster |
For |
Sandwich covariance. With few clusters the G1 cluster sandwich
is biased down for between-cluster terms, and tests based
on it are anti-conservative. For vcov_type = "sandwich" the
default reference is therefore the bias-reduced CR2 sandwich
(vcov_sandwich; Bell and McCaffrey 2002) with the
approximate Hotelling -test of Pustejovsky and Tipton (2018),
on
degrees of freedom. It is much closer to its nominal level but can stay
somewhat anti-conservative for between-cluster terms; see the section
“Small-sample references for the cluster sandwich” of
vignette("simulationStudies"). CR2 needs every random effect
nested within one cluster; otherwise the G1 sandwich with
is used, and the table heading says which. A term
that a single cluster informs on its own (e.g. a factor level present
in one cluster) gets no reference distribution (NA), with a
warning.
Bootstrap test (experimental). test = "boot" holds its
level under error contamination but is anti-conservative when whole
groups are contaminated, because the bootstrap is generated from
estimates those groups have biased. It warns when the null fit heavily
downweights a group (smallest random-effect robustness weight below
0.5). null = "robust" generates from a refit without the heavily
downweighted groups instead (the observed statistic still uses all the
data); it falls back to the plain generator when no group or more than
half of the groups are flagged, when removing them would leave the
design rank-deficient, or when the refit fails, and the table heading
says which applied. When both fits put the tested variance on the
boundary, sits at its atom at zero and the p-value is 1. See the
ANOVA section of vignette("robustlmm-3.5.0").
Score test (experimental). test = "score" uses the
robust null fit only. Per cluster it forms the bounded
contribution , with the whitened marginal residuals,
the whitened design column of a tested component and
, and calibrates
the statistic by a parametric bootstrap that refits only the null
model. One added
component is tested with the self-normalised sum of the ;
several are tested jointly with the cone (chi-bar) statistic of the
whitened score vector, with a warning when the tested directions are
nearly collinear. Scope: both fits on the same data, with the same fixed
effects and a single grouping factor, and an alternative that adds
uncorrelated scalar components; otherwise it stops with an error.
Contamination aligned with a tested direction cannot be distinguished
from the alternative; the per-cluster contributions
(attr(., "boot")$s_j) show which clusters drive a rejection.
Simulation results: the ANOVA section of
vignette("robustlmm-3.5.0") (one component) and the section
“Several added variance components” of
vignette("simulationStudies").
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.
The Wald tests are stable. Experimental, and subject to change:
test = "score", and for test = "boot" the options
null = "robust" and nuisance = "threshold" (its
default).
Self SG, Liang KY (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association 82(398), 605–610.
Guedon T, Baey C, Kuhn E (2024). Bootstrap test procedure for variance components in nonlinear mixed effects models in the presence of nuisance parameters and a singular Fisher information matrix. arXiv:2306.10779.
Bell RM, McCaffrey DF (2002). Bias reduction in standard errors for linear regression with multi-stage samples. Survey Methodology 28(2), 169–181.
Pustejovsky JE, Tipton E (2018). Small-sample methods for cluster-robust variance estimation and hypothesis testing in fixed effects models. Journal of Business & Economic Statistics 36(4), 672–683.
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.
asymptoticEfficiency computes the theoretical asymptotic efficiency
for an M-estimator for various types of equations.
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), ... )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), ... )
psi |
object of class psi_func |
equation |
equation to base computations on. |
dimension |
dimension for the multivariate location and scale problem. |
desiredEfficiency |
scalar, specifying the desired asymptotic efficiency, needs to be between 0 and 1. |
interval |
interval in which to do the root search, passed on to
|
... |
passed on to |
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.
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)
This method can be used to bind multiple datasets generated using different random generators into one large dataset. The underlying dataset needs to be the same.
bindDatasets(..., datasetList = list(...))bindDatasets(..., datasetList = list(...))
... |
multiple datasets to be bound together |
datasetList |
list of datasets created with one of the generate dataset functions |
merged list with generators and the contents of the prepared
dataset. See 'prepareMixedEffectDataset and
generateAnovaDatasets for a description of the contents.
Manuel Koller
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)))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)))
c = 4.685.Pre-built psi_func_rcpp from makeBisquarePsi with
the conventional tuning constant c = 4.685 (
efficiency at the normal model). The bisquare redescends comparatively
fast; lqqPsi is the recommended redescender.
rlmerMod object.Returns (the Hampel
object built from the score value ), reusing the
Jacobian 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 .
caseweightIF(fit, idx = NULL, use.expected = FALSE)caseweightIF(fit, idx = NULL, use.expected = FALSE)
fit |
An |
idx |
Optional integer indices selecting observations to compute (default: all observations). |
use.expected |
Passed to |
are held fixed (partial IF).
A list with IF_beta (), IF_u
(), IF_b, and the per-observation
score-contribution columns Gb_cw, Gu_cw.
Change the default arguments for a psi_func_rcpp object
## S4 method for signature 'psi_func_rcpp' chgDefaults(object, ...)## S4 method for signature 'psi_func_rcpp' chgDefaults(object, ...)
object |
instance to convert |
... |
arguments to change |
Note that names of named arguments are ignored. Only the order of the arguments considered when assigning new arguments.
sPsi <- chgDefaults(smoothPsi, k=2) curve(sPsi@psi(x), 0, 3) curve(smoothPsi@psi(x), 0, 3, col="blue", add=TRUE)sPsi <- chgDefaults(smoothPsi, k=2) curve(sPsi@psi(x), 0, 3) curve(smoothPsi@psi(x), 0, 3, col="blue", add=TRUE)
Use compare to quickly compare the estimated parameters of the fits
of multiple lmerMod or rlmerMod objects.
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, ...)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, ...)
... |
objects to compare, or, for the |
digits |
number of digits to show in output |
dnames |
names of objects given as arguments (optional) |
show.rho.functions |
whether to show rho functions in output. |
object |
object |
x |
object of class "comparison.table" or "xtable.comparison.table" |
caption |
see |
label |
see |
align |
see |
display |
see |
add.hlines |
replace empty lines in comparison table by hlines.
Supersedes |
latexify.namescol |
replace “sigma” and “x” in the first column by latex equivalents. |
include.rownames |
include row numbers (the object returned by
|
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.
getInfo returns a list with estimated coefficients, estimated
variance components, sigma, deviance and parameter configuration used to
fit.
## 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)## 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)
Two routes.
method = "Wald"the default: with from vcov(object, type = vcov_type)
and a normal or quantile (see df).
method = "boot" or "BCa"
delegated to
confintROB (Mason, Cantoni and Ghisletta
2021, 2024), returning its fixed-effect rows; vcov_type =
"sandwich", df and cluster cannot be combined with
it.
## S3 method for class 'rlmerMod' confint( object, parm = NULL, level = 0.95, method = c("Wald", "boot", "BCa"), vcov_type = c("default", "sandwich"), df = c("auto", "none", "satterthwaite", "cluster"), boot.type = c("wild", "parametric"), nsim = 1000L, seed = NULL, cluster = NULL, ... )## S3 method for class 'rlmerMod' confint( object, parm = NULL, level = 0.95, method = c("Wald", "boot", "BCa"), vcov_type = c("default", "sandwich"), df = c("auto", "none", "satterthwaite", "cluster"), boot.type = c("wild", "parametric"), nsim = 1000L, seed = NULL, cluster = NULL, ... )
object |
An |
parm |
Either |
level |
Coverage level; default 0.95. |
method |
One of |
vcov_type |
Covariance for |
df |
Quantile for |
boot.type |
Bootstrap kind passed to |
nsim |
Bootstrap replicates; default 1000. |
seed |
Optional RNG seed for reproducibility of the bootstrap. |
cluster |
For |
... |
Additional arguments forwarded to |
The Wald limit is adequate for about 20 or more groups; the bootstrap
earns its cost mainly with fewer (Mason et al. 2024). boot.type =
"wild" is robust to misspecification of the response covariance,
"parametric" relies on the fitted central model, and
method = "BCa" corrects the percentile interval for bias and
skewness. The samples of both bootstraps are also available from
simulate.
Sandwich covariance. With few clusters, normal quantiles with
the G1 sandwich under-cover, and for between-cluster coefficients the G1
sandwich itself is biased down. For vcov_type = "sandwich" the
default (df = "auto", i.e. "cluster") is therefore the
bias-reduced CR2 sandwich of vcov_sandwich with a
Bell–McCaffrey quantile per coefficient (the degrees of freedom
are in attr(., "df")); see the section “Small-sample
references for the cluster sandwich” of vignette("simulationStudies").
When the random effects are not nested within the clusters, the G1
sandwich with quantiles is used. A coefficient that a
single cluster informs on its own (e.g. a factor level present in one
cluster) gets an NA interval, with a warning.
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; the small-J sandwich
reference adds "df" (per selected coefficient) and
"correction" ("CR2", or "G1" when CR2 is not
defined).
Stable.
Bell RM, McCaffrey DF (2002). Bias reduction in standard errors for linear regression with multi-stage samples. Survey Methodology 28(2), 169–181.
Mason F, Cantoni E, Ghisletta P (2021). Parametric and semi-parametric bootstrap-based confidence intervals for robust linear mixed models. Methodology 17(4), 271–295.
Mason F, Cantoni E, Ghisletta P (2024). Linear mixed models and latent growth curve models for group comparison studies contaminated by outliers. Psychological Methods 30(1), 155–173.
Tests for a fitted rlmer model. A
single contrast gives a test, several contrasts a joint
test. The reference distributions are those of summary,
anova and the emmeans support.
contrast_test( object, L, rhs = 0, ddf = c("auto", "satterthwaite", "cluster", "none"), vcov_type = c("default", "sandwich"), cluster = NULL, ... )contrast_test( object, L, rhs = 0, ddf = c("auto", "satterthwaite", "cluster", "none"), vcov_type = c("default", "sandwich"), cluster = NULL, ... )
object |
an |
L |
the contrasts: a numeric vector (one contrast) or a matrix with
one contrast per row. Either one column per fixed-effect coefficient,
in the order of |
rhs |
the hypothesised value of |
ddf |
the denominator degrees of freedom. |
vcov_type |
|
cluster |
for |
... |
not used. |
If the Satterthwaite degrees of freedom are not available for the fit (for example at a variance-component boundary that leaves them unidentified), the normal or chi-square reference is used, with a warning.
A data frame of class "anova". For one contrast: the
columns Estimate (of ), Std. Error,
df, t value and Pr(>|t|). For several contrasts:
one row with NumDF, DenDF, F value and
Pr(>F). The contrast matrix is in attr(., "L").
Stable.
fm <- rlmer(Reaction ~ Days + (1 | Subject), sleepstudy) ## expected reaction time on day 5 contrast_test(fm, c(1, 5)) ## is it 300 ms? contrast_test(fm, c("(Intercept)" = 1, Days = 5), rhs = 300) ## both coefficients zero, jointly contrast_test(fm, diag(2))fm <- rlmer(Reaction ~ Days + (1 | Subject), sleepstudy) ## expected reaction time on day 5 contrast_test(fm, c(1, 5)) ## is it 300 ms? contrast_test(fm, c("(Intercept)" = 1, Days = 5), rhs = 300) ## both coefficients zero, jointly contrast_test(fm, diag(2))
Joint Mahalanobis influence on the fitted . With groups = NULL (default) the
unit is the observation: the per-observation influence vectors are
stacked into a matrix , and the
result is with .
With groups set the unit is the cluster: the
per-cluster influence functions (
the per-cluster score contributions from .scoreByCluster),
restricted to the rows, are scored the
same way with . 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 is
singular (a variance-component boundary) the Moore-Penrose
pseudo-inverse is used.
## S3 method for class 'rlmerMod' cooks.distance(model, groups = NULL, IF = NULL, ...)## S3 method for class 'rlmerMod' cooks.distance(model, groups = NULL, IF = NULL, ...)
model |
An |
groups |
Cluster grouping for cluster-level Cook's distance:
|
IF |
Optional pre-computed |
... |
Currently unused. |
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.
Named numeric vector: one entry per observation
(groups = NULL) or per cluster level.
implicitIF_full,
influence,
hatvalues
Convert a list of datasets to a dataset list similar to the ones created by
generateAnovaDatasets and
generateMixedEffectDatasets.
createDatasetsFromList( datasetList, formula, trueBeta, trueSigma, trueTheta, ... )createDatasetsFromList( datasetList, formula, trueBeta, trueSigma, trueTheta, ... )
datasetList |
list of data objects, usually of type |
formula |
formula to fit the model using |
trueBeta |
scalar or vector with the true values of the fixed effects coefficients. Can be of length one in which case it will be replicated to the required length if needed. |
trueSigma |
scalar with the true value of the error scale. |
trueTheta |
scalar or vector with the true values for the variance component coefficients, not including sigma. Can be of length one in which case it will be replicated to the required length if needed. |
... |
all additional arguments are added to the returned list. |
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.
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.
generateAnovaDatasets and
generateMixedEffectDatasets
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)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)
Convenience function to create rho-functions with custom tuning parameter.
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, ... )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, ... )
tuningParameter |
argument passed on to
|
which |
string specifying which tuning parameter should be extracted. |
rho.e |
|
rho.sigma.e |
|
rho.b.diagonal |
|
rho.sigma.b.diagonal |
|
rho.b.blockDiagonal |
|
rho.sigma.b.blockDiagonal |
|
... |
passed on to |
'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.
Manuel Koller
createRhoFunction(c(1.345, 2.28, 1.345, 2.28, 5.14, 5.14), "rho.sigma.e")createRhoFunction(c(1.345, 2.28, 1.345, 2.28, 5.14, 5.14), "rho.sigma.e")
Breakdown ceiling implied by the fixed-effects design.
design_breakdown(object, max_subsets = 2e+05)design_breakdown(object, max_subsets = 2e+05)
object |
an |
max_subsets |
budget for the exact enumeration of hyperplanes;
above it a random search is used and |
Computes , the largest number of rows of the fixed-effects
model matrix that lie on one hyperplane through the origin, and from it
max_bp: the largest finite-sample breakdown point any
regression-equivariant estimator of the fixed effects can have,
(Mili and Coakley 1996);
mm_bp: the breakdown point that an MM-type refinement is
guaranteed to inherit from its start,
at (following Yohai 1987), and b_opt
, the scale constant at which that guarantee
reaches max_bp.
In general position and the ceiling is about 1/2.
Categorical predictors raise : with an intercept and a balanced
two-level factor, , the ceiling is about 1/4 and the MM
guarantee at is empty. No start, including
rlmer_ransac, can do better than this ceiling for the
fixed effects. b_opt is not used by ransac_lme4;
see the section “Ranking RANSAC candidates” of
vignette("simulationStudies").
a list with n, p, k, k_share
(), max_bp, mm_bp, b_opt and
exact.
Experimental. The interface and the results may change in a later version.
Mili L, Coakley CW (1996). Robust estimation in structured linear regression. The Annals of Statistics 24(6), 2593–2607.
Yohai VJ (1987). High breakdown-point and high efficiency robust estimates for regression. The Annals of Statistics 15(2), 642–656.
d <- data.frame(y = rnorm(40), g = gl(2, 20), x = rnorm(40)) design_breakdown(model.matrix(~ g, d)) # dummy: ceiling ~ 1/4 design_breakdown(model.matrix(~ x, d)) # general positiond <- data.frame(y = rnorm(40), g = gl(2, 20), x = rnorm(40)) design_breakdown(model.matrix(~ g, d)) # dummy: ceiling ~ 1/4 design_breakdown(model.matrix(~ x, d)) # general position
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.
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)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)
tuningParameter |
vector of tuning parameters. The vector is expected to be of length 6, containing the tuning parameters for rho.e, rho.sigma.e, rho.b.diagonal, rho.sigma.b.diagonal, rho.b.blockDiagonal and rho.sigma.b.blockDiagonal. 'rho.b.diagonal' denotes the tuning parameter to be used for 'rho.b' for models with diagonal random effects covariance matrix. Names are optional. |
which |
string specifying which tuning parameter should be extracted. |
label |
label or vector of labels in results. Only predefined labels of the form 'fitDatasets_rlmer_...' are supported (for others NA is returned). |
scalar tuning parameter
Manuel Koller
extractPredefinedTuningParameter("fitDatasets_rlmer_DAStau", "rho.e")extractPredefinedTuningParameter("fitDatasets_rlmer_DAStau", "rho.e")
Methods to fit various mixed effects estimators to all generated datasets.
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")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")
datasets |
Datasets list to be used to generate datasets. |
control |
a list (of correct class for the respective fitting function) containing control parameters to be passed through. |
label |
a string used to identify which fits have been created by which function. |
postFit |
a function, taking one argument, the resulting fit. This makes it easy to add an additional step after fitting. |
datasetIndices |
optional vector of dataset indices to fit, useful to try only a few datasets instead of all of them. |
method |
argument passed on to |
tuningParameter |
argument passed on to
|
... |
argument passed on to |
formula |
optional model formula to fit; defaults to the formula
stored in |
init |
optional argument passed on to |
K |
number of random subsamples used by the RANSAC initial
estimator ( |
sub_frac |
fraction of the data per RANSAC subsample. Only used by
|
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.
list of fitted models. See also lapplyDatasets which
is called internally.
Manuel Koller
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)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)
First-order contamination fraction that would drive each variance parameter of a robust LMM to its bound.
fragility(object, ...) ## S3 method for class 'rlmerMod' fragility(object, ...)fragility(object, ...) ## S3 method for class 'rlmerMod' fragility(object, ...)
object |
an |
... |
ignored. |
For every variance parameter with a lower bound of zero
(the diagonal entries of the relative Cholesky factor), fragility
reports
where is the influence of cluster on
at the fitted model. Replacing a fraction
of the clusters by copies of the most influential one moves
by , so it
reaches its bound at to first order: a zero variance
for the first component of a random-effects block, a singular
covariance matrix (for a block a correlation of
) for the later ones.
Only the observed clusters are searched, so is an upper
bound on the worst case: a small value flags a fragile parameter, a
large one does not certify robustness, and values above 1/2 are only
indicative. clusters, , is
the number of copies of the driving cluster that, added to the
observed clusters, make up the fraction . Crossed
designs are not supported.
a data frame with one row per bounded variance parameter:
parameter, bound, theta (estimate, relative to
the residual scale), eps (), clusters
() and driver (the most
influential cluster).
Experimental. The interface and the results may change in a later version.
fit <- rlmer(Reaction ~ Days + (Days | Subject), sleepstudy, method = "DASvar") fragility(fit)fit <- rlmer(Reaction ~ Days + (Days | Subject), sleepstudy, method = "DASvar") fragility(fit)
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.
generateAnovaDatasets( numberOfDatasetsToGenerate, numberOfLevelsInFixedFactor, numberOfSubjects, numberOfReplicates, errorGenerator = rnorm, randomEffectGenerator = rnorm, trueBeta = 1, trueSigma = 4, trueTheta = 1, ..., arrange = FALSE )generateAnovaDatasets( numberOfDatasetsToGenerate, numberOfLevelsInFixedFactor, numberOfSubjects, numberOfReplicates, errorGenerator = rnorm, randomEffectGenerator = rnorm, trueBeta = 1, trueSigma = 4, trueTheta = 1, ..., arrange = FALSE )
numberOfDatasetsToGenerate |
number of datasets to generate. |
numberOfLevelsInFixedFactor |
scalar or vector with the number of levels per fixed factor or grouping variable. |
numberOfSubjects |
scalar or vector with the number of levels per variance component. |
numberOfReplicates |
number of replicates per unique combination of fixed factor and variance component. |
errorGenerator |
random number generator used for the errors. |
randomEffectGenerator |
random number generator used for the spherical random effects. |
trueBeta |
scalar or vector with the true values of the fixed effects coefficients. Can be of length one in which case it will be replicated to the required length if needed. |
trueSigma |
scalar with the true value of the error scale. |
trueTheta |
scalar of vector with the true values for the variance component coefficients, not including sigma. Can be of length one in which case it will be replicated to the required length if needed. |
... |
all additional arguments are added to the returned list. |
arrange |
If |
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 descending 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.
list with generators and the original arguments
generateData: |
function to generate data taking one argument, the dataset index. |
createXMatrix: |
function to generate X matrix taking one argument,
the result of |
createZMatrix: |
function to generate Z matrix taking one argument,
the result of |
createLambdaMatrix: |
function to generate Lambda matrix taking one
argument, the result of |
randomEffects: |
function to return the generated random effects taking one argument, the dataset index. |
sphericalRandomeffects: |
function to return the generated spherical random effects taking one argument, the dataset index. |
errors: |
function to return the generated errors taking one argument, the dataset index. |
allRandomEffects: |
function without arguments that returns the matrix of all generated random effects. |
allErrors: |
function without arguments that returns the matrix of all generated errors. |
numberOfDatasets: |
|
numberOfLevelsInFixedFactor: |
|
numberOfSubjects: |
|
numberOfReplicates: |
|
numberOfRows: |
number of rows in the generated dataset |
trueBeta: |
true values used for beta |
trueSigma: |
true value used for sigma |
trueTheta: |
true values used for theta |
formula: |
formula to fit the model using |
...: |
additional arguments passed via |
Manuel Koller
generateMixedEffectDatasets and
createDatasetsFromList
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)))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 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.
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, ... )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, ... )
numberOfDatasetsToGenerate |
number of datasets to generate. |
numberOfSubjects |
number of subjects per dataset. |
numberOfTimepoints |
number of observation time points per subject. |
numberOfTreatmentLevels |
number of treatment levels. Default: 1 (no treatment effect, intercept and time only). |
timeRange |
numeric vector of length 2, range of time values (min, max).
Default: |
errorGenerator |
random number generator used for the errors. Called as
|
randomEffectGenerator |
random number generator used for the spherical
random effects. Called as |
trueBeta |
scalar or vector with the true values of the fixed effects coefficients. Can be of length one in which case it will be replicated to the required length. The order is: intercept, time, treatment contrasts (if any), treatment-by-time interactions (if any). |
trueSigma |
scalar with the true value of the error scale. |
trueTheta |
numeric vector of length 3 with the true values for the
Cholesky factor of the random effects covariance matrix (lme4 convention).
Default: |
contamFun |
optional contamination function. If provided, it receives the full dataset (a data frame with columns id, time, treatment, y) and an info list, and must return the (possibly modified) data frame. This allows arbitrary contamination including changing group assignments. See Details for the contents of the info list. |
... |
all additional arguments are added to the returned list. |
The generated data follows the model:
where is the number of treatment groups,
with 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:
The Cholesky factor is:
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
list with generators and the original arguments
generateData: |
function to generate data taking one argument, the dataset index. |
createXMatrix: |
function to generate X matrix taking one
argument, the result of |
createZMatrix: |
function to generate Z matrix taking one
argument, the result of |
createLambdaMatrix: |
function to generate Lambda matrix taking
one argument, the result of |
randomEffects: |
function to return the generated random effects taking one argument, the dataset index. |
sphericalRandomEffects: |
function to return the generated spherical random effects taking one argument, the dataset index. |
errors: |
function to return the generated errors taking one argument, the dataset index. |
allRandomEffects: |
function without arguments that returns the matrix of all generated random effects. |
allErrors: |
function without arguments that returns the matrix of all generated errors. |
numberOfDatasets: |
|
numberOfSubjects: |
|
numberOfTimepoints: |
|
numberOfTreatmentLevels: |
|
numberOfRows: |
number of rows in the generated dataset |
trueBeta: |
true values used for beta |
trueSigma: |
true value used for sigma |
trueTheta: |
true values used for theta |
formula: |
formula to fit the model using |
...: |
additional arguments passed via |
Manuel Koller
generateAnovaDatasets, generateMixedEffectDatasets
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))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))
Generates mixed effects datasets using parametric bootstrap.
generateMixedEffectDatasets( numberOfDatasetsToGenerate, preparedDataset, errorGenerator = rnorm, randomEffectGenerator = rnorm )generateMixedEffectDatasets( numberOfDatasetsToGenerate, preparedDataset, errorGenerator = rnorm, randomEffectGenerator = rnorm )
numberOfDatasetsToGenerate |
number of datasets to generate. |
preparedDataset |
dataset as prepared by
|
errorGenerator |
random number generator used for the errors. |
randomEffectGenerator |
random number generator used for the spherical random effects. |
list with generators and the contents of the prepared dataset. See
prepareMixedEffectDataset and
generateAnovaDatasets for a description of the contents.
Manuel Koller
generateAnovaDatasets,
prepareMixedEffectDataset and
createDatasetsFromList
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))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 balanced repeated-measures datasets with a structured random-effect covariance.
generateRepeatedMeasuresDatasets( numberOfDatasetsToGenerate, numberOfSubjects, numberOfVisits, numberOfReplicates = 1L, structure = c("unstructured", "cs", "ar1", "diag"), marginalSd = 1, correlation = 0, trueBeta = 0, trueSigma = 1, errorGenerator = rnorm, randomEffectGenerator = rnorm )generateRepeatedMeasuresDatasets( numberOfDatasetsToGenerate, numberOfSubjects, numberOfVisits, numberOfReplicates = 1L, structure = c("unstructured", "cs", "ar1", "diag"), marginalSd = 1, correlation = 0, trueBeta = 0, trueSigma = 1, errorGenerator = rnorm, randomEffectGenerator = rnorm )
numberOfDatasetsToGenerate |
number of datasets to generate. |
numberOfSubjects |
number of subjects (grouping levels). |
numberOfVisits |
number of levels of the within-subject factor (the dimension of the random effect). |
numberOfReplicates |
number of replicates per subject-by-visit cell. |
structure |
random-effect covariance structure, one of
|
marginalSd |
marginal standard deviation(s) of the random effects,
recycled to length |
correlation |
the structure's single correlation parameter (the
common correlation for |
trueBeta |
the true intercept (the only fixed effect). |
trueSigma |
the true residual standard deviation. |
errorGenerator, randomEffectGenerator
|
functions used to draw the
errors and spherical random effects, see
|
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: (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.
A list with the same structure as the return value of
generateMixedEffectDatasets.
generateMixedEffectDatasets,
generateLongitudinalDatasets
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) }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) }
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.
generateSensitivityCurveDatasets( data, observationsToChange, shifts, scales, center, formula, ... )generateSensitivityCurveDatasets( data, observationsToChange, shifts, scales, center, formula, ... )
data |
dataset to be modified. |
observationsToChange |
index or logical vector indicating which observations should be modified. |
shifts |
vector of shifts that should be applied one by one to each of the modified observations. |
scales |
vector scales that should be used to scale the observations around their original center. |
center |
optional scalar used to define the center from which the observations are scaled from. If missing, the mean of all the changed observations is used. |
formula |
formula to fit the model using |
... |
all additional arguments are added to the returned list. |
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.
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.
oneWay <- generateAnovaDatasets(1, 1, 10, 5) datasets <- generateSensitivityCurveDatasets(oneWay$generateData(1), observationsToChange = 1:5, shifts = -10:10, formula = oneWay$formula) datasets$generateData(1)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”) “components” – in a generalized
sense – from a fitted mixed-effects model, i.e. from an object
of class rlmerMod or merMod.
## 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)## 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)
object |
a fitted mixed-effects model of class
|
name |
a character string specifying the name of the
“component”. Possible values are:
|
... |
potentially further arguments; not here. |
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.
Unspecified, as very much depending on the name.
getCall();
more standard methods for rlmerMod objects, such as ranef,
fixef, vcov, etc.:
see methods(class="rlmerMod")
## 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")))## 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")))
The self-leverage 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 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).
## S3 method for class 'rlmerMod' hatvalues(model, groups = NULL, ...)## S3 method for class 'rlmerMod' hatvalues(model, groups = NULL, ...)
model |
An |
groups |
Leverage aggregation: |
... |
Currently unused. |
Named numeric vector: one entry per observation
(groups = NULL) or per cluster level.
cooks.distance,
rlmer (the design.weights argument bounds
high-leverage design points)
rlmerMod object.Extends implicitIF by adding the -DAS and
-DAS scoring equations to the implicit-derivative linear
system. The Jacobian of the full score wrt is computed via numDeriv::jacobian
with Richardson extrapolation, reusing the converged DAS scales on
the non- columns (they depend on only); the
Jacobian wrt the response is closed-form. Only DASvar and
DAStau methods are supported (DAStau additionally requires
block sizes ).
implicitIF_full(fit, eps = 1e-06, use.cache = TRUE)implicitIF_full(fit, eps = 1e-06, use.cache = TRUE)
fit |
An |
eps |
Numerical step for |
use.cache |
If |
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 ) 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.
A list with components IF_beta (),
IF_u (), IF_sigma (),
IF_theta (), and the model-based delta-method
covariance vcov_model_delta.
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.
## S3 method for class 'rlmerMod' influence(model, do.coef = TRUE, ...)## S3 method for class 'rlmerMod' influence(model, do.coef = TRUE, ...)
model |
An |
do.coef |
Ignored (kept for compatibility with the
|
... |
Currently unused. |
A list with IF_beta (),
IF_sigma (), IF_theta (), and the full implicitIF_full object as full.
implicitIF_full,
cooks.distance
Apply function for all generated datasets.
lapplyDatasets(datasets, FUN, ..., label, POST_FUN, datasetIndices = "all")lapplyDatasets(datasets, FUN, ..., label, POST_FUN, datasetIndices = "all")
datasets |
Datasets list to be used to generate datasets. |
FUN |
the function to be applied to each generated dataset. The
function will be called like |
... |
optional arguments to |
label |
optional parameter, if present, each result is added an
attribute named label with the value of |
POST_FUN |
function to be applied to the result of |
datasetIndices |
optional vector of dataset indices to fit, useful to
try only a few datasets instead of all of them. Use |
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.
Manuel Koller
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)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)
Convenience function that loads the results stored in each of the files and
then calls mergeProcessedFits to merge them.
loadAndMergePartialResults(files)loadAndMergePartialResults(files)
files |
vector of filenames (including paths) of files containing the processed results |
Manuel Koller
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.
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.
makeRobustbasePsi, bisquarePsi.
Construct a redescending Tukey bisquare psi function.
makeBisquarePsi(c = 4.685)makeBisquarePsi(c = 4.685)
c |
tuning cutoff. The default |
Returns a psi_func_rcpp object whose psi, rho, weight and
derivative slots implement Tukey's bisquare (biweight) function
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).
psi_func_rcpp object usable in rho.e, rho.b.
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.
pf <- bisquarePsi pf@psi(c(-6, -3, 0, 3, 6))pf <- bisquarePsi pf@psi(c(-6, -3, 0, 3, 6))
Construct a redescending psi-function from one of robustbase's
compiled Mpsi families.
makeRobustbasePsi( family = c("bisquare", "lqq", "optimal", "hampel", "ggw"), cc = NULL )makeRobustbasePsi( family = c("bisquare", "lqq", "optimal", "hampel", "ggw"), cc = NULL )
family |
robustbase |
cc |
tuning constant (a scalar for |
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 ( 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.
psi_func_rcpp object usable as rho.e (or rho.b,
rho.sigma.e, rho.sigma.b) in rlmer.
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.
lqqPsi, bisquarePsi,
Mpsi, lmrob.control.
pf <- makeRobustbasePsi("lqq") pf@psi(c(-6, -3, 0, 3, 6))pf <- makeRobustbasePsi("lqq") pf@psi(c(-6, -3, 0, 3, 6))
S-type scale psi-function for high-breakdown scale estimation.
makeSScalePsi(c = 1.547645)makeSScalePsi(c = 1.547645)
c |
tuning constant of the bisquare |
Returns a psi_func_rcpp object for the scale arguments
rho.sigma.e and (for one-dimensional random-effect blocks)
rho.sigma.b of rlmer. With it the scale equations
become M-scale equations with a bounded , as in
S-estimation: with
, the Tukey bisquare
scaled to with tuning constant c,
and for standard normal (Fisher consistency).
The breakdown point of the scale equation is : 1/2
for the default c = 1.547645, against 0.295 for the default
psi2propII(smoothPsi) scale and 0.156 for
psi2propII(smoothPsi, k = 2.28) (setting = "RSEa").
It is implemented as a weight function for the existing equations,
with , so the design-adaptive (DAS) standardisation, the consistency
constant and everything built on the scale estimates (standard errors,
Satterthwaite degrees of freedom) carry over unchanged.
A high-breakdown scale does not make the fit high-breakdown by itself:
the fixed effects also need a redescending rho.e (e.g.
bisquarePsi) and, when the start is contaminated, a
high-breakdown start (init = "robust"); see rlmer.
The S-scale is less efficient than the default scale estimate: its
asymptotic efficiency for the scale at the normal is 0.54 with
c = 1.547645; a larger c trades breakdown point for
efficiency (c = 2.937015 gives , breakdown point
1/4 and efficiency 0.87). It is meant for rho.sigma.e and simple
variance components; rlmer refuses it for blocks of correlated
random effects and as rho.e or rho.b. To change c,
use makeSScalePsi(c = ...), not chgDefaults.
a psi_func_rcpp object; attribute "sscale" holds
b, kappa0 and c.
Experimental. The interface and the results may change in a later version.
Rousseeuw PJ, Yohai VJ (1984). Robust regression by means of S-estimators. In Robust and Nonlinear Time Series Analysis, Lecture Notes in Statistics 26, 256–272. Springer.
ss <- makeSScalePsi() attr(ss, "sscale") ## the bounded chi(t) = W(t) (t^2 - kappa0) of the scale equation t <- c(0, 0.5, 1, 2, 5, 50) ss@wgt(t) * (t^2 - attr(ss, "sscale")$kappa0)ss <- makeSScalePsi() attr(ss, "sscale") ## the bounded chi(t) = W(t) (t^2 - kappa0) of the scale equation t <- c(0, 0.5, 1, 2, 5, 50) ss@wgt(t) * (t^2 - attr(ss, "sscale")$kappa0)
Combine list of processed fits into one list in matrix form.
mergeProcessedFits(processedFitList)mergeProcessedFits(processedFitList)
processedFitList |
list of processed fits as produced by
|
similar list as returned by processFit just with
matrix entries instead of vectors.
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)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 miscellaneous utilities for instances of the PsiFunction class.
## 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)## 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)
object |
instance of class |
show(smoothPsi)show(smoothPsi)
Computes a partial moment for the standard normal distribution. This is the
expectation taken not from -Infinity to Infinity but just to z.
partialMoment_standardNormal(z, n)partialMoment_standardNormal(z, n)
z |
partial moment boundary, the expectation is taken from -Inf to z. |
n |
which moment to compute, needs to be >= 2. |
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)
partialMoment_standardNormal(0, 2)partialMoment_standardNormal(0, 2)
The plot method objects of class
PsiFunction simply visualizes the
, , and weight functions and their
derivatives.
## 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", ...)## 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", ...)
x |
instance of class |
y |
(optional) vector of abscissa values (to plot object at). |
which |
|
main |
string or logical indicating the kind of plot title;
either |
col |
colors to be used for the different slots |
leg.loc |
legend placement, see also |
... |
passed to |
If you want to specify your own title, use main=FALSE, and a
subsequent title(...) call.
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(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"))
Diagnostic plots for objects of class rlmerMod and lmerMod.
## 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, ...)## 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, ...)
x |
an object as created by |
y |
currently ignored. |
which |
integer number between 1 and 4 to specify which plot is desired. |
title |
Titles for the different plots. The fourth item can be a format
string passed to |
multiply.weights |
multiply the residuals / random effects with the robustness weights when producing the Q-Q plots. |
add.line |
add reference line to plots, use |
... |
passed on to |
ask |
waits for user input before displaying each plot. |
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").
a list of plots of class ggplot that can be
used for further modification before plotting (using print).
## Not run: rfm <- rlmer(Yield ~ (1|Batch), Dyestuff) plot(rfm) fm <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) plot.rlmerMod(fm) ## End(Not run)## Not run: rfm <- rlmer(Yield ~ (1|Batch), Dyestuff) plot(rfm) fm <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy) plot.rlmerMod(fm) ## End(Not run)
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.
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 )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 )
data |
A data frame containing longitudinal data. Must have columns for subject ID, time, treatment group, and response variable. |
formula |
A formula for the mixed-effects model. Default is
|
idVar |
Character string naming the subject ID column in |
timeVar |
Character string naming the time column in |
treatmentVar |
Character string naming the treatment column in
|
responseVar |
Character string naming the response column in
|
rlmerArgs |
A list of additional arguments passed to
|
lineAlpha |
Numeric in [0, 1]. Transparency of subject lines. Default: 0.6. |
fixedLineWidth |
Numeric. Width of fixed-effect overlay lines. Default: 1.2. |
lowColor |
Color for low robustness weights (potential outliers).
Default: |
highColor |
Color for high robustness weights (typical observations).
Default: |
fixedLineColor |
Color for the fixed-effect prediction lines.
Default: |
fixedLinetype |
Linetype for the fixed-effect prediction lines. Can be
a single value (e.g., |
title |
Optional plot title. |
xlab |
Label for x-axis. If |
ylab |
Label for y-axis. If |
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.
A ggplot object.
rlmer, generateLongitudinalDatasets
## 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)## 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)
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.
## 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, ... )## 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, ... )
object |
An |
newdata, re.form, ReForm, REForm, REform, terms, type, allow.new.levels, na.action, ...
|
See the lme4 |
interval |
One of |
level |
Coverage level for the interval; default 0.95. |
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 (the local-shift
sensitivity). For interval = "prediction" the additional
residual variance is added. Intervals are
fit +/- z * se with ; for
small numbers of clusters the normal approximation may under-cover,
and a bootstrap CI (via confint(..., method = "boot")) is
preferable.
A numeric vector (when interval = "none") or a data
frame with columns fit, lwr, upr, se.
Stable.
This function runs lmer and extracts all information needed to
generate new datasets using parametric bootstrap later.
prepareMixedEffectDataset( formula, data, REML = TRUE, overrideBeta, overrideSigma, overrideTheta, ... )prepareMixedEffectDataset( formula, data, REML = TRUE, overrideBeta, overrideSigma, overrideTheta, ... )
formula |
passed on to |
data |
passed on to |
REML |
passed on to |
overrideBeta |
use to override beta used to simulate new datasets, by
default |
overrideSigma |
use to override sigma used to simulate new datasets, by
default |
overrideTheta |
use to override theta used to simulate new datasets, by
default |
... |
all additional arguments are added to the returned list. |
List that can be passed to
generateMixedEffectDatasets.
data: |
the original dataset |
X: |
the X matrix as returned by
|
Z: |
the Z matrix as returned by
|
Lambda: |
the Lambda matrix as returned
by |
numberOfFixedEffects: |
the number of fixed effects coefficients |
numberOfRandomEffects: |
the number of random effects |
numberOfRows: |
number of rows in the generated dataset |
trueBeta: |
true values used for beta |
trueSigma: |
true value used for sigma |
trueTheta: |
true values used for theta |
formula: |
formula to fit the model using |
...: |
additional arguments passed via |
Manuel Koller
preparedDataset <- prepareMixedEffectDataset(Reaction ~ Days + (Days|Subject), sleepstudy) str(preparedDataset)preparedDataset <- prepareMixedEffectDataset(Reaction ~ Days + (Days|Subject), sleepstudy) str(preparedDataset)
Convenience function to run simulation study in parallel on a single machine.
processDatasetsInParallel( datasets, path, baseFilename, fittingFunctions, chunkSize, saveFitted = FALSE, checkProcessed = FALSE, createMinimalSaveFile = FALSE, ncores = 1, clusterType = "PSOCK", ... )processDatasetsInParallel( datasets, path, baseFilename, fittingFunctions, chunkSize, saveFitted = FALSE, checkProcessed = FALSE, createMinimalSaveFile = FALSE, ncores = 1, clusterType = "PSOCK", ... )
datasets |
dataset list generated by one of the generate functions. |
path |
path to save the datasets to. |
baseFilename |
filename to use, without extension. |
fittingFunctions |
vector of |
chunkSize |
number of datasets to process together in a single job. |
saveFitted |
logical, if true, the raw fits are also stored. |
checkProcessed |
logical, if true, will check whether the contents of the processed output is reproduced for the first dataset. This is useful to ensure that everything is still working as expected without having to re-run the whole simulation study. |
createMinimalSaveFile |
logical, if true, will create a file with the processed results of the first three datasets. This is helpful if one wants to store only the final aggregated results but still wants to make sure that the full code works as expected. |
ncores |
number of cores to use in processing, if set to 1, datasets
are processed in the current R session. Use
|
clusterType |
type of cluster to be created, passed to
|
... |
passed on to |
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.
The list of all processed results merged together.
To help reproducibility, the output of toLatex(sessionInfo(),
locale = FALSE) is stored in the sessionInfo attribute.
Manuel Koller
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.
processFile( file, fittingFunctions, saveFitted = FALSE, checkProcessed = FALSE, createMinimalSaveFile = FALSE, datasets, ... )processFile( file, fittingFunctions, saveFitted = FALSE, checkProcessed = FALSE, createMinimalSaveFile = FALSE, datasets, ... )
file |
file saved by |
fittingFunctions |
vector of |
saveFitted |
logical, if true, the raw fits are also stored. |
checkProcessed |
logical, if true, will check whether the contents of the processed output is reproduced for the first dataset. This is useful to ensure that everything is still working as expected without having to re-run the whole simulation study. |
createMinimalSaveFile |
logical, if true, will create a file with the processed results of the first three datasets. This is helpful if one wants to store only the final aggregated results but still wants to make sure that the full code works as expected. |
datasets |
optional, datasets as stored in |
... |
passed on to |
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”.
The list of all processed results merged together.
To help reproducibility, the output of toLatex(sessionInfo(),
locale = FALSE) is stored in the sessionInfo attribute.
Manuel Koller
Methods to process fitted objects and convert into a data structure that is useful in post-processing.
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, ... )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, ... )
obj |
object returned by the fitting method. |
all |
logical, shorthand to enable all exports. |
coefs |
logical, if true coefficients are added to export. |
stdErrors |
logical, if true, standard errors are added to export. |
tValues |
logical, if true, t-values are added to export. |
sigma |
logical, if true, sigma is added to export. |
thetas |
logical, if true, thetas are added to export. |
b |
scalar logical or index vector, if true, all random effects are
added to export. If an index vector is given, then only the corresponding
random effects are added to the export. The same order as in |
meanB |
logical, if true, the mean of the random effects is added to the export. |
meanAbsB |
logical, if true, the mean of the absolute value of the random effects is added to the export. |
residuals |
scalar logical or index vector, similar to argument
|
converged |
logical, if true, convergence code is added to export. |
numWarnings |
logical, if true, the number of warnings generated during the fitting process is added to export. |
procTime |
logical, if true, time needed to fit object is added to export. |
... |
optional parameters used for some implementations. |
isInterceptCorrelationSlopeModel |
optional logical, can be used to override the assumption that a model with three variance components can be interpreted as having intercept, correlation and slope. |
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.
List with extracted values, most items can be suppressed to save disk space.
label: |
Name of fitting method used to create the fit |
datasetIndex: |
Index of the dataset in the dataset list |
coefficients: |
Vector of estimated fixed-effects coefficients of the fitted model |
standardErrors: |
Vector of estimated standard errors of the fixed-effects coefficients |
tValues: |
Vector of t-Values (or z-Values depending on fitting method) of the fixed-effects coefficients |
sigma: |
Estimated residual standard error |
thetas: |
Vector of random-effects parameter estimates. As parameterized as by
|
b: |
Vector of requested predicted random-effects. |
meanB: |
Vector of means of the predicted random-effects. |
meanAbsB: |
Vector of means of the absolute values of the predicted random-effects. |
residuals: |
Vector of requested residuals. |
converged: |
Convergence status as reported by the fitting method. |
numberOfWarnings: |
the number of warnings generated during the fitting process. |
proc.time: |
Vector of times (user, system, elapsed) as reported by |
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)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)
-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.
## see examples## see examples
The “classical” -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 -function huberPsi is identical to
the one in the package robustbase. The psi slot equals the
identity function within (where is the tuning
parameter). Outside this interval it is equal to . The
rho slot equals the quadratic function within and a
linear function outside.
The smoothed Huber -function is very similar to the
regular Huber -function. Instead of a sharp bend like the
Huber function, the smoothed Huber function bends smoothly. The first tuning
constant, k, can be compared to the tuning constant of the original Huber
function. The second tuning constant, s, determines the smoothness of the
bend.
chgDefaults and psi2propII for changing
tuning parameters; psi_func-class for a more detailed
description of the slots;
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")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")
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.
psi2propII(object, ..., adjust = FALSE) ## S4 method for signature 'psi_func_rcpp' psi2propII(object, ..., adjust = FALSE)psi2propII(object, ..., adjust = FALSE) ## S4 method for signature 'psi_func_rcpp' psi2propII(object, ..., adjust = FALSE)
object |
instance of Rcpp_PsiFunction class to convert |
... |
optional, new default arguments passed to chgDefaults. |
adjust |
logical, whether tuning parameters should be adjusted automatically, such that the scale estimate has the same asymptotic efficiency as the location estimate. |
par(mfrow=c(2,1)) plot(smoothPsi) plot(psi2propII(smoothPsi))par(mfrow=c(2,1)) plot(smoothPsi) plot(psi2propII(smoothPsi))
Half-weight (basin) radius for a fixed-effects start.
ransac_basin_radius(object, cc = NULL, rho = NULL)ransac_basin_radius(object, cc = NULL, rho = NULL)
object |
a fitted |
cc |
without |
rho |
an optional |
Computes the half-weight radius
of the ball around the initial fixed-effects estimate on which the
population Hessian of a redescending M-estimator stays positive
definite (the RANSAC-RSE basin theorem, Step 1). Write
, , for the expected
redescender weight of an observation whose residual carries a
deterministic mean-shift of error-scale units; is
even in and decreasing in near the origin. The
half-weight shift is the first with
. By Cauchy-Schwarz every observation's shift
within the ball is at most , so its expected weight is at least and
throughout the ball. For the
bisquare at the default cutoff , ; for sleepstudy this gives
in the units of .
Up to version 3.5.0-2 this function returned the larger
support-preservation radius
(6.6 for sleepstudy). That radius only kept every observation
inside the redescender's support; it did not bound the Hessian away
from singularity, which is what the basin guarantee needs, so it was
replaced by the half-weight radius.
When rho is supplied, is computed from its
Dpsi slot by numerical integration, so the radius applies to
any redescender (e.g. lqqPsi); cc is then only
the upper end of the search interval for (the rejection
point of rho, found numerically, when cc is
NULL; 10 if rho has none). Without rho,
the bisquare with cutoff cc (default 4.685)
is used. If is not monotone on the search interval, the first
crossing of from the origin is taken.
the half-weight radius in the units of
(a scalar), with attributes "max_xnorm",
"s_half" (the half-weight shift ) and
"cc".
Experimental. The interface and the results may change in a later version.
lm0 <- lme4::lmer(Reaction ~ Days + (Days | Subject), lme4::sleepstudy) ransac_basin_radius(lm0) # bisquare, c = 4.685 ransac_basin_radius(lm0, rho = lqqPsi) # from lqq's psi'lm0 <- lme4::lmer(Reaction ~ Days + (Days | Subject), lme4::sleepstudy) ransac_basin_radius(lm0) # bisquare, c = 4.685 ransac_basin_radius(lm0, rho = lqqPsi) # from lqq's psi'
RANSAC-style random subsample initial estimator for linear mixed-effects models.
ransac_lme4( formula, data, K = 200L, sub_frac = 0.5, scale_fn = NULL, adaptive = TRUE, patience = 50L, K_min = 50L, tol = 0.001, stratify = TRUE, n_keep = 1L, seed = NULL, verbose = FALSE )ransac_lme4( formula, data, K = 200L, sub_frac = 0.5, scale_fn = NULL, adaptive = TRUE, patience = 50L, K_min = 50L, tol = 0.001, stratify = TRUE, n_keep = 1L, seed = NULL, verbose = FALSE )
formula |
model formula in |
data |
full data frame. |
K |
maximum number of random subsamples (default 200). With
|
sub_frac |
fraction of the data per subsample (default 0.5). |
scale_fn |
function from a numeric residual vector to a
scalar scale, used to rank the candidate fits. The default
( |
adaptive |
logical (default |
patience |
number of consecutive non-improving draws that triggers the adaptive stop (default 50). |
K_min |
minimum draws before the adaptive stop can fire (default 50). |
tol |
relative-improvement threshold for the adaptive stop (default 1e-3). |
stratify |
logical (default |
n_keep |
number of distinct best-scoring candidate starts to
return in |
seed |
optional RNG seed for reproducibility. |
verbose |
logical; print progress every 50 subsamples. |
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 residuals are conditional: clusters in the subsample keep their
predicted random effects, the others get none. A contaminated cluster
inside the subsample is therefore absorbed by its own random effect, and
under cluster-level contamination the start's between-cluster
coefficients can be poor (see the section “Ranking RANSAC
candidates” of vignette("simulationStudies")).
The breakdown point any start can give the fixed effects is limited
by the design: with categorical predictors many observations lie on one
hyperplane through the origin, and the ceiling drops (to about 1/4 for
a balanced two-level factor with an intercept); see
design_breakdown.
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.
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).
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 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 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.
Experimental. The interface and the results may change in a later version.
Koller, M. and Stahel, W. A. (2017) Nonsingular subsampling for regression S estimators with categorical predictors. Computational Statistics 32(2), 631–646.
set.seed(1) res <- ransac_lme4(Reaction ~ Days + (Days | Subject), data = sleepstudy, K = 30) res$scaleset.seed(1) res <- ransac_lme4(Reaction ~ Days + (Days | Subject), data = sleepstudy, K = 30) res$scale
Fits the model of object again with the response replaced by
newresp, keeping the psi-functions, the method, the prior weights
and the offset stored in object; the original data is not needed.
The fit is started from the estimates of object. This is the refit
generic of lme4; together with
simulate it gives a bootstrap with
robust refits.
## S3 method for class 'rlmerMod' refit(object, newresp, ...)## S3 method for class 'rlmerMod' refit(object, newresp, ...)
object |
an |
newresp |
the new response: a numeric vector with one value per
observation used in the fit, or a one-column data frame such as a
column of |
... |
not used. |
An rlmerMod object. Its call is that of object.
Stable.
The generics package has a generic of the same name. If it
is attached after lme4, call lme4::refit.
fm <- rlmer(Reaction ~ Days + (1 | Subject), sleepstudy) ystar <- simulate(fm, seed = 1) fixef(refit(fm, ystar$sim_1))fm <- rlmer(Reaction ~ Days + (1 | Subject), sleepstudy) ystar <- simulate(fm, seed = 1) fixef(refit(fm, ystar$sim_1))
Residuals of a fitted rlmer model: the difference of the
observations and the fitted values including the random effects, or a
weighted or whitened version of it.
## S3 method for class 'rlmerMod' residuals( object, type = c("response", "pearson", "weighted", "normalized"), scaled = FALSE, ... )## S3 method for class 'rlmerMod' residuals( object, type = c("response", "pearson", "weighted", "normalized"), scaled = FALSE, ... )
object |
rlmerMod object |
type |
type of residuals, see Details. |
scaled |
divide the |
... |
ignored |
"response", as for
lmer.
"pearson"the response residuals multiplied by the
square root of the prior weights and divided by
. These are the “scaled residuals” of
summary.
"weighted"the response residuals multiplied by the
square root of the prior weights and by the robustness weights
getME(object, "w_e").
"normalized"the marginal residuals (minus the offset) multiplied by the inverse
of the lower Cholesky factor of the fitted marginal covariance
, with the prior weights. At the
model they are uncorrelated with unit variance, so they can be used
to check the fitted within-cluster covariance structure. The
Cholesky factor is taken in the order of the observations in the
data, so each residual is adjusted for the earlier observations of
its cluster.
Stable.
## Not run: fm <- rlmer(Yield ~ (1|Batch), Dyestuff) stopifnot(all.equal(resid(fm, type="weighted"), resid(fm) * getME(fm, "w_e"))) ## End(Not run)## Not run: fm <- rlmer(Yield ~ (1|Batch), Dyestuff) stopifnot(all.equal(resid(fm, type="weighted"), resid(fm) * getME(fm, "w_e"))) ## End(Not run)
Robust estimation of linear mixed effects models, for hierarchical nested and non-nested, e.g., crossed, datasets.
rlmer( formula, data, ..., method = c("DAStau", "DASvar"), setting, rho.e, rho.b, rho.sigma.e, rho.sigma.b, rel.tol = 1e-08, max.iter = 40 * (r + 1)^2, verbose = 0, doFit = TRUE, init, size_obr = FALSE, design.weights = NULL, control = rlmerControl() ) lmerNoFit(formula, data = NULL, ..., initTheta)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, control = rlmerControl() ) lmerNoFit(formula, data = NULL, ..., initTheta)
formula |
a two-sided linear formula object describing the
fixed-effects part of the model, with the response on the left of a
|
data |
an optional data frame containing the variables named in
|
... |
Additional parameters passed to lmer to find the initial
estimates. See |
method |
method to be used for estimation of theta and sigma, see Details. |
setting |
a string specifying suggested choices for the arguments
|
rho.e |
object of class psi_func, specifying the functions to use for the huberization of the residuals. |
rho.b |
object of class psi_func or list of such objects (see Details), specifying the functions to use for the huberization of the random effects. |
rho.sigma.e |
object of class psi_func, specifying the weight functions
to use for the huberization of the residuals when estimating the variance
components, use the |
rho.sigma.b |
(optional) object of class psi_func or list of such
objects, specifying the weight functions to use for the huberization of
the random effects when estimating the variance components (see Details).
Use |
rel.tol |
relative tolerance used as criteria in the fitting process. |
max.iter |
maximum number of iterations allowed. |
verbose |
verbosity of output. Ranges from 0 (none) to 3 (a lot of output) |
doFit |
logical scalar. When |
init |
optional lmerMod- or rlmerMod-object to use for starting values,
a list with elements ‘fixef’, ‘u’, ‘sigma’,
‘theta’, the string
An lmerMod or rlmerMod passed here supplies the model, not only
the starting values: the working object is built from it, so
When |
size_obr |
logical scalar; if |
design.weights |
Mallows-type design weights for robustness to
high-leverage design points: |
control |
an |
initTheta |
parameter to initialize theta with (optional) |
This function implements the Robust Scoring Equations estimator for linear
mixed effect models. It can be used much like the function
lmer in the package lme4. The supported models
are the same as for lmer (gaussian family only). The
robust approach used is based on the robustification of the scoring
equations and an application of the Design Adaptive Scale approach.
Example analyses and theoretical details on the method are available in the
vignette (see vignette("rlmer")).
Models are specified using the formula argument, using the same
syntax as for lmer. Additionally, one also needs to
specify what robust scoring or weight functions are to be used (arguments
starting with rho.). By default a smoothed version of the Huber
function is used. Furthermore, the method argument can be used to
speed up computations at the expense of accuracy of the results.
Currently, there are two different methods available for fitting models. They only differ in how the consistency factors for the Design Adaptive Scale estimates are computed. Available fitting methods for theta and sigma.e:
DAStau (default): For this method, the consistency factors are
computed using numerical quadrature. This is slower but yields more accurate
results. This is the direct analogue to the DAS-estimate in robust linear
regression.
DASvar: This method computes the consistency factors using a
direct approximation which is faster but less accurate. For complex models
with correlated random effects with more than one correlation term, this is
the only method available.
DAStau supports blocks of random effects of dimension at most 2;
for fits containing a larger block it falls back to DASvar with a
warning. An experimental option,
options(robustlmm.dastau.mc = TRUE), lifts this restriction via a
Monte-Carlo calibration of the consistency factors; it is
simulation-validated but not backed by a finite-sample theorem, see
robustlmm-options.
The tuning parameters of the weight functions “rho” can be used to
adjust robustness and efficiency of the resulting estimates (arguments
rho.e, rho.b, rho.sigma.e and rho.sigma.b).
Better robustness will lead to a decrease of the efficiency. With the default
setting, setting = "RSEn", the tuning parameters are set to yield
estimates with approximately 95% efficiency for the fixed effects. The
variance components are estimated with a lower efficiency but better
robustness properties.
One has to use different weight functions and tuning parameters for simple
variance components and for such including correlation parameters. By
default, they are chosen appropriately to the model at hand. However, when
using the rho.sigma.e and rho.sigma.b arguments, it is up to
the user to specify the appropriate function. See
asymptoticEfficiency for methods to find tuning parameters
that yield a given asymptotic efficiency.
For simple variance components and the residual error scale use the
function psi2propII to change the tuning parameters. This is
similar to Proposal 2 in the location-scale problem (i.e., using the
squared robustness weights of the location estimate for the scale estimate;
otherwise the scale estimate is not robust).
For multi-dimensional blocks of random effects modeled, e.g.,
a model with correlated random intercept and slope, (referred to as
block diagonal case below), use the chgDefaults function to
change the tuning parameters. The parameter estimation problem is
multivariate, unlike the case without correlation where the problem was
univariate. For the employed estimator, this amounts to switching from
simple scale estimates to estimating correlation matrices. Therefore
different weight functions have to be used. Squaring of the weights (using
the function psi2propII) is no longer necessary. To yield
estimates with the same efficiency, the tuning parameters for the
block diagonal are larger than for the simple case. Tables of tuning parameters
are given in Table 2 and 3 of the vignette (vignette("rlmer")).
The residual scale is estimated from
,
, with the weights of
rho.sigma.e. For the psi2propII weights of the
default settings is bounded by its limit
, and the equation breaks down () at the fraction of
outliers: 0.295 for "RSEn" and 0.156 for "RSEa". The scale
is inflated well before that, and a redescending rho.e then no
longer rejects the outliers, whatever the start. The simple variance
components behave in the same way at the level of the random effects.
setting = "RSEhb" is a high-breakdown configuration: the S-type
scale sScalePsi (breakdown point 1/2) for
rho.sigma.e and for simple variance components, a bisquare
rho.b for simple random effects, the start init =
"robust", and a bisquare rho.e with c = 3.444 (85%
efficiency). Its fixed effects break down under considerably more
contamination than those of the default settings, at a cost in
efficiency and in the precision of the variance components. For
inference use the cluster sandwich (vcov_type = "sandwich" in
confint and
anova) and the default method =
"DAStau": "DASvar" under-estimates the residual scale in this
setting (and "DAStau" falls back to it, with a warning, for
random-effect blocks of dimension above 2). Blocks of correlated random
effects keep their multivariate rho.b and rho.sigma.b and
get no protection against outlying clusters. See the section
“Studies for the high-breakdown configuration” of
vignette("simulationStudies") for the breakdown, efficiency and
coverage results and their limits.
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 for rho.e and
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 for both rho.b and
rho.sigma.b. Tables of tuning parameter are given in Table 2 and 3 of
the vignette (vignette("rlmer")).
If custom weight functions are specified using the argument rho.b
(rho.e) but the argument rho.sigma.b (rho.sigma.e) is
missing, then the squared weights are used for simple variance components
and the regular weights are used for variance components including
correlation parameters. The same tuning parameters will be used when
setting = "RSEn" is used. To get
higher efficiency either use setting = "RSEa" (and only set arguments
rho.e and rho.b). Or specify the tuning parameters by hand
using the psi2propII and chgDefaults functions.
To specify separate weight functions rho.b and rho.sigma.b for
different variance components, it is possible to pass a list instead of a
psi_func object. The list entries correspond to the groups as shown by
VarCorr(.) when applied to the model fitted with lmer. A set
of correlated random effects count as just one group.
lmerNoFit:The lmerNoFit function can be used to get trivial starting values.
This is mainly used to verify the algorithms to reproduce the fit by
lmer when starting from trivial initial values.
object of class rlmerMod.
weights is passed through to
lmer and has the usual meaning,
. Weighted
method = "DASvar" fits made with robustlmm 3.2-0 to 3.4-5
are wrong and should be refitted: the Design Adaptive Scale
counted the prior weights twice, so , the
variance components and every standard error derived from them
were biased by an amount that depended on the scale the weights
happened to be supplied on. The default method = "DAStau"
was affected only through its starting value, and in practice
only at rho.e = rho.b = cPsi. Unweighted fits are
unaffected in every version. See the 3.5.0 entry of
news(package = "robustlmm") for the magnitudes and the
measured exposure of each method.
rlmer warns if the fit did not converge, and print and
summary then show a “Convergence code” block. The code
is in object@optinfo$conv$opt: 0 converged; 1 the iterations
stopped at max.iter; 2 converged, but the estimating equations
are not satisfied; 3 an inner iteration of the last step did not
converge. The messages are in object@optinfo$warnings.
Stable. Experimental, and subject to change: setting = "RSEhb",
init = "robust" and init = "ransac"; structured covariances
(cs(), ar1()) in the formula; size_obr; and the
Monte-Carlo DAS-tau calibration (rlmerControl).
Manuel Koller, with thanks to Vanda Lourenço for improvements.
lmer, vignette("rlmer")
## 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)## 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)
Fit rlmer with a RANSAC-derived initial estimator.
rlmer_ransac( formula, data, K = 200L, sub_frac = 0.5, n_starts = 1L, phony_threshold = 0.99, seed = NULL, max_tries = 5L, scale_fn = NULL, ... )rlmer_ransac( formula, data, K = 200L, sub_frac = 0.5, n_starts = 1L, phony_threshold = 0.99, seed = NULL, max_tries = 5L, scale_fn = NULL, ... )
formula, data
|
passed to |
K, sub_frac, seed
|
passed to |
n_starts |
number of distinct RANSAC starts for the multi-start consensus (default 1 = single best start, the previous behaviour). |
phony_threshold |
a fit is treated as phony (non-interior) when
its maximum |
max_tries |
maximum number of RANSAC re-seeds when a redescending
|
scale_fn |
passed to |
... |
other arguments to |
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 one with the
lowest residual scale (by default the centred bisquare M-scale about
zero, see scale_fn) whose random-effects covariance is interior
( phony_threshold). This samples several
basins of a redescending and so recovers the interior
solution when the single best start happens to fall into the phony
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 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
(e.g. rho.e = smoothPsi) or a different start.
What the start can change. A start decides which solution
rlmer reaches only where a redescending fit has several, as with
the phony attractor above. It does not raise the
breakdown point of the final fit: once the scale equations of
rlmer break down (at a fraction of 0.295 outliers for the default
settings and 0.156 for setting = "RSEa"), a redescending
rho.e accepts the outliers from any start. See the section
“Ranking RANSAC candidates” of vignette("simulationStudies").
The RANSAC start is also not high-breakdown itself under contamination
(every half-sample contains outliers); rlmer(setting = "RSEhb")
combines high-breakdown scale equations, redescending
-functions and the start init = "robust", see
rlmer.
rlmerMod object.
Experimental. The interface and the results may change in a later version.
fit <- rlmer_ransac(Reaction ~ Days + (Days | Subject), data = sleepstudy, K = 30)fit <- rlmer_ransac(Reaction ~ Days + (Days | Subject), data = sleepstudy, K = 30)
Collects the tuning parameters of rlmer's fitting
algorithm in one object for its argument control.
rlmerControl( rel.tol = 1e-08, max.iter = NULL, dastau.mc = NULL, dasmc.nsim = NULL, dasmc.seed = NULL, retry = FALSE, retry.steps = c("iterations", "DASvar", "robust"), retry.max.iter.factor = 5, lmerControl = NULL )rlmerControl( rel.tol = 1e-08, max.iter = NULL, dastau.mc = NULL, dasmc.nsim = NULL, dasmc.seed = NULL, retry = FALSE, retry.steps = c("iterations", "DASvar", "robust"), retry.max.iter.factor = 5, lmerControl = NULL )
rel.tol |
relative tolerance used as criterion in the fitting
algorithms. Default |
max.iter |
maximum number of iterations. Default |
dastau.mc, dasmc.nsim, dasmc.seed
|
settings of the experimental
Monte-Carlo calibration of the DAS-tau method for this fit. Default
|
retry |
fit again if the fit does not converge? Default
|
retry.steps |
the fallbacks to try, in this order; a subset of
|
retry.max.iter.factor |
factor by which |
lmerControl |
a control object of |
The arguments rel.tol and max.iter of rlmer keep
working. A setting may be given either there or in control, not
in both. To change a setting of control with update,
update control itself.
A list of class "rlmerControl".
With retry = TRUE, a fit that does not converge (see the
section “Convergence” of rlmer) is fitted again
along a chain of fallbacks, and the first converged fit is returned:
"iterations": continue from the unconverged estimate;
"DASvar": fit with method = "DASvar", and start
the requested method from that fit (skipped if the method already
is "DASvar");
"robust": start from the high-breakdown start,
init = "robust".
Every attempt may use retry.max.iter.factor times as many
iterations as the first fit.
Which attempt succeeded is recorded in object@optinfo$retry
and shown in the convergence block of print and
summary. If no attempt converges, the first fit is returned
with its warnings.
Stable. Experimental, and subject to change: the retry
(retry, retry.steps, retry.max.iter.factor) and the
Monte-Carlo DAS-tau settings.
rlmerControl(max.iter = 200, retry = TRUE) ## Not run: rlmer(Reaction ~ Days + (Days | Subject), sleepstudy, control = rlmerControl(retry = TRUE)) ## End(Not run)rlmerControl(max.iter = 200, retry = TRUE) ## Not run: rlmer(Reaction ~ Days + (Days | Subject), sleepstudy, control = rlmerControl(retry = TRUE)) ## End(Not run)
Class "rlmerMod" of Robustly Fitted Mixed-Effect Models
A robust mixed-effects model as returned by rlmer.
Stable, including the emmeans support.
Objects are created by calls to
rlmer.
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:
deviance (disabled, see below)
extractAIC (disabled, see below)
logLik (disabled, see below)
ranef (only partially implemented)
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.
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.
rlmer; corresponding class in package lme4:
merMod
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)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)
Methods for the tidy, glance and augment generics
of the generics package (re-exported by broom and
broom.mixed) for models fitted by rlmer. The column
names follow broom.mixed. The methods are available once one of
these packages is loaded. They return a tibble if the tibble
package is installed and a data frame otherwise.
tidy.rlmerMod( x, effects = c("fixed", "ran_pars"), conf.int = FALSE, conf.level = 0.95, vcov_type = c("default", "sandwich"), df = c("auto", "none", "satterthwaite"), cluster = NULL, ... ) glance.rlmerMod(x, threshold = 0.5, ...) augment.rlmerMod( x, data = model.frame(x), newdata = NULL, se_fit = FALSE, influence = FALSE, ... )tidy.rlmerMod( x, effects = c("fixed", "ran_pars"), conf.int = FALSE, conf.level = 0.95, vcov_type = c("default", "sandwich"), df = c("auto", "none", "satterthwaite"), cluster = NULL, ... ) glance.rlmerMod(x, threshold = 0.5, ...) augment.rlmerMod( x, data = model.frame(x), newdata = NULL, se_fit = FALSE, influence = FALSE, ... )
x |
an |
effects |
which rows to return: |
conf.int |
add confidence limits for the fixed effects
( |
conf.level |
their confidence level. |
vcov_type |
covariance of the fixed effects: |
df |
for |
cluster |
for |
... |
|
threshold |
robustness weights below this value count as
downweighted in |
data |
the data to augment; by default the model frame. The original data can be given if it has one row per observation used in the fit, or additionally the rows dropped for missing values. |
newdata |
new data to predict for. If given, |
se_fit |
add the standard errors of the fitted values
( |
influence |
add Cook's distances ( |
augment returns the robustness weights of the fit next
to the data: .wt_e is the weight of each observation and
.wt_b the weight of the random effect of its cluster. A weight of
1 means the observation or cluster is treated as in the classical fit;
a weight well below 1 marks it as outlying and downweighted. With
several grouping factors there is one column per factor,
.wt_b.<factor>. If a cluster has several random effects, the
smallest of their weights is reported.
tidyone row per parameter, with columns effect,
group, term, estimate, std.error,
statistic, df, p.value and, if requested,
conf.low and conf.high. Random-effect parameters are
named as in broom.mixed (sd__(Intercept),
cor__(Intercept).x, sd__Observation) and have an
estimate only.
glanceone row: nobs, sigma,
converged (NA for an object that was not fitted),
the number of levels of each grouping factor
(ngrps.<factor>), and the number of observations
(nobs.downweighted) and of clusters per grouping factor
(ngrps.downweighted.<factor>) with a robustness weight below
threshold.
augmentthe data with the columns .fitted,
.resid (response residuals), .wt_e, .wt_b (see
above), .hat (hatvalues)
and, if requested, .se.fit and .cooksd.
glance has no logLik, AIC, BIC or
deviance columns: the robust estimator solves estimating
equations that are not the derivative of a likelihood, so these
quantities are not defined for an rlmerMod.
Stable.
fm <- rlmer(Reaction ~ Days + (Days | Subject), sleepstudy) if (requireNamespace("generics", quietly = TRUE)) { print(generics::tidy(fm, conf.int = TRUE)) print(generics::glance(fm)) ## the least trusted observations and subjects aug <- generics::augment(fm) print(aug[order(aug$.wt_e)[1:5], ]) }fm <- rlmer(Reaction ~ Days + (Days | Subject), sleepstudy) if (requireNamespace("generics", quietly = TRUE)) { print(generics::tidy(fm, conf.int = TRUE)) print(generics::glance(fm)) ## the least trusted observations and subjects aug <- generics::augment(fm) print(aug[order(aug$.wt_e)[1:5], ]) }
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).
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.maxNumeric, default 5000.
The size cutoff for computing the Satterthwaite df under
df = "auto". The cost of the df is a deterministic
dimension-only “workload” – one score
evaluation per parameter-Jacobian column: for method "DASvar" and for
"DAStau", where is the number of observations,
the number of fixed effects, the number of random
effects and the number of variance parameters. If
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 for a single random intercept fit with
"DASvar" ( 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")).
For a single fit these settings can also be given in
rlmerControl, which takes precedence over the options.
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 ;
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 . 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 the
classical quadrature path remains the default and the Monte-Carlo
path offers no improvement there.
robustlmm.dastau.mcLogical, default FALSE.
Master switch. When TRUE, method = "DAStau" uses
the Monte-Carlo calibration for all non-diagonal random-effect
blocks of dimension (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.allLogical, default
FALSE. Research switch: also use the Monte-Carlo path
for blocks of dimension , replacing the Gauss-Hermite
quadrature. Intended only for comparing the two calibration
paths; it offers no improvement over the quadrature.
robustlmm.dasmc.nsimInteger, 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.seedInteger, 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.
Some computations can take minutes on large fits. Before starting one,
robustlmm compares a dimension-only measure of its cost with a
threshold and, above it, prints a message that names the cost and
how to avoid it. Nothing is skipped and no question is asked. The
cost is measured in the workload units of
robustlmm.summary.df.max above.
robustlmm.message.if.maxNumeric, default
50000. Threshold on for computing the influence
function of a fit (implicitIF_full), which the
Satterthwaite degrees of freedom, the cluster sandwich and
cooks.distance need. It is computed once and cached.
robustlmm.message.boot.maxNumeric, default
1e7. Threshold on nsim times the number of robust
refits per replicate times , for anova(test = "boot")
(two refits), anova(test = "score") (one) and
confint(method = "boot") (counted as 0.02, because
confintROB refits with lmer).
Set an option to Inf to silence its message, or wrap the call
in suppressMessages.
robustlmm.check_rhs_optimisationLogical, default
FALSE. When TRUE, rlmer cross-checks
the vectorised right-hand-side computation in the block-diagonal
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.
## show the df on larger fits (raise the size cutoff) ## Not run: options(robustlmm.summary.df.max = 20000) ## End(Not run)## show the df on larger fits (raise the size cutoff) ## Not run: options(robustlmm.summary.df.max = 20000) ## End(Not run)
Saves dataset to one or more files.
saveDatasets(datasets, path = getwd(), file, chunkSize)saveDatasets(datasets, path = getwd(), file, chunkSize)
datasets |
dataset list generated by one of the generate functions. |
path |
path to save the datasets to. |
file |
filename to use, without extension. |
chunkSize |
if provided, datasets are split into |
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.
filename or vector of filenames.
Manuel Koller
Shorten labels created by the various fitDatasets functions,
for use in plotting, etc.
shortenLabelsKS2022(labels)shortenLabelsKS2022(labels)
labels |
vector of labels as assigned by |
The labels are shortened as they are in the simulation study published in Koller and Stahel (2022).
Vector of shortened labels
Manuel Koller
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.
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)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)
Generates new responses for the design of a model fitted by
rlmer: parametric samples from the fitted model, or wild
bootstrap samples built from its residuals. These are the samples of
the parametric and wild bootstraps of confintROB; use them to run
a bootstrap of your own.
## S3 method for class 'rlmerMod' simulate( object, nsim = 1, seed = NULL, type = c("parametric", "wild"), use.u = FALSE, cluster = NULL, wild.weights = c("mammen", "rademacher"), ... )## S3 method for class 'rlmerMod' simulate( object, nsim = 1, seed = NULL, type = c("parametric", "wild"), use.u = FALSE, cluster = NULL, wild.weights = c("mammen", "rademacher"), ... )
object |
an |
nsim |
number of responses to simulate. |
seed |
|
type |
|
use.u |
for |
cluster |
for |
wild.weights |
for |
... |
not used. |
type = "parametric"draws from the central Gaussian model
at the robust estimates: random effects with the fitted covariance
and errors with variance for
prior weights , added to the fitted fixed part (including any
offset). The samples follow the model the robust fit describes for
the bulk of the data; they do not reproduce contamination, so
they contain no outliers. With use.u = TRUE the fitted
random effects are kept and only the errors are redrawn. This leaves
out the between-cluster variation and is not suitable for confidence
intervals of the model parameters.
type = "wild"the cluster wild bootstrap of Modugno and
Giannerini (2015). The marginal residuals are rescaled within each cluster for the
leverage of the fixed-effects design, multiplied by one random
two-point weight per cluster, and added back to the fixed part. It
makes no assumption on the distribution or the covariance of the
response within a cluster.
A data frame with nsim columns sim_1, ..., each
a simulated response for the observations used in the fit. The seed
is stored in attr(, "seed"), as by
simulate for lmer fits.
All observations of a cluster get the same random weight, so the
clusters have to be independent of each other. With one grouping
factor these are its levels. For nested grouping factors the default
is the outermost (coarsest) factor. Crossed grouping factors have no
independent clusters; cluster must then be given, and the
dependence across its levels is ignored.
Stable.
Mammen E (1993). Bootstrap and wild bootstrap for high dimensional linear models. The Annals of Statistics 21(1), 255–285.
Mason F, Cantoni E, Ghisletta P (2021). Parametric and semi-parametric bootstrap-based confidence intervals for robust linear mixed models. Methodology 17(4), 271–295.
Modugno L, Giannerini S (2015). The wild bootstrap for multilevel models. Communications in Statistics – Theory and Methods 44(22), 4812–4825.
confint with method =
"boot", which runs the bootstrap of confintROB.
fm <- rlmer(Reaction ~ Days + (Days | Subject), sleepstudy) ysim <- simulate(fm, nsim = 3, seed = 1) head(ysim) ## Not run: ## a wild bootstrap of the slope by hand: refit each sample. ## confintROB refits with lmer(), which is much faster: ystar <- simulate(fm, nsim = 200, seed = 2, type = "wild") slopes <- sapply(ystar, function(y) fixef(lmer(y ~ Days + (Days | Subject), sleepstudy))["Days"]) quantile(slopes, c(0.025, 0.975)) ## the same with robust refits slopes <- sapply(ystar[1:20], function(y) { d <- sleepstudy; d$Reaction <- y fixef(update(fm, data = d))["Days"] }) ## End(Not run)fm <- rlmer(Reaction ~ Days + (Days | Subject), sleepstudy) ysim <- simulate(fm, nsim = 3, seed = 1) head(ysim) ## Not run: ## a wild bootstrap of the slope by hand: refit each sample. ## confintROB refits with lmer(), which is much faster: ystar <- simulate(fm, nsim = 200, seed = 2, type = "wild") slopes <- sapply(ystar, function(y) fixef(lmer(y ~ Days + (Days | Subject), sleepstudy))["Days"]) quantile(slopes, c(0.025, 0.975)) ## the same with robust refits slopes <- sapply(ystar[1:20], function(y) { d <- sleepstudy; d$Reaction <- y fixef(update(fm, data = d))["Days"] }) ## End(Not run)
Method that splits up dataset objects into smaller chunks, so that they can be processed separately.
splitDatasets(datasets, chunkSize = 50)splitDatasets(datasets, chunkSize = 50)
datasets |
dataset object to split into chunks |
chunkSize |
number of datasets to keep in one chunk |
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:
chunkIndex: |
index of the chunk |
Manuel Koller
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)))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)))
Pre-built psi_func_rcpp from makeSScalePsi with
c = 1.547645 (breakdown point 1/2), for rho.sigma.e and
one-dimensional rho.sigma.b.
Experimental. The interface and the results may change in a later version.
rlmerMod object.Computes the robust score sandwich , where is the Schur-complement
(marginal) Jacobian of the profiled -score and sums the per-cluster -score
contributions
( the Mallows design weights, 1 unless
design.weights is used). Equal to the user-facing
vcov(object, type = "sandwich").
vcov_sandwich(fit, cluster = NULL, correction = c("G1", "none", "CR2", "CR3"))vcov_sandwich(fit, cluster = NULL, correction = c("G1", "none", "CR2", "CR3"))
fit |
|
cluster |
Cluster specification; see |
correction |
One of |
Exact when every grouping factor is nested within the clusters (a
single grouping factor, or cluster naming the top level of a
nested design); approximate for crossed factors, with a warning.
Small-sample corrections. "G1" (default) scales by
. This removes the small-sample bias for the intercept but
not for between-cluster coefficients, for which the sandwich stays
biased down when there are few clusters. "CR2" is the
bias-reduced sandwich of Bell and McCaffrey (2002) adapted to the robust
estimating equations: with split into
cluster contributions, each is replaced by
. For a Gaussian fit it is
unbiased under the fitted model. It takes the symmetric square root
after whitening by the clusters' marginal covariances, whereas the CR2
of Pustejovsky and Tipton (2018) uses their Cholesky factors; the two
coincide when these covariances are multiples of the identity and
otherwise differ slightly. "CR3" uses
: times the one-step
leave-one-cluster-out jackknife. CR2 and CR3 need every random effect
nested within one
cluster. anova and
confint use CR2 with Bell–McCaffrey /
approximate Hotelling degrees of freedom (Pustejovsky and Tipton 2018)
for vcov_type = "sandwich".
and are held fixed (partial
sandwich).
A covariance matrix for ,
with dimnames from the fixed-effect coefficient names and attribute
"n.clusters".
Stable.
Bell RM, McCaffrey DF (2002). Bias reduction in standard errors for linear regression with multi-stage samples. Survey Methodology 28(2), 169–181.
Pustejovsky JE, Tipton E (2018). Small-sample methods for cluster-robust variance estimation and hypothesis testing in fixed effects models. Journal of Business & Economic Statistics 36(4), 672–683.
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.
## S3 method for class 'rlmerMod' vcov( object, type = c("default", "sandwich"), cluster = NULL, correction = c("G1", "none", "CR2", "CR3"), ... )## S3 method for class 'rlmerMod' vcov( object, type = c("default", "sandwich"), cluster = NULL, correction = c("G1", "none", "CR2", "CR3"), ... )
object |
An |
type |
|
cluster |
When |
correction |
When |
... |
Additional arguments passed to the default |
The cluster sandwich is exact when every grouping factor is nested
within the clusters and approximate for crossed factors (with a
warning). With few clusters the default G1 sandwich is biased down for
between-cluster coefficients; correction = "CR2" is the
bias-reduced version, which anova and
confint use for vcov_type =
"sandwich" (see vcov_sandwich). vcov_sandwich
warns for with "G1" or "none".
A covariance matrix for .
Stable.
This is a convenience function to make it simple to access the simulation study script files that are shipped with robustlmm.
viewCopyOfSimulationStudy( study = c("sensitivityCurves.R", "consistencyAndEfficiencyDiagonal.R", "consistencyAndEfficiencyBlockDiagonal.R", "breakdown.R", "breakdownMC.R", "convergence.R", "robustnessDiagonal.R", "robustnessBlockDiagonal.R", "robustnessBlockDiagonalExtended.R"), destinationPath = getwd(), overwrite = FALSE )viewCopyOfSimulationStudy( study = c("sensitivityCurves.R", "consistencyAndEfficiencyDiagonal.R", "consistencyAndEfficiencyBlockDiagonal.R", "breakdown.R", "breakdownMC.R", "convergence.R", "robustnessDiagonal.R", "robustnessBlockDiagonal.R", "robustnessBlockDiagonalExtended.R"), destinationPath = getwd(), overwrite = FALSE )
study |
Name of the script file, partial matching is supported via
|
destinationPath |
optional path to directory in which the copy of the script should be created. By default the current working directory is used. |
overwrite |
logical; should existing destination files be overwritten? |
The function creates a copy of the script file that can be safely edited without changing the original file.
## Not run: viewCopyOfSimulationStudy("sensitivityCurves") ## End(Not run)## Not run: viewCopyOfSimulationStudy("sensitivityCurves") ## End(Not run)