vcov_sandwich() /
vcov(type = "sandwich") gain
correction = "CR2" and "CR3": the bias-reduced
cluster sandwich of Bell and McCaffrey
(2002) and its approximate-jackknife version, adapted to the robust estimating equations. With few clusters the default G1 sandwich is biased down for between-cluster coefficients.
anova() and confint() with
vcov_type = "sandwich" now use the CR2 sandwich by default,
with the approximate Hotelling F-test of Pustejovsky and Tipton (2018)
and Bell-McCaffrey t quantiles (new ddf / df
values "auto", the default, and "cluster").
ddf = "none" / df = "none" restore the
chi-square table and normal intervals of earlier versions; with the
sandwich, "satterthwaite" now falls back to the new
reference. Both functions gain cluster for models with
several grouping factors. A coefficient that a single cluster informs on
its own gets NA, with a warning. The test is much closer to
its level than the chi-square table but can stay somewhat
anti-conservative for between-cluster terms; see the section
“Small-sample references for the cluster sandwich” of
vignette("simulationStudies").
anova(test = "boot") gains nuisance: by
default ("threshold") a nuisance variance component of the
null model estimated close to zero is set to zero when generating the
bootstrap samples, which keeps the bootstrap consistent when that
component is zero (Guedon, Baey and Kuhn 2024); "raw"
restores the previous generator. The simulation study could not show a
finite-sample benefit; see the section “Near-zero nuisance variance
components in the bootstrap test” of
vignette("simulationStudies").
anova(test = "score") tests several added
independent variance components jointly, with the cone (chi-bar)
statistic of the whitened robust score; see the section “Several added
variance components” of
vignette("simulationStudies").
ransac_lme4() and rlmer_ransac() rank
the candidates by a centred residual scale (the bisquare M-scale about
zero) instead of robustbase::Qn, which cannot tell a
candidate with a shifted intercept from the correct one.
scale_fn = robustbase::Qn restores the previous ranking;
rlmer_ransac() gains scale_fn, which also
scores its multi-start consensus. See the section “Ranking RANSAC
candidates” of vignette("simulationStudies").
New fragility(): the first-order fraction of
contaminating clusters that would drive each variance parameter to its
bound. summary() prints it when it computes the
Satterthwaite degrees of freedom.
New design_breakdown(): the breakdown ceiling of the
fixed-effects design (Mili and Coakley 1996) and the breakdown point an
MM-type refinement is guaranteed to inherit from its start (following
Yohai 1987).
New makeSScalePsi() and the pre-built
sScalePsi: an S-type scale function for
rho.sigma.e and one-dimensional rho.sigma.b.
The scale equations of rlmer() break down at a fraction of
0.295 outliers for the default settings and 0.156 for
setting = "RSEa", and a redescending rho.e
stops rejecting the outliers once the scale is inflated.
sScalePsi turns the same equations into the bounded-rho
M-scale equations of S-estimation (breakdown point 1/2), keeping the
consistency constant, the design-adaptive standardisation and everything
built on the scale estimates. rlmer() refuses it for
random-effect blocks of dimension above 1 and as rho.e or
rho.b.
New rlmer(init = "robust"): a high-breakdown start
for any random-effects structure. An S-estimate of the fixed-effects
part (lmrob(method = "S")) screens the observations and
lmer() on the retained observations supplies the starting
values. It is computed for the model as specified (contrasts, subset,
prior weights and offset) and is reproducible (fixed seed, random number
stream untouched). The start matters only with a redescending
rho.e or rho.b.
New setting = "RSEhb" for rlmer():
sScalePsi for the residual scale and simple variance
components, a bisquare rho.b for simple random effects,
init = "robust", and a bisquare rho.e with
c = 3.444 (85% efficiency). Its fixed effects break down
under considerably more contamination than with the default setting, at
a cost in efficiency and variance-component precision; use
vcov_type = "sandwich" for inference and the default
method = "DAStau". ?rlmer has a new section on
the breakdown point; the simulation results are in the section “Studies
for the high-breakdown configuration” of
vignette("simulationStudies").
New rlmerControl() and
rlmer(..., control = ): one object for
rel.tol, max.iter, the settings of the
experimental Monte-Carlo DAS-tau calibration (so far available only
through options()), and an optional retry. The arguments
rel.tol and max.iter and the options keep
working; giving a setting both as an argument and in
control is an error. An lmerControl() object
passed as control still goes to the classical
start.
rlmerControl(retry = TRUE): a fit that does not
converge is fitted again along a chain of fallbacks (continue with more
iterations; start from a "DASvar" fit; start from
init = "robust"), and the first converged fit is returned.
The attempt that succeeded is recorded in fit@optinfo$retry
and shown by print() and summary(). Off by
default.
vignette("robustlmm-3.5.0") has a new section “A
robust sensitivity analysis for an MMRM”: a repeated-measures trial (the
example data of the mmrm package) analysed with
rlmer() next to the primary MMRM, with the treatment
contrast through emmeans and contrast_test(),
the robustness weights from augment() and cluster-level
Cook’s distances, and a statement of what such an analysis does not show
(the models differ; an unstructured covariance over the visits cannot be
fitted by rlmer(); validity under dropout has not been
studied).
The help pages state the lifecycle of each function in a section
“Lifecycle”: stable, or experimental (interface and results may change).
Experimental are anova(test = "score"),
null = "robust" and nuisance = "threshold";
cs() / ar1(); the Monte-Carlo DAS-tau
calibration; sScalePsi / makeSScalePsi(),
setting = "RSEhb" and init = "robust";
fragility() and design_breakdown(); the RANSAC
functions; and the retry of rlmerControl().
A message now announces computations that can take long: the
influence function of a large fit (needed for Satterthwaite degrees of
freedom, the cluster sandwich and Cook’s distances) and bootstraps with
many robust refits (anova(test = "boot") with the default
nsim = 1000 on most data). The message names the cost and
how to avoid it; nothing is skipped. The thresholds are the new options
robustlmm.message.if.max and
robustlmm.message.boot.max (Inf silences
them); see ?"robustlmm-options".
Earlier and clearer errors for invalid arguments in
rlmer(), vcov(), confint(),
anova(), predict() and
contrast_test(). Combinations in which an argument had no
effect are now errors instead of being ignored:
confint(method = "boot") with
vcov_type = "sandwich" (was a warning), with
df or with cluster; cluster (and
for vcov(), correction) without the sandwich
in vcov(), confint() and anova();
null or nuisance with
anova(test = "score") or an explicit
test = "Wald". Remove the argument to get the previous
result. A max.iter that is not a whole number is an error
as well.
New simulate() method for rlmerMod
objects. type = "parametric" draws from the central
Gaussian model at the robust estimates (random effects redrawn by
default, use.u = TRUE keeps the fitted ones; prior weights
and offsets are honoured); type = "wild" gives the cluster
wild bootstrap samples of Modugno and Giannerini (2015), with Mammen or
Rademacher weights and a cluster argument for nested and
crossed designs. These are the samples of the bootstraps in
confintROB; for a single grouping factor without prior
weights or offset the wild samples are identical to its under the same
seed. confint(method = "boot") is unchanged.
New refit() method (the lme4 generic):
refits a model to a new response with the psi-functions, method, prior
weights and offset stored in the fit (the original data is not needed),
started from the original fit. With simulate() it makes
lme4::bootMer() work for rlmerMod
objects.
New tidy(), glance() and
augment() methods for rlmerMod objects
(generics of the generics package, as re-exported by
broom and broom.mixed, whose column names they
follow). augment() returns the robustness weights next to
the data: .wt_e per observation and .wt_b for
the random effect of its cluster. glance() counts the
observations and clusters downweighted below 0.5; it has no
logLik / AIC / BIC because the
robust estimator is not based on a likelihood. They return a tibble if
tibble is installed. generics and
tibble are in Suggests; robustlmm now requires
R >= 3.6.0 (as lme4 does).
New contrast_test(): Wald test of user-specified
linear contrasts of the fixed effects, L beta = rhs. One
contrast gives a t test, several a joint F test, with the references of
summary(), anova() and the emmeans support:
robust Satterthwaite degrees of freedom for the model-based covariance,
CR2 with Bell-McCaffrey / approximate-Hotelling degrees of freedom for
vcov_type = "sandwich".
residuals() gains type = "pearson"
(divided by sigma and whitened by the prior weights; the “scaled
residuals” of summary()) and
type = "normalized" (marginal residuals whitened by the
Cholesky factor of the fitted marginal covariance), for checking the
fitted within-cluster covariance structure.
anova(test = "boot") and
anova(test = "score") for models with prior weights or an
offset. The bootstrap samples were generated with a common error
variance and without the offset, so they did not come from the fitted
null model; they now use the error variance sigma^2 / w and add the
offset, as simulate() does. The refits also failed for such
models unless the weights variable happened to be visible outside the
data (the test then stopped with “only 0 of … replicates succeeded”);
the bootstrap data now carries the variables of weights and
offset. Models without prior weights and offset are
unaffected: their seeded results are identical.
Documentation: the references to Mason, Cantoni and Ghisletta in
?confint.rlmerMod and the vignettes gave a wrong title for
the 2024 paper and wrong end pages for the 2021 paper; corrected. The
statement that the parametric bootstrap “is exact under the fitted
central model” is replaced by “relies on the fitted central
model”.
residuals(fit) (type = "response") of a
fit with prior weights returned the residuals multiplied by the square
root of the weights. It now returns y - fitted(fit), as
lme4 does. Affected: prior-weighted fits only, including
scaled = TRUE; multiply
residuals(fit, type = "pearson") by sigma(fit)
to get the old values. summary()’s scaled residuals, the
plots and all estimates are unchanged.
rlmer(..., na.action = na.exclude) failed with
“length(v_e) == n is not TRUE” when the data had missing values. It now
fits, and residuals() and predict() pad the
excluded rows with NA.
rlmer() now reports convergence problems. The fit
collected its own warnings and never showed them, so a fit that stopped
at max.iter, or whose estimating equations were not
satisfied at convergence, was returned silently. These warnings are now
raised once after the fit, stored in fit@optinfo$warnings,
and shown in a “Convergence code” block by print() and
summary(); fit@optinfo$conv$opt is 0 for a
converged fit. Warnings of the inner iterations are reported only when
they concern the returned estimate. A fit with max.iter = 0
(the start itself) is reported as not converged. Affected: any fit that
did not converge now warns; converged fits and all estimates are
unchanged. Wrap the call in suppressWarnings() to silence
it; the bootstrap and RANSAC loops of the package do so for their
refits.
The effects solve (fitEffects()) failed with a
singular LU factorisation when, under redescending psi-functions, a
random effect and all of its observations had zero weight. A ridge of
relative size 1e-10 on the random-effect block now sets such random
effects to zero; a fixed-effects design without information still
errors.
rlmer(..., subset = ) warned that init
had been fitted on fewer observations than formula and
data give: the check meant for a user-supplied
lmer() fit also ran on the internal start.
rlmer(init = list(...)) discarded the supplied fixed
effects, and with them the fixed effects of every RANSAC start
(init = "ransac", rlmer_ransac()); this
affects robustlmm 3.5.0-2 as released on 2026-07-30. The supplied values
were applied, and (beta, u) re-solved at the supplied
sigma and theta, on the object coerced from
the internal lmer() start, before the fit’s psi-functions
were set; the re-solve therefore used the classical psi and returned the
GLS fit. The values are now applied to the fully configured fit, so the
re-solve uses rho.e and rho.b. For monotone
psi-functions (the default) the converged fit does not depend on the
start and is unchanged; for a redescending rho.e
(e.g. bisquarePsi) the start now selects the solution, as
documented. Re-runs of the affected RANSAC simulation studies changed no
conclusion. A supplied u now also updates the fitted values
(with doFit = FALSE the residuals were those of the
internal start’s u).
vcov_sandwich() /
vcov(type = "sandwich") and caseweightIF()
omitted the Mallows design weights (design.weights) from
the observations’ score contributions, although the estimating equations
weight them (sum_i eta_i x_i psi(r_i)). For fits with
design weights below 1 the sandwich standard errors were too large, by
nearly a factor of two in an example checked against a
leave-one-cluster-out jackknife. Fits without design weights are
unchanged.
anova(test = "boot") p-values were arbitrary when
the observed quasi-deviance difference sat at its atom at zero, that is
when both fits put the tested variance component on the boundary. The
computed difference is then rounding noise of either sign, and the
p-value compared that noise with the bootstrap replicates’ noise; under
a null tested component this produced most of the test’s rejections in a
simulation. Differences up to 0.01 (genuine positive values are of
order
test = "score" does not depend on the alternative fit.anova(test = "boot") failed (no successful bootstrap
replicate) when the larger model used a variable that the smaller one
does not, such as the covariate of an added random slope
(0 + z | g): the bootstrap data sets were built from the
smaller model’s frame only. They now include every variable either model
uses, also one that enters only transformed
((0 + log(z) | g): the model frame holds
log(z), not z), which is taken from the fit’s
data.
ransac_basin_radius() now returns the
half-weight radius
r* = s_1/2 sigma / max_j ||x_j||, where s_1/2
is the residual mean-shift (in error-scale units) at which the expected
redescender weight m(s) = E psi'(Z + s),
Z ~ N(0, 1), has halved (s_1/2 = 1.410 for the
bisquare at c = 4.685), in place of the
support-preservation radius c sigma / (2 max_j ||x_j||)
returned up to 3.5.0-2. The old radius only kept every observation
inside the redescender’s support; it did not keep the population Hessian
away from singularity, which is what the basin guarantee needs, whereas
on the new ball -H(beta) >= -H_0 / 2 > 0 (the
corrected Step 1 of the RANSAC-RSE basin theorem). On
sleepstudy the radius drops from 6.6 to 4.0, so the
fixed-effect basin check of rlmer_ransac() fires more
often. With rho supplied the shift is computed from its
Dpsi slot, so lqq and the other redescenders get their own
radius (lqq: s_1/2 = 1.423); the return value carries the
new attributes "s_half" and "cc". Re-running
inst/simulationStudy/ransacBasin.R with the new gate (500
replicates x 6 generators) raises the pooled warning rates from
P(warn | phony) = 0.23,
P(warn | not phony) = 0.33 to 0.44 and
0.56: the smaller ball is left about twice as often, but
the check remains a non-discriminating detector of the phony
random-effects correlation (it is a fixed-effect Hessian condition; the
phony failure is caught by the separate |rho-hat| -> 1
check), and the fits themselves are unchanged (phony rates and slope
RMSE identical).
The GitHub Actions workflow gains a test-coverage job (report
uploaded as an artefact; 78% of the R and C++ code is covered by the
test scripts), a check that the help pages are up to date with roxygen2,
a spell check against inst/WORDLIST (seven typos in help
pages and vignettes fixed), and a weekly run against R-devel.
The scripts of the simulation studies for the methods new in
3.5.1 (the small-sample and variance-component studies and the breakdown
study of vignette("simulationStudies")) are now shipped in
the package, in
system.file("simulationStudy", package = "robustlmm").
vignette("simulationStudies") is shipped as a
pre-built PDF (engine R.rsp::asis, R.rsp added
to Suggests) and the slowest chunks of
vignette("robustlmm-3.5.0") show precomputed output, to
keep the check time within the CRAN limit.
calcTau() no longer refuses to solve its fixed point
when the starting value exceeds 1. The tau <= 1
constraint sat in the loop condition, so a start above 1 skipped the
iteration entirely and the post-loop clamp returned exactly 1; it now
starts from min(tau, 1). This is what let the
weighted-tau bug of 3.5.0 reach
method = "DAStau" (see that release’s entry for who should
refit; that part was fixed in 3.5.0, this guard was not). No fit
produced by the current code reaches the condition – the largest
starting value observed over a randomised search is 0.99948 – so the
change is defensive and is a no-op on current behaviour;
tests/tau.R now pins the fixed point to be independent of
the starting value.
init = "ransac" estimated the model on the
RANSAC subsample <<<<<<< HEAD instead of on your
data, in 3.4-5 and 3.5.0-1. If you have used
init = "ransac", rlmer_ransac() or
fitDatasets_rlmer_ransac/_ransac_bisq with
either version, refit: the returned object was an estimate on roughly
half the observations. Nothing else is affected, and no other
init form is affected. ======= instead of on your data, in
3.4-5 and 3.5.0-1 (3.5.0 itself was auto-rejected and never reached
CRAN).** If you have used init = "ransac",
rlmer_ransac() or
fitDatasets_rlmer_ransac/_ransac_bisq with
either version, refit: the returned object was an estimate on roughly
half the observations. Nothing else is affected, and no other
init form is affected. >>>>>>>
master
An lmerMod passed as init supplies the
model, not only the starting values – the working object is
built from it, so formula and data stop
mattering from that point on. That is correct and efficient for a
starting fit on the same data, which is what the argument was designed
for. But ransac_lme4()$fit is fitted on the winning RANSAC
subsample, so resolving init = "ransac" to that object
silently dropped the rest of the data. On a two-dimensional
random-effect block with m = 50 clusters and
n = 1000 observations, the fit used 510 observations and 27
of the 50 clusters.
rlmer() now takes fixef, sigma
and theta from the RANSAC fit and estimates the model on
the full data. u is deliberately not carried
over: it has one entry per random-effect coefficient of the fit it came
from, and zeroing it instead is not a neutral choice – combined with a
theta from elsewhere it drives the random-effect covariance
onto the boundary (theta = 0). Instead the effects are
re-solved at the supplied sigma and theta.
That solve is the robust IRLS step for beta and
u jointly, run to convergence, so the RANSAC
fixef acts as a warm start rather than as a value that is
held fixed; sigma and theta are the parameters
actually held.
Two supporting changes. The components of the list form of
init are now optional; anything omitted keeps the value
from the internal lmer initialisation, which is what makes
the above expressible. And passing an lmerMod fitted on a
different number of observations than data has rows now
warns, since that is the same silent data-loss trap in user code. It is
a warning rather than an error because dropping rows with
NAs is a legitimate reason for the counts to differ.
Four call sites were affected, in three independent copies of the
pattern: init = "ransac", rlmer_ransac(), and
the inline start in fitDatasets_rlmer_ransac and
fitDatasets_rlmer_ransac_bisq (which do not route through
rlmer_ransac()). All four are fixed.
Found while running the estimator against the simulation design of
Sugasawa, Hui and Welsh (2026, JCGS 35(2):720-733). Regression
test: tests/ransac-init-fullwdata.R, which covers all four
entry points.
Several foundations of the new inference toolchain – the cluster
sandwich vcov, Wald and bootstrap confint,
prediction intervals, the observation-level influence diagnostics and
the base RANSAC initial estimator – were first released to CRAN in the
3.4-4 / 3.4-5 series (see below). 3.5.0 adds the features listed here
and extends several of the 3.4-5 foundations.
anova.rlmerMod. Single-fit per-term
robust Wald table; pairwise Wald restriction test for fixed-effects-only
nested models; parametric-bootstrap quasi-deviance test
(test = "boot") for variance-component comparisons on the
PSD-cone boundary, with a common-scale convention (use
fit0’s \(\hat\sigma\)
throughout) to avoid the sign-flip artefact when scales differ between
models. Type-I matches anova(lmer0, lmer1) LRT under all
tested conditions; under error-distribution contamination at moderate
\(J\) rlmer-Wald is 20-30 percentage
points more powerful. The bootstrap path gains an experimental
null = "robust" option: it generates the parametric
bootstrap from a contamination-cleaned null fit (trimming the clusters
the robust fit heavily downweights) instead of fit0’s raw
estimates, which reduces the anti-conservativeness of the
variance-component bootstrap under group contamination at larger \(J\) while preserving power, at the cost of
mild conservatism on clean data. Falls back to the plain parametric null
when no (or more than half the) clusters are flagged, when trimming
would leave the design rank-deficient (a dropped contrast level, a
collapsed grouping factor, or a constant random-slope covariate —
checked with the same nonsingular-subsampling test as the RANSAC initial
estimator), or when the trimmed refit fails. Validated by simulation,
not by a finite-sample theorem.
anova(fit0, fit1, test = "score")
(experimental). One-sided robust score test for the special case of a
single added independent scalar variance component
(e.g. (1|g) vs (1|g) + (0 + x|g), or diagonal
structures adding one component; anything else stops with an error). The
statistic is a self-normalised sum of psi-bounded per-cluster score
contributions computed from the robust null fit only,
calibrated by a score-only parametric bootstrap that refits just the
null model per replicate (~4x cheaper than test = "boot";
default nsim = 199). In simulation (Gaussian balanced
designs, one scalar tested component) it fixes the deviance bootstrap’s
contamination-driven anti-conservativeness without losing power:
contaminated-null Type-I 0.035 vs 0.115 (test = "boot") at
\(J = 50\) with 10% of clusters
shifted, clean-null 0.045, power 0.920 vs 0.900; an adversarial sweep
(small \(J\), uncentered \(x\), observation-level outliers) stayed
within 0.015-0.075 Type-I. No null cleaning is needed
(null = "robust" is ignored) and no downweighted-group
warning is issued for this path. Contamination aligned with the tested
direction is indistinguishable from the alternative for any test with
power; the attached per-cluster contributions
attr(., "boot")$s_j identify the driving clusters (pair
with cooks.distance(., groups = ) /
hatvalues). Simulation-validated, not proven; the deviance
bootstrap test = "boot" remains the default
variance-component path.
confint(., method = "boot" / "BCa") default
boot.type is now "wild" (matching
confintROB’s own recommendation);
boot.type = "parametric" remains available. The bootstrap
delegation itself first shipped in 3.4-5, with "parametric"
as the default.
Cluster-level influence and robust leverage.
cooks.distance(fit, groups = TRUE) (or
groups = "<factor>") gives a per-cluster Cook’s
distance – the Mahalanobis norm of the per-cluster influence function
\(\Xi = -J_{par}^{-1} S\) on \((\hat\beta, \hat\sigma,
\hat\theta)\) – flagging whole-group outliers that no single
observation reveals, as a general influence diagnostic. (It is not a
reliable pre-screen for the bootstrap variance-component test of
anova, which is anti-conservative under group
contamination: the robust fit absorbs a contaminated group into a
downweighted random effect, so the bounded influence saturates exactly
when contamination is present; that test’s automatic guard therefore
screens with the random-effect robustness weights instead.)
hatvalues(fit) returns the robust leverage (the
self-leverage \(A_{ii}\) in the
random-effect-whitened convolution; reduces to the classical
lme4 leverage at rho = cPsi), with
groups to sum it per cluster.
cooks.distance(fit) per observation is unchanged. Cluster
level requires a single or nested grouping structure (crossed designs
error).
RANSAC refinements (extends the
ransac_lme4() / rlmer_ransac() /
init = "ransac" estimator first released in 3.4-5).
ransac_lme4() gains adaptive K (early-stop
once the best Q_n score plateaus; K is now a maximum
budget, adaptive = FALSE restores the fixed loop) and
stratified subsampling (stratify = TRUE, default: each
subsample covers every level of the finest grouping factor, so the
per-subsample lmer stays fittable on designs with many
small clusters). It also gains nonsingular subsampling
(Koller and Stahel 2017, Algorithm 1, the method behind the
setting = "KS2014" default of robustbase):
with categorical predictors a rank-deficient subsample does not always
make lmer error — it silently drops the aliased columns of
a dropped factor level — so, rather than draw-then-repair, each
subsample is now drawn full-rank by construction. A new
C++ routine (nonsingularSubsampleLU) runs a Gaxpy-variant
LU decomposition with partial pivoting and column skipping over a random
permutation of the observations, greedily selecting the first
p linearly-independent observations and skipping any that
are collinear with those already chosen; the whole clusters carrying
that core seed the draw, and further clusters are added until the
retained data are identifiable (full-rank fixed design, ≥ 2 grouping
levels, varying random slopes) and at the target size. The result is
guaranteed identifiable whenever the full design is, so no degenerate
candidate ever enters the scale competition. The new
n_singular return element now counts the collinear
candidate observations skipped by the LU across all draws. A purely
continuous, full-rank design reduces exactly to a uniform random cluster
subsample (the paper’s guarantee). Following Koller and Stahel (2017,
Remark 2), rlmer_ransac() guards the refinement: when a
redescending rho.e zero-weights observations back into a
rank-deficient positive-weight design it automatically re-seeds and
re-fits from a fresh nonsingular RANSAC start, up to the new
max_tries argument (default 5); only if every attempt still
collapses is the last fit returned with a warning (suggesting a
positive-weight rho.e such as smoothPsi). New
exported ransac_basin_radius() returns the
support-preservation radius
r*(c) = c sigma / (2 max_j ||x_j||), where c
is the rejection point of the redescender — the bisquare cutoff, or
found numerically for lqq and the other redescenders.
rlmer_ransac() now runs three post-fit diagnostics: a
fixed-effect basin check (warns when a redescending rho.e
left r* of the high-breakdown start), a random-effects
phony-correlation check (warns, for any rho, when the
fitted RE covariance is near-singular, |rho-hat| -> 1),
and the refinement-singularity check above (Remark 2). A simulation
study (inst/simulationStudy/ransacBasin.R) shows the basin
radius is a fixed-effect condition and the phony failure is a distinct
random-effects attractor, so the two checks are reported separately.
rlmer_ransac(n_starts = m) adds a multi-start consensus: it
fits from the m best distinct RANSAC starts and returns the
lowest-residual- scale fit whose random-effects covariance is interior,
recovering the good solution when the single best start falls into the
phony |rho-hat| -> 1 attractor (per-start summary in
attr(fit, "consensus")); n_starts = 1 is the
previous behaviour.
Redescending psi-functions from robustbase,
and lqqPsi. makeBisquarePsi() /
bisquarePsi now delegate their psi, rho, weight and
derivative evaluations to robustbase’s compiled bisquare
family (Mpsi, Mwgt, Mchi,
normalised via MrhoInf), so they match lmrob()
exactly; the returned values are identical (to numerical tolerance,
verified to 1e-10) to the previous hand-coded implementation. New
exported general constructor makeRobustbasePsi(family, cc)
builds a psi_func_rcpp from any of the
robustbase redescenders — "bisquare",
"lqq", "optimal", "hampel",
"ggw" — defaulting the tuning to
robustbase::.Mpsi.tuning.default(family). The bisquare
redescends comparatively fast; the lqq
(linear-quadratic-quadratic) psi of Koller and Stahel (2011), the
recommended redescender used by
robustbase::lmrob.control(setting = "KS2014"), is now
pre-built and exported as lqqPsi and is recommended over
bisquarePsi for redescending fits (pair it with
init = "ransac" / rlmer_ransac(), which a
redescender needs for a good start).
Satterthwaite degrees of freedom and p-values.
Robust influence-function-based Satterthwaite df, consistent across all
four inference surfaces: summary adds df and
Pr(>|t|) columns;
confint(fit, df = "satterthwaite") uses t-quantiles;
anova(fit, ddf = "satterthwaite") reports an F-test; and
emmeans/emtrends/contrast return
finite df instead of Inf. Unlike lmerTest, the
df derive from the robust IF-based covariance of the variance
parameters, so they stay honest under contamination.
summary now shows the df by default
(df = "auto") when it is cheap to compute — either the
underlying influence function is already cached on the fit, or a
deterministic, dimension-only size workload is within the cutoff
getOption("robustlmm.summary.df.max", 5000) — and otherwise
falls back to the historic t-value table with a one-line note on how to
request it (df = "satterthwaite" forces it,
df = "none" disables it). The cutoff is dimensionless and
machine-independent, so the same fit behaves identically everywhere;
raise robustlmm.summary.df.max to show the df on larger
fits. The expensive influence function is now cached on the fit, so
summary, anova, emmeans,
confint, cooks.distance and the sandwich
vcov share a single computation. Requires the default
vcov. Supports single-factor, nested
((1 | school/class)) and crossed
((1 | subject) + (1 | item)) designs with diagonal random
effects. Nested designs use a one-way cluster sandwich over the coarsest
grouping factor (into which all nested random effects aggregate);
crossed designs use a Cameron-Gelbach-Miller multiway cluster-robust
covariance (V_a + V_b - V_{a:b}), projected to the nearest
positive-semidefinite matrix. In both cases the analytic standard errors
match the empirical ones in Monte-Carlo (ratios ~1). Designs not covered
(e.g. crossed random slopes) fall back to the t-value table. When a
variance component is estimated at 0 (a boundary fit), the df is now
computed conditional on that component being held at the boundary – it
equals the df of the model with that component dropped, and Monte-Carlo
coverage of the resulting intervals is nominal; previously such fits
suppressed the df entirely. Only a genuinely non-identifiable (singular)
fit now suppresses it. The package options are documented under
?\robustlmm-options``.
rlmer(., design.weights = ...).
Mallows-type design weights \(\eta_i \in (0,
1]\) for robustness to high-leverage design points (the classical
gap of the RSE / GM-estimator: the \(\psi\)-functions bound the influence of
\(y\)-outliers but not of extreme
covariates). "mcd" computes \(\eta\) from robust Mahalanobis distances of
the non-constant columns of X; a numeric vector or
NULL (default, exact previous behaviour) are also accepted.
The weight enters the e-side of every estimating equation (\(\beta, u, \sigma,
\theta\)), the model vcov, the Satterthwaite degrees
of freedom and the influence diagnostics; summary reports
how many observations are downweighted. In a leverage- contamination
study (10% of points moved to a contaminating line at $\(7 robust SD) the plain RSE breaks down with `lmer`
while
Mallows-RSE stays near the clean estimate (\)\(4x less slope bias
at 5% contamination), at a clean-model efficiency cost of about 1% at
the simulation-backed default tuning (\)= 1$, cutoff \(0.975\)). Active design weights are
supported for a single grouping factor; rlmer stops with an
error on models with more than one.diag(f | g),
cs(f | g) and
ar1(f | g) random-effect covariance
structures are now supported (both the heterogeneous and homogeneous /
equal-variance forms of cs and ar1;
diag() was already accepted in 3.4-5, cs() and
ar1() are new). diag() (uncorrelated random
slopes, equivalent to (f || g)) is fitted directly.
cs() (compound symmetric: one common correlation) and
ar1() (autoregressive:
Cor(i, j) = rho^|i - j|) are fitted by projecting each
block’s unstructured scoring-equation update onto the structured
manifold every iteration – the marginal variances are kept and the
structure’s single correlation parameter is estimated from the block’s
robust correlation matrix: the average off-diagonal correlation for
cs (the constrained estimator for that linear structure,
Anderson 1973), and for ar1 the conditional maximiser of
the working-Gaussian objective over the manifold – an ECME /
conditional-maximisation step (Meng & Rubin 1993; Liu & Rubin
1994) that makes the projected fixed point a stationary point of the
constrained objective and, unlike the former log-correlation
least-squares fit, recovers a negative rho. So the
algorithm converges to a genuinely structured fit. For a 2x2 block (a
single correlation) cs() reproduces the unstructured fit to
numerical precision (1e-7 in the regression test).
VarCorr() reports the structured (sd, corr)
parametrisation. DAStau is limited to blocks of size <=
2 (it falls back to DASvar for larger blocks, as before);
cs/ar1 blocks of size >= 3 are therefore
fitted with DASvar. cs/ar1 remain
experimental (the projection has no consistency proof) and are validated
only for a single grouping factor: rlmer now
warns when a cs/ar1 term is
combined with more than one grouping factor, where neither the fit nor
its inference is validated. Inference degrades gracefully –
summary(df = "satterthwaite") reports the robust df where
the constrained covariance admits it and otherwise falls back to the
t-table. A projected cs/ar1 block of dimension
nc requires more than nc observations per
cluster (within-cluster replication); on a saturated block (cluster size
<= nc, e.g. one observation per visit for a 4-level
ar1) the structured fit is not identifiable, and
rlmer now stops up front with an
informative error instead of failing cryptically deep in the
projection.
Experimental Monte-Carlo DAS-tau calibration
(options(robustlmm.dastau.mc = TRUE)).
method = "DAStau" computes its consistency factors by
Gauss-Hermite quadrature, which is limited to random-effect blocks of
dimension <= 2; fits containing a larger block – including
cs/ar1/unstructured blocks of size >= 3 –
fall back to DASvar with a warning. The new option instead
computes the same self-consistent DAS-tau fixed point by plain
Monte-Carlo integration for those blocks. The MC sample is drawn once
per fit (common random numbers, deterministically seeded via
robustlmm.dasmc.seed, moment-matched), so fits are
reproducible and the caller’s .Random.seed is untouched;
the number of draws is robustlmm.dasmc.nsim (default
1e5). Simulation-validated on clean Gaussian data (in the
companion ar1 study it removes the small DASvar calibration
residual of ~0.004 in the fitted correlation at 4-level ar1
blocks, J = 200), but not backed by a finite-sample theorem. Default
behaviour is unchanged: with the option off (the default) the
DASvar fallback fires exactly as before, and blocks of
dimension <= 2 keep the Gauss-Hermite path even with the option on
(the MC path offers no improvement there). See
?"robustlmm-options".
vcov_sandwich emits a warning when
n.clusters < 20. In a simulation study at \(J = 8\), the sandwich CI undercovers (~0.89
vs nominal 0.95) and the Wald anova Type-I inflates to
~0.15-0.20 (3-4x nominal). G1 correction is necessary but not sufficient
at small \(J\); prefer
type = "default" or pair the sandwich CI with the
confintROB bootstrap. Caveat paragraph added to
?vcov.rlmerMod, ?vcov_sandwich,
?confint.rlmerMod, and
?anova.rlmerMod.
anova(., test = "boot") is well-calibrated under
clean Gaussian and t3 heavy-tailed errors but
anti-conservative under subject-level contamination (Type-I 0.13-0.24 vs
nominal 0.05 at 10% subjects shifted +5 sigma). It now warns
automatically when the null fit heavily downweights a
random-effects group (smallest robustness weight below 0.5) – a direct
signal that the fit treated that group as an outlier, so the bootstrap
null built from it may be poisoned. The robust fit absorbs a
contaminated group into a downweighted random effect, so
cooks.distance (an influence measure) does not reliably
flag it for this purpose; the random-effects weights are the effective
screen. The path remains experimental.
method = "DASvar" fits that use prior
weights were wrong, in every release from 3.2-0 to
3.4-5. If you have fitted a weighted model with
method = "DASvar" and any of those versions, refit it.
Unweighted fits are not affected, in any version. The default
method = "DAStau" is affected far more narrowly — see the
end of this entry.
Some history, because it bears on who is exposed.
rlmer() refused weights outright from 1.1
(2013) to 3.1-2 (2022) — “Argument weights is unsave to use at the
moment” — and 1.0 accepted them only with an explicit “Argument
weights is untested” warning. 3.2-0 removed that guard and enabled
the argument, but never documented it: weights is not in
rlmer()’s formals and has no entry in ?rlmer;
it reaches lmer() through .... So from 3.2-0
on, weighted fits were accepted in silence, with
tests/weights.R pinning them against lmer —
and were wrong. They are documented, and correct, from 3.5.0.
The Design Adaptive Scale tau_e counted the prior
weights twice. The residuals entering it are already whitened by
U_e = diag(sqrt(w)), so the leading term of
tau^2 is 1 by construction; the formula used w
there instead. Because prior weights are identified only up to a common
factor, this was not a fixed bias: multiplying weights by a
constant c no longer multiplied sigma by
sqrt(c), so the answer depended on the arbitrary scale the
weights happened to be supplied on. On Dyestuff with
weights = 0.01 throughout, sigma came out 12
times too large; with weights = 100, 11 times too small.
With weights varying over the more usual range [0.25, 4] on a
250-observation one-way design the error was 14% low for
method = "DASvar" and 5% low for
method = "DAStau", and the default vcov
covered the slope 79.8% of the time at a nominal 95%. Nothing was
reported to the user: the fitting warnings this produced were collected
into fit@optinfo$warnings and never printed.
What changes on a refit: sigma, VarCorr,
every standard error, and therefore summary,
confint, vcov and the anova
tables. Under a robust rho the fixed-effect estimates also
move slightly, because the robustness weights are evaluated at
r / (sigma * tau). Repaired coverage on the design above is
94.5%.
How much of this reaches method = "DAStau", the
default. Far less. DAStau does not use the DASvar formula as
its answer, only as the starting value of its fixed-point iteration, and
the starting value mattered only because calcTau refused to
iterate at all when it came in above 1 (fixed in 3.5.0-2, not here).
Whenever the fit re-enters calcTau even once more, it
warm-starts below 1 and converges to the correct value, so the damage
undoes itself. Over a grid of 3 designs x 6 weight regimes x 2
rho x 2 starting values, method = "DASvar" is
wrong in 70 of 70 weighted arms (by up to 109%) while
method = "DAStau" is wrong in 10 of 72 — all ten of them at
rho.e = rho.b = cPsi, which is the diagnostic setting that
makes rlmer() reproduce lmer(), not a setting
for robust estimation. With the default smoothPsi there was
no affected arm at all (36 of 36 clean, worst deviation 1.3e-07): the
outer loop re-enters calcTau 13 times on the check dataset,
and only the first entry sees an inflated start. So in practice: refit
weighted DASvar fits unconditionally, and weighted
DAStau fits if you ran them at cPsi. This is
measured over a finite grid, not proved — if a weighted
DAStau fit matters to you, refitting under 3.5.0 and
comparing costs less than reasoning about it.
The fix is pinned by tests/weights.R, which now runs
both method = "DAStau" and method = "DASvar"
against lmer at rho.e = rho.b = cPsi – the
setting in which rlmer() reduces to REML and must reproduce
lmer exactly. Previously it exercised the default method
only, which is why this went unnoticed. Versions before 3.2-0 are
unaffected because they refused weights outright; 3.2-0 is
where weights became a supported argument.
Random-effect robustness weights
(getME(., "w_b_vector"), "w_sigma_b_vector")
are now aligned with the spherical random effects ("b_s").
Previously they were concatenated per covariance block type, which
silently permuted them for RE terms expanding into several block types –
notably heterogeneous diag() terms (lme4 >= 2.0-0).
Consequences fixed: fitEffects() applied the permuted
weights inside the fitting iterations (estimates for diag()
fits may change slightly), getME(., "w_b") /
"w_sigma_b" attributed weights to the wrong levels,
plot() / qq() colored points by the wrong
weights, and the anova(test = "boot") group guard and
null = "robust" trimming could flag/trim the wrong
clusters.
anova(test = "boot", null = "robust") now rebuilds
the bootstrap generator’s U_b = Lambda(theta) at the full
cleaned theta vector (exact for any number of variance components;
previously only single-theta fits were rescaled). The trimmed refit is
checked for parameter compatibility with the full fit, and the anova
table heading now states which fallback applied (layout unresolved /
>50% flagged / refit unusable) instead of implying no cluster was
flagged.
tests/anova.R: single-fit Wald shape; pairwise Wald
matches single-fit Wald for the same restriction; VC-test warning steers
to bootstrap; reproducibility under seed.robustlmm-3.5.0 (“New Features in
robustlmm 3.5.0”): a feature tour with worked examples of the sandwich
vcov, Wald / bootstrap confint, robust
anova, prediction intervals, Satterthwaite degrees of
freedom, observation- and cluster-level influence diagnostics, the
RANSAC initial estimator, Mallows design weights and the structured
(diag/cs/ar1) covariances. Every
example runs on a small built-in dataset.lemon dropped from Suggests. Its
geom_pointline was soft-deprecated under ggplot2 4.0; the
simulation-study plots now use the equivalent
ggh4x::geom_pointpath (ggh4x was already a
dependency).ggplot2 deprecation cleanup for ggplot2 4.0:
plot.rlmerMod no longer uses the deprecated
aes_string() (now the .data[[]] tidy-eval
idiom), and the simulation-study scripts use
ggh4x::geom_pointpath in place of the soft-deprecated
lemon::geom_pointline. The vignettes build without ggplot2
deprecation warnings.
inst/simulationStudy/inferenceStudy.R and its
results writeup: the inference simulation study (CI coverage, Type-I,
power under contamination, and the variance-component bootstrap test).
Kept on the development branch as the research record and excluded from
the package build via .Rbuildignore; the user-facing
conclusions are distilled into the ?anova.rlmerMod /
vcov_sandwich caveats.
generateRepeatedMeasuresDatasets(): new data
generator for balanced repeated-measures datasets with a within-subject
factor and a structured
(cs/ar1/diag/unstructured)
random-effect covariance, matching the interface of
generateMixedEffectDatasets().
fitDatasets_rlmer_cs() /
fitDatasets_rlmer_ar1(): simulation-study fitting wrappers
that force a cs / ar1 random-effects
covariance structure (rewriting the single random-effects term
regardless of how the data were generated), plus
inst/simulationStudy/structuredCovariances.R, a
consistency-and-robustness study showing that under random-effects
contamination the classical (lmer) structured correlation
collapses toward zero while the robust structured fit stays near the
truth. Study script and results are excluded from the build via
.Rbuildignore.
anova(f0, f1) with random-effects difference)
requires test = "boot"; the parametric bootstrap is
currently the only exposed VC path. The residual-bootstrap variant
explored in a companion simulation study is sim-only at this stage.anova(f0, f1, f2, ...) are not yet
supported; only pairwise comparison.The features below shipped to CRAN in 3.4-4 (2026-06-21) and its build fix 3.4-5 (2026-06-21). They are listed here so the changelog matches what CRAN users of those versions already have; 3.5.0 extends several of them (see above).
vcov(., type = "sandwich"). Robust
cluster-sandwich covariance of the fixed effects, \(\hat V = \hat A^{-1} \hat B \hat A^{-T}\)
with \(\hat A\) the Schur-complement
(marginal) \(\beta\)-Jacobian and \(\hat B = \sum_j s_j s_j^T\) the per-cluster
\(\beta\)-score contributions. Defaults
are preserved byte-for-byte (type = "default" returns the
pre-existing lme4-inherited linearised vcov).
correction = "G1" (default) applies the \(J/(J-1)\) small-sample scaling.
cluster auto-detects the sole grouping factor; on crossed
designs cluster = "factor-name" warns that the sandwich is
approximate.
confint(., method = "Wald").
Closed-form per-coefficient Wald CI from
vcov(., type = vcov_type). The default
vcov_type = "default" matches
confint(lmer_fit, method = "Wald") on clean data and is
sharper under error-distribution contamination (20-40% shorter intervals
at matched coverage under t3 and CN_e). Bare
confint(fit) previously errored on dpoMatrix
vcov; now returns the Wald interval.
confint(., method = "boot" / "BCa").
Delegates to confintROB::confintROB (Mason, Cantoni &
Ghisletta 2021, 2024) with boot.type = "parametric" (the
3.4-5 default; 3.5.0 changes the default to "wild") or
"wild". confintROB is in
Suggests; requireNamespace guard emits a clear
error if it is not installed.
predict(., interval = "confidence" / "prediction", level = 0.95).
Returns (fit, lwr, upr, se) data.frame. Fixed- effect SE
from vcov(., type = "sandwich"); RE contribution from the
partial influence function of \(\hat
u\). Default interval = "none" preserves the
pre-existing numeric vector return.
cooks.distance.rlmerMod.
Per-observation joint Mahalanobis influence on \((\hat\beta, \hat\sigma, \hat\theta)\),
computed from implicitIF_full. Backed by
numDeriv::jacobian (added to Imports). The most influential
Penicillin observations are 85, 86, 122, 123 and 13.
influence.rlmerMod. Per-observation
IF matrix wrapper around implicitIF_full.
caseweightIF,
vcov_sandwich,
implicitIF_full exported as user-facing
helpers (the IF substrate is also useful outside the vcov /
confint / anova paths).
emmeans support:
recover_data / emm_basis methods are
registered for rlmerMod, so emmeans /
emtrends / contrast work on robust fits
(asymptotic, Inf degrees of freedom; 3.5.0 adds finite
Satterthwaite df).
init = "ransac" for
rlmer(). String form for the RANSAC high-breakdown start:
equivalent to init = ransac_lme4(formula, data)$fit with
default K = 200, sub_frac = 0.5. Errors with a
clear message if no subsample fit succeeds.
size_obr = FALSE for
rlmer(), an opt-in experiment, off by default. When
TRUE, the size-controlling weight w_delta in
the block-diagonal variance-components scoring equation is replaced by
the Hampel-OBR form \(w_\tau(d) = \min(1,
b_\tau / |d - s\kappa - a|)\) (Stahel 1987; Hampel et al. 1986),
with \(b_\tau\) from the
rho.sigma.b tuning constant and \(a\) fixed by Fisher consistency under \(\chi^2_s\). It has no effect for diagonal
\(V_b\) (block size 1). Documented in
?rlmer, tested by tests/stahel-obr.R, and
available to simulation studies as
fitDatasets_rlmer_DAStau_sizeOBR.
Setting it TRUE is not a
recommendation. The default finite-difference
w_delta is not Hampel-OBR-optimal at the central model, and
this option does recover asymptotic OBR efficiency for the
variance-component magnitudes, but the gain at matched gross-error
sensitivity is a typically modest 1-2 percentage points. The default is
unchanged and Proposal-II remains the default everywhere. The option
exists so the comparison can be made rather than argued about; the shape
weight, by contrast, is already OBR-optimal and was left alone.
tests/influence.R: default-vcov-unchanged guard;
caseweightIF vs reweight-FD agrees to
< 1e-3 on Dyestuff/sleepstudy; sandwich finite + PD;
crossed-design warning; predict intervals; cooks.distance
top-5 match.tests/bootstrapWald.R: Wald closed-form match;
confintROB delegation including the parm-NULL fix
(confintROB returns a 0-row matrix on explicit parm = NULL;
the wrapper omits parm in the delegated call and subsets
rows itself).tests/cache-invalidation.R: .scoreVec’s
.tau_e / .Tbk cache-clearing guard for the
full IF.tests/vcov-jsweep-smoke.R: minimal sandwich-coverage
smoke at \(J = 18\).tests/sandwich-vs-bootstrap.R: sandwich SE vs
Monte-Carlo SE cross-check on sleepstudy.numDeriv promoted from absent to Imports (used by
implicitIF_full).