NEWS.md

In robustlmm: Robust Linear Mixed Effects Models

robustlmm NEWS

robustlmm 3.5.0-2 (2026-07-30)

Bug fixes

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

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.

robustlmm 3.5.0-1 (2026-07-10)

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.

New user-facing inference methods

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

Robustness to high-leverage design points

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

Structured random-effect covariances (lme4 >= 2.0-0)

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

Small-J caveats and warnings

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

Bug fixes and improvements

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

Regression suite

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

Documentation

• New vignette 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.

Dependencies

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

Internal / infrastructure

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

Known limitations

• VC test (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.
• Multi-fit chains anova(f0, f1, f2, ...) are not yet supported; only pairwise comparison.

robustlmm 3.4-4 / 3.4-5 (2026-06-21, released to CRAN)

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

New user-facing inference methods

• 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). Top-flagged Penicillin observations match the IF-thread1 prototype's documented set {85, 86, 122, 123, 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.

Regression suite

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

Dependencies

• numDeriv promoted from absent to Imports (used by implicitIF_full).

robustlmm documentation built on July 30, 2026, 5:11 p.m.