anova.rlmerMod: 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 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
)

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.

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.

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.