anova.rlmerMod: Analysis of variance for an rlmer fit.

In robustlmm: Robust Linear Mixed Effects Models

Analysis of variance for an rlmer fit.

Description

Three modes.

Single fit anova(fit)

returns a per-term robust Wald chi-square table from vcov(fit, type = vcov_type). Each term is tested marginally as T_t = \hat\beta_t^\top \hat V_t^{-1} \hat\beta_t \sim \chi^2_{k_t} under H_0: \beta_t = 0.

Nested fits differing only in fixed effects (default test = "Wald")

robust Wald restriction test on the extra coefficients. The same V used by the single-fit table is used here; vcov_type = "sandwich" carries through.

Nested fits differing in random-effects structure

the Wald and score asymptotics are invalid (boundary problem; see Self-Liang 1987 and Koller 2026 paper 2). test = "Wald" warns and switches to the parametric bootstrap automatically; test = "boot" is the default valid path. The bootstrap simulates nsim datasets from the fitted central LMM at fit0's estimates, refits both rlmer models per replicate, and uses the quasi-deviance difference D = 2 \sum [\rho_e(r_i^{(0)}/\hat\sigma_0) - \rho_e(r_i^{(1)}/\hat\sigma_1)] as the discrepancy statistic (Heritier and Ronchetti 1994; Cantoni and Ronchetti 2001; Heritier, Cantoni, Copt and Victoria-Feser 2009 sec. 5). For the special case of a single added variance component the experimental test = "score" offers a contamination-robust alternative (see Details).

Usage

## S3 method for class 'rlmerMod'
anova(
  object,
  ...,
  test = c("Wald", "boot", "score"),
  null = c("parametric", "robust"),
  vcov_type = c("default", "sandwich"),
  ddf = c("none", "satterthwaite"),
  nsim = 1000L,
  seed = NULL,
  verbose = FALSE
)

Arguments

object	An rlmerMod object.
...	A second rlmerMod object for pairwise comparison.
test	One of "Wald" (default; closed form for nested fixed-effects-only tests), "boot" (parametric bootstrap quasi-deviance; the default valid path for variance-component tests), or "score" (experimental one-sided robust score test for a single added variance component, calibrated by a score-only parametric bootstrap; see Details for its scope and validation record).
null	Bootstrap null generation for test = "boot". "parametric" (default) generates from fit0's robust estimates (the exact current behaviour). "robust" is an experimental contamination-robust null that generates from a cleaned null fit (see Details); it only affects the variance-component bootstrap path and is ignored for the Wald paths.
vcov_type	Forwarded to vcov for the Wald paths; ignored for the bootstrap.
ddf	Denominator degrees of freedom for the Wald paths. "none" (default) reports the chi-square table as before; "satterthwaite" reports an F-test with a Satterthwaite denominator df (the multivariate lmerTest::contestMD generalisation, built on the robust IF-based covariance of the variance parameters). It requires vcov_type = "default" and a single grouping factor; otherwise it warns and falls back to the chi-square table.
nsim	Bootstrap replicates when test = "boot" or test = "score"; default 1000 for "boot" and 199 for "score" (whose null-only refits are ~4x cheaper).
seed	Optional RNG seed for reproducibility of the bootstrap.
verbose	Bootstrap progress messages.

Details

Only pairwise comparison is implemented; chains anova(fit0, fit1, fit2, ...) are not yet supported.

Small-J caveat for the sandwich. vcov_type = "sandwich" is markedly anti-conservative at small J: in a simulation study the pairwise Wald Type-I rate reached 0.15-0.20 (vs. nominal 0.05) at J = 8 and was still elevated at J = 18, returning to nominal only by J \gtrsim 50. vcov_sandwich emits a warning for J < 20; prefer vcov_type = "default" for hypothesis tests at small J.

Subject-contamination caveat for the bootstrap (experimental). The parametric-bootstrap quasi-deviance path (test = "boot") is experimental — it is currently the only exposed variance-component path. It calibrates correctly under clean Gaussian and heavy-tailed errors but is markedly anti-conservative when a subset of subjects is contaminated. In a simulation study with 10% of subjects shifted by 5 standard errors (shift_subj), Type-I climbed to 0.13-0.24 (2.5-5x nominal) across J \in \{8, 18, 50\}. The mechanism: the bootstrap simulates from the fitted central LMM at fit0's estimates, so any bias fit0's \hat\sigma_0, \hat\theta_0 picked up from the contaminated subjects propagates into a too-narrow bootstrap null. As an automatic guard, anova(test = "boot") now inspects the null fit's own random-effects robustness weights and warns when a group is heavily downweighted (smallest weight below 0.5) — a direct signal that the fit treated that group as an outlier, so the bootstrap null may be poisoned. (The robust fit absorbs a contaminated group into a downweighted random effect, which is why cooks.distance — an influence measure — does not reliably flag it here: the downweighting that removes the bias also removes the influence.) When the warning fires, treat the p-value with caution.

