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
$\sim$7 robust SD) the plain RSE breaks down with lmer while
Mallows-RSE stays near the clean estimate ($\sim$4x less slope bias
at 5% contamination), at a clean-model efficiency cost of about 1% at
the simulation-backed default tuning ($\gamma = 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).