# 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 `NA`s 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`).