Experimental robust null (null = "robust"). As a mitigation for the subject-contamination anti-conservativeness above, the bootstrap path accepts null = "robust": instead of generating the parametric bootstrap from fit0's raw estimates, it generates from a contamination-cleaned null fit. Clusters the robust fit heavily downweights (smallest per-cluster random-effect weight below a tuned threshold, the same signal the automatic guard uses) are flagged, fit0 is refitted with those clusters removed, and the bootstrap is generated from the resulting de-biased (\hat\sigma, \hat\theta, \hat\beta) (with the generator's U_b = \Lambda(\theta) rebuilt at the full cleaned \hat\theta vector). The observed discrepancy D is still taken from the original untrimmed fits, so the test uses all the data; only the bootstrap null generation is cleaned. If no cluster is flagged, more than half the clusters are flagged, trimming would leave the design rank-deficient (a dropped contrast level, a collapsed grouping factor, or a constant random-slope covariate — checked by the same nonsingular-subsampling test used by the RANSAC initial estimator), or the trimmed refit fails, the call falls back to the plain parametric null (so on clean data null = "robust" closely matches null = "parametric"); the table heading states which case applied. In simulation this reduces the contamination-induced Type-I inflation at larger J while preserving power, at the cost of mild conservatism on clean data. It is experimental, validated by simulation rather than a finite-sample theorem; the table heading notes when the robust null was applied.

Experimental robust score test (test = "score"). For the common special case where the alternative adds exactly one independent scalar variance component relative to the null — e.g. (1|g) vs (1|g) + (0 + x|g), or diagonal structures adding one component — test = "score" runs a one-sided robust score test computed from the robust null fit only and calibrated by a score-only parametric bootstrap that refits just the null model per replicate (roughly 4x cheaper than test = "boot"'s double refits; hence the smaller default nsim = 199). Per cluster j, the whitened marginal residuals \tilde r_j = V_j^{-1/2}(y_j - X_j \hat\beta_0) and the whitened tested direction v_j = V_j^{-1/2} z_j (with V_j the null fit's marginal covariance and z_j the added term's design column) give the bounded contribution s_j = (v_j^\top \psi(\tilde r_j))^2 - \kappa_1 \|v_j\|^2 (\psi the fit's rho.e psi-function, \kappa_1 = E[\psi(Z)^2]); the statistic is the self-normalised cluster sum S = \sum_j s_j / \sqrt{\sum_j (s_j - \bar s)^2}, and the one-sided p-value is (1 + \#\{S^* \ge S\})/(n_{\mathrm{eff}} + 1). Because each s_j is psi-bounded, whole-cluster contamination shifts the statistic and its bootstrap reference law together instead of inflating the test: in simulation (Gaussian balanced designs, one scalar tested component; not a theorem) contaminated-null Type-I was 0.035 vs 0.115 for test = "boot" with 10% of clusters shifted at J = 50, clean-null Type-I 0.045, and power 0.920 vs 0.900 — with no null cleaning needed (raw vs cleaned generating parameters changed at most 1-2 decisions in 600+ paired replicates, so null = "robust" is ignored and no downweighted-group warning is issued for this path). Across an adversarial sweep the pattern held: J = 18 clean/contaminated 0.055/0.015, uncentered-x designs clean/contaminated 0.045/0.075, 5% single-observation outliers 0.060. Scope: both fits on the same data with identical fixed effects and a single shared grouping factor, and the alternative adds exactly one uncorrelated scalar component (length(theta) differs by 1); anything else — multi-component or correlated-slope alternatives, crossed or nested factors — stops with an error, use test = "boot" there. Caveat: contamination aligned with the tested direction is indistinguishable from the alternative for any test with power; inspect the attached per-cluster contributions (attr(., "boot")$s_j; the largest values identify the clusters driving the statistic) and cluster diagnostics (cooks.distance(fit, groups = ), hatvalues) when a rejection is suspect. Experimental: simulation-validated, not proven; the deviance bootstrap test = "boot" remains the default variance-component path.

Value

An "anova" data.frame; the bootstrap path attaches attr(., "boot") = list(D_boot, n_fail, nsim, D_obs). The score path attaches attr(., "boot") = list(S_boot, n_fail, nsim, S_obs, s_j) with s_j the named per-cluster contributions to the observed statistic.

References

Heritier S, Ronchetti E (1994). Robust bounded-influence tests in general parametric models. JASA 89(427): 897–904.

Cantoni E, Ronchetti E (2001). Robust inference for generalized linear models. JASA 96(455): 1022–1030.

Heritier S, Cantoni E, Copt S, Victoria-Feser MP (2009). Robust Methods in Biostatistics. Wiley.

See Also

vcov, confint

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