% !Rnw weave = Sweave
\documentclass[nojss]{jss}
\usepackage{rlmer}
\usepackage{booktabs}

%\VignetteIndexEntry{New features in robustlmm 3.5.0: robust inference, diagnostics and structured covariances}
%\VignetteDepends{robustlmm, lme4, emmeans}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% declarations for jss.cls %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\author{Manuel Koller}
\title{New Features in \pkg{robustlmm} 3.5.0: Robust Inference,
  Diagnostics and Structured Covariances}
\Plaintitle{New Features in robustlmm 3.5.0: Robust Inference,
  Diagnostics and Structured Covariances}
\Shorttitle{New Features in \pkg{robustlmm} 3.5.0}

\Abstract{The 3.5.0 release of \pkg{robustlmm} adds a layer of
  inference, diagnostic and modelling tools on top of the Robust
  Scoring Equations estimator introduced in the main package vignette.
  This companion vignette is a feature tour. It shows how to obtain a
  robust cluster-sandwich covariance for the fixed effects, Wald and
  bootstrap confidence intervals, robust Wald and bootstrap ANOVA
  tables, prediction intervals, and honest Satterthwaite degrees of
  freedom that stay valid under contamination. It then covers the
  influence diagnostics
  (observation- and cluster-level Cook's distances and robust
  leverages), a high-breakdown \textsc{ransac} initial estimator,
  Mallows-type design weights that bound the influence of
  high-leverage design points, and support for structured
  (\code{diag}, \code{cs}, \code{ar1}) random-effect covariances. Every
  example runs on a small built-in dataset so that the methods can be
  reproduced interactively. For the model and the estimation algorithm
  itself, see the main vignette, \code{vignette("rlmer")}.}

\Keywords{robust statistics, mixed-effects model, inference,
  influence, Cook's distance, sandwich, bootstrap, Satterthwaite,
  \textsc{ransac}, design weights, compound symmetry, \proglang{R}}
\Plainkeywords{robust statistics, mixed-effects model, inference,
  influence, Cook's distance, sandwich, bootstrap, Satterthwaite,
  RANSAC, design weights, compound symmetry, R}

\Address{
  Manuel Koller\\
  E-mail: \email{kollerma@proton.me}\\
  URL: \url{https://github.com/kollerma/robustlmm}
}

%% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\SweaveOpts{prefix.string=text-350}

\begin{document}
\SweaveOpts{concordance=TRUE}

\baselineskip 15pt

<<preliminaries, echo=FALSE, results=hide>>=
options(prompt = "R> ", continue = "+  ", width = 76,
        useFancyQuotes = FALSE, digits = 4)
suppressPackageStartupMessages(library("robustlmm"))
set.seed(20260623)
## Numbers quoted from the simulation studies are precomputed into this
## RDS by vignettes/precompute-3.5.0.R and pulled in with \Sexpr{} below,
## so they always match the studies and are never hand-typed.
nums <- readRDS("robustlmm-3.5.0-numbers.rds")
## numbers of the MMRM section, written by vignettes/precompute-mmrm.R
mm <- readRDS("mmrm-sensitivity-numbers.rds")
fmt <- function(x, digits = 2) formatC(x, format = "f", digits = digits)
## The slowest chunks are not evaluated (eval=FALSE, echo=FALSE); their
## output is precomputed by vignettes/precompute-chunks.R (to keep the
## check time down) and showChunk() writes code and output as Sweave
## would. It stops if a chunk's code no longer matches the code the
## output was computed from.
chunkOutput <- if (file.exists("robustlmm-3.5.0-chunks.rds"))
    readRDS("robustlmm-3.5.0-chunks.rds")
showChunk <- function(label) {
    src <- readLines("robustlmm-3.5.0.Rnw")
    s <- grep(paste0("^<<", label, "[,>]"), src)
    e <- s + match(TRUE, grepl("^@", src[-seq_len(s)]))
    x <- chunkOutput[[label]]
    if (!identical(src[(s + 1):(e - 1)], x$code))
        stop("the stored output of chunk '", label, "' is out of date; ",
             "re-run vignettes/precompute-chunks.R")
    out <- "\\begin{Schunk}"
    input <- FALSE
    for (st in x$steps) {
        if (!input) out <- c(out, "\\begin{Sinput}")
        out <- c(out, paste0(c(getOption("prompt"),
                               rep(getOption("continue"),
                                   length(st$src) - 1L)), st$src))
        input <- TRUE
        if (length(st$output)) {
            out <- c(out, "\\end{Sinput}", "\\begin{Soutput}", st$output,
                     "\\end{Soutput}")
            input <- FALSE
        }
    }
    if (input) out <- c(out, "\\end{Sinput}")
    writeLines(c(out, "\\end{Schunk}"))
}
## Sweave does not echo R warnings into the vignette the way an
## interactive session does; this helper prints them next to the call
## that triggers them.
echoWarnings <- function(expr)
    withCallingHandlers(expr, warning = function(w) {
        writeLines(strwrap(paste("Warning:", conditionMessage(w)),
                           width = 74))
        invokeRestart("muffleWarning")
    })
@

\section{Introduction}\label{sec:intro}

The main \pkg{robustlmm} vignette \citep{vignette} introduces the
Robust Scoring Equations (RSE) estimator and shows how to fit a model
with \code{rlmer}, tune its robustness and efficiency, and read the
robustness weights. This companion vignette documents the features
added for the 3.5.0 release. Several of the foundations---the sandwich
covariance, the Wald and bootstrap confidence intervals, the prediction
intervals, the observation-level influence diagnostics and the base
\textsc{ransac} start---first shipped to CRAN in releases 3.4-4/3.4-5;
this vignette documents the full series. The features fall into four
groups:

\begin{enumerate}
\item \textbf{Inference for the fixed effects}
  (Section~\ref{sec:inference}): a robust cluster-sandwich covariance
  (\code{vcov(., type = "sandwich")}), Wald and bootstrap confidence
  intervals (\code{confint}), robust Wald and bootstrap ANOVA tables
  (\code{anova}), prediction intervals (\code{predict}), and honest
  Satterthwaite degrees of freedom and $p$~values
  (\code{summary(., df = ...)}).
\item \textbf{Influence diagnostics} (Section~\ref{sec:influence}):
  observation- and cluster-level Cook's distances
  (\code{cooks.distance}), the full per-observation influence matrix
  (\code{influence}) and robust leverages (\code{hatvalues}).
\item \textbf{A high-breakdown initial estimator}
  (Section~\ref{sec:ransac}): a \textsc{ransac} start
  (\code{init = "ransac"}, \code{ransac\_lme4}, \code{rlmer\_ransac}).
\item \textbf{Robustness in the design space}
  (Section~\ref{sec:design}): Mallows-type design weights
  (\code{design.weights}) and structured random-effect covariances
  (\code{diag}, \code{cs}, \code{ar1}, Section~\ref{sec:structured}).
\end{enumerate}

The \code{cs}/\code{ar1} covariances (group~4) are \emph{experimental}:
they are validated by simulation but rest on newer, partly heuristic
machinery, and their interface or defaults may still change; their
section title is marked accordingly. The \textsc{ransac} start (group~3)
and the Mallows design weights (group~4) are stable: settled interfaces
with simulation-backed defaults.

All examples use the small built-in datasets \code{sleepstudy},
\code{Dyestuff} and \code{Penicillin} from \pkg{lme4} so that they run
in seconds. Throughout we fit with \code{method = "DASvar"}, the fast
direct approximation, to keep the vignette quick to build; for final
analyses the default \code{method = "DAStau"} is more accurate (see the
main vignette). We start from a robust fit of the \code{sleepstudy}
data, a random-slope model for reaction time as a function of days of
sleep deprivation:

<<base-fit>>=
fm <- rlmer(Reaction ~ Days + (Days | Subject), sleepstudy,
            method = "DASvar")
@

\section{Inference for the fixed effects}\label{sec:inference}

\subsection{Robust cluster-sandwich covariance}\label{sec:sandwich}

The fixed-effect covariance that \code{summary} reports by default is
the linearised (model-based) covariance inherited from the \pkg{lme4}
machinery. \code{vcov} now takes a \code{type} argument that selects
between this default and a robust \emph{cluster-sandwich} estimator,
$\hat V = \hat A\inv \hat B \hat A^{-\top}$, where $\hat A$ is the
marginal (Schur-complement) $\beta$-Jacobian and
$\hat B = \sum_j \vec s_j \vec s_j\tr$ sums the per-cluster
$\beta$-score contributions. The sandwich is heteroscedasticity- and
misspecification-robust and does not assume the fitted covariance is
correct.

<<vcov-sandwich>>=
vcov(fm, type = "default")
echoWarnings(vcov(fm, type = "sandwich"))
@

The default is preserved byte-for-byte, so existing code is
unaffected. The small-sample correction \code{correction = "G1"}
(default) applies the $J / (J - 1)$ scaling. On a design with a single
grouping factor the cluster variable is detected automatically; with
several grouping factors pass \code{cluster = "factor-name"} (also to
\code{confint} and \code{anova}): with the top level of a nested design
the sandwich is exact, on crossed designs only approximate.

\emph{A small-$J$ caveat.} The sandwich is a large-$J$ approximation.
\code{sleepstudy} has only $J = \Sexpr{nlevels(sleepstudy$Subject)}$
subjects --- fewer than the $20$ clusters below which \code{vcov_sandwich}
warns --- which is why the call above prints the warning. With few
clusters the G1 sandwich is biased down for between-cluster coefficients
and has fewer effective degrees of freedom than $J - 1$, so normal or
$\chi^2$ references on it are anti-conservative. \code{correction =
"CR2"} gives the bias-reduced sandwich \citep{bell2002bias} adapted to the
robust estimating equations, and \code{confint} and \code{anova} use it
by default for \code{vcov\_type = "sandwich"}, with Bell--McCaffrey /
approximate-Hotelling degrees of freedom \citep{pustejovsky2018small}
(Sections~\ref{sec:confint} and~\ref{sec:anova}); the simulation
evidence is in the section ``Small-sample references for the cluster
sandwich'' of \code{vignette("simulationStudies")}.

\subsection{Confidence intervals}\label{sec:confint}

\code{confint} gains a \code{method} argument. The default
\code{method = "Wald"} returns the closed-form per-coefficient interval
$\hat\beta_k \pm q\sqrt{\hat V_{kk}}$ from \code{vcov(., type =
  vcov\_type)}. Bare \code{confint(fit)} used to error on the
\code{dpoMatrix} covariance; it now returns the Wald interval.

<<confint-wald>>=
confint(fm, method = "Wald")
confint(fm, method = "Wald", vcov_type = "sandwich")
@

For the default covariance the quantile $q$ is the normal
$z_{1 - \alpha/2}$; for the sandwich it is a Bell--McCaffrey $t$
quantile per coefficient with the CR2 covariance (its degrees of freedom
are in \code{attr(., "df")}). Passing
\code{df = "satterthwaite"} substitutes the robust Satterthwaite
$t$-quantile (Section~\ref{sec:df}), which widens the interval at small
$J$:

<<confint-df, eval=FALSE, echo=FALSE>>=
confint(fm, method = "Wald", df = "satterthwaite")
@
<<confint-df-output, echo=FALSE, results=tex, strip.white=false>>=
showChunk("confint-df")
@

When the asymptotic $\chi^2_p$ limit is suspect---typically at small
$J$---a bootstrap is preferable. \code{confint} delegates
\code{method = "boot"} and \code{method = "BCa"} to the \pkg{confintROB}
package \citep{mason2021confintrob, mason2024confintrob}, which
implements the parametric and the wild bootstrap (\code{boot.type},
default \code{"wild"}), with optional bias-correction-and-acceleration
(\code{method = "BCa"}). \pkg{confintROB} is in \code{Suggests}, so
these paths require it to be installed; \code{method = "Wald"} does
not. The bootstrap is costly, so it is not run here:

<<confint-boot, eval=FALSE>>=
## parametric BCa bootstrap (needs the 'confintROB' package)
confint(fm, method = "BCa",  nsim = 1000, seed = 20260601)
## wild bootstrap
confint(fm, method = "boot", boot.type = "wild", nsim = 1000,
        seed = 20260601)
@

Practical guidance from those papers: the $\chi^2_p$ Wald limit is
adequate for $J \gtrsim 20$ groups, the wild bootstrap is robust to
misspecification of the response covariance, and BCa is preferable when
the bootstrap distribution is skewed.

\subsection{ANOVA: robust Wald and bootstrap tests}\label{sec:anova}

\code{anova} now works on \code{rlmer} fits. Called on a single fit it
prints a per-term robust Wald table:

<<anova-single>>=
anova(fm)
@

Pass \code{vcov\_type = "sandwich"} to base the table on the
cluster-sandwich covariance instead; the table is then the small-$J$
$F$-test of the CR2 sandwich with approximate-Hotelling denominator
degrees of freedom (\code{ddf = "none"} restores the $\chi^2$ table).
Called on two nested
fixed-effects models it performs a Wald restriction test:

<<anova-two>>=
fm0 <- rlmer(Reaction ~ 1 + (Days | Subject), sleepstudy,
             method = "DASvar")
anova(fm0, fm)
@

With \code{ddf = "satterthwaite"} the table reports an $F$-test with
robust Satterthwaite denominator degrees of freedom:

<<anova-ddf>>=
anova(fm, ddf = "satterthwaite")
@

For variance-component comparisons on the PSD-cone boundary,
\code{test = "boot"} runs a parametric-bootstrap quasi-deviance test
(using a common scale to avoid a sign-flip artefact). This refits the
model many times and so is not run here:

<<anova-boot, eval=FALSE>>=
anova(fm0, fm, test = "boot", nsim = 1000, seed = 20260601)
@

The bootstrap variance-component test is \emph{experimental}: it is
currently the only exposed variance-component path, and a simulation
study shows it becomes anti-conservative when whole groups are
contaminated. The test guards against this on its own: it warns
automatically when the null fit heavily downweights a random-effects
group (smallest random-effect robustness weight below 0.5), a direct
signal that the bootstrap null built from that fit may be poisoned;
\code{null = "robust"} (below) is the mitigation. The per-term Wald
table and the fixed-effect restriction test above are not affected.

Under group contamination the bootstrap null generation itself can be
poisoned by the contaminated estimates. The experimental
\code{null = "robust"} option,
\code{anova(fm0, fm, test = "boot", null = "robust")}, generates the
bootstrap samples from a contamination-cleaned null fit: clusters whose
smallest random-effect robustness weight falls below a tuned threshold
are trimmed for the generation step only, while the observed statistic
still uses all the data. In simulation this reduces --- but does not
fully remove --- the contamination-induced Type-I inflation at larger
$J$ while preserving power. It is validated by simulation.

For the common special case where the alternative adds exactly one
independent scalar variance component --- as in the example above,
\code{(1|Subject)} vs \code{(1|Subject) + (0 + Days|Subject)}, or
diagonal structures adding one component --- the experimental
\code{test = "score"}, \code{anova(fm0, fm, test = "score")}, is the
contamination-robust alternative. It computes a one-sided score
statistic from psi-bounded per-cluster contributions of the robust
\emph{null} fit only and calibrates it by a score-only parametric
bootstrap that refits just the null model per replicate, about four
times cheaper than \code{test = "boot"} (default \code{nsim = 199}).
In simulation the bounded statistic held its size where the deviance
bootstrap did not: with $\Sexpr{nums$scoreTest$contamPct}\%$ of the
clusters shifted at $J = \Sexpr{nums$scoreTest$J}$, its Type-I error was
$\Sexpr{fmt(nums$scoreTest$contamScore, 3)}$ against the deviance
bootstrap's $\Sexpr{fmt(nums$scoreTest$contamBoot, 3)}$, with clean-null
Type-I error $\Sexpr{fmt(nums$scoreTest$cleanScore, 3)}$
($\Sexpr{fmt(nums$scoreTest$cleanBoot, 3)}$) and power
$\Sexpr{fmt(nums$scoreTest$powerScore, 3)}$
($\Sexpr{fmt(nums$scoreTest$powerBoot, 3)}$), over
$\Sexpr{nums$scoreTest$nreps}$ samples per setting
(\code{inst/simulationStudy/scoreTest.R}, which reproduces these
settings end-to-end through the package interface). With 18 groups,
an uncentred covariate or $5\%$ outlying observations its Type-I error
stayed between $\Sexpr{fmt(nums$scoreTest$sweepSize[1], 3)}$ and
$\Sexpr{fmt(nums$scoreTest$sweepSize[2], 3)}$. Contamination aligned
with the tested direction is indistinguishable from the alternative for
any test with power (with shifted slopes in the contaminated clusters the
test rejected in
$\Sexpr{fmt(100 * nums$scoreTest$sweepSlopeDirection, 0)}\%$ of the
samples); the returned per-cluster contributions (\code{attr(., "boot")\$s\_j})
identify the driving clusters, to be inspected together with the
cluster diagnostics of Section~\ref{sec:influence} when a rejection is
suspect. Several added independent components are tested jointly (see
the section ``Several added variance components'' of
\code{vignette("simulationStudies")}). Outside its scope (correlated-slope
alternatives, crossed or nested grouping factors) it stops with an
error; \code{test = "boot"} remains the default variance-component
path. It is validated by simulation.

Under error-distribution contamination the robust Wald test can be
substantially more powerful than the classical likelihood-ratio test
while matching its Type-I error on clean data; a simulation study
quantifies the power gain (see the inference study in
\code{vignette("simulationStudies")}).

\subsection{Satterthwaite degrees of freedom}\label{sec:df}

By default \code{summary} of an \code{rlmer} fit reports a $t$-value
without a reference distribution, because the small-sample distribution
of the robust estimator is not pivotal. The 3.5.0 release adds robust
Satterthwaite degrees of freedom \citep{satterthwaite1946}. Unlike
\pkg{lmerTest} \citep{kuznetsova2017lmertest}, the degrees of freedom
derive from the robust influence-function-based covariance of the
variance parameters, so they stay honest under contamination.
\code{summary} gains a \code{df} argument; \code{df = "satterthwaite"}
adds \code{df} and \code{Pr(>|t|)} columns:

<<summary-df>>=
summary(fm, df = "satterthwaite")$coefficients
@

The degrees of freedom are consistent across all four inference
surfaces: \code{summary}, \code{confint(df = "satterthwaite")}
(Section~\ref{sec:confint}), \code{anova(ddf = "satterthwaite")}
(Section~\ref{sec:anova}), and \pkg{emmeans}, which now returns finite
degrees of freedom instead of \code{Inf}:

<<emmeans-df>>=
emmeans::emtrends(fm, ~ 1, var = "Days")
@

\code{summary} shows the degrees of freedom \emph{by default}
(\code{df = "auto"}) whenever it is cheap---either the underlying
influence function is already cached on the fit, or a deterministic,
dimension-only size workload is below the cutoff
\code{getOption("robustlmm.summary.df.max", 5000)} (machine-independent
by construction).
Otherwise it falls back to the historic $t$-value table with a one-line
note; \code{df = "satterthwaite"} forces the computation and
\code{df = "none"} disables it. The (expensive) influence function is
cached on the fit and shared by \code{summary}, \code{anova},
\code{confint}, \code{emmeans}, \code{cooks.distance} and the sandwich
\code{vcov}, so it is computed once. Single-factor, nested
(\code{(1 | school/class)}) and crossed (\code{(1 | s) + (1 | item)})
designs are supported; nested designs use a one-way sandwich over the
coarsest grouping factor and crossed designs a
Cameron--Gelbach--Miller multiway cluster-robust covariance
\citep{cameron2011robust}. When a variance component is estimated at
the boundary (zero), the degrees of freedom are computed conditional on
that component being held at zero. The package options are documented
on the \code{help("robustlmm-options")} help page.

The \emph{robust} Satterthwaite degrees of freedom are a new
construction. Their nominal coverage has been validated by the
package's own Monte-Carlo simulations --- for single-factor, nested,
crossed and variance-component-boundary designs --- rather than by an
external published study (see the inference study in
\code{vignette("simulationStudies")}). On a design they do not yet cover
(e.g.\ crossed random slopes) \code{summary} falls back to the plain
$t$-value table rather than reporting an unvalidated number.

\subsection{Prediction intervals}\label{sec:predict}

\code{predict} gains an \code{interval} argument. The default
\code{interval = "none"} preserves the previous numeric-vector return;
\code{"confidence"} and \code{"prediction"} return a data frame with
columns \code{fit}, \code{lwr}, \code{upr} and \code{se}. The
fixed-effect contribution to the standard error comes from the sandwich
\code{vcov}; the random-effect contribution from the partial influence
function of $\hat u$.

<<predict-interval>>=
nd <- data.frame(Days = 0:2, Subject = "308")
suppressWarnings(predict(fm, newdata = nd, interval = "confidence"))
suppressWarnings(predict(fm, newdata = nd, interval = "prediction"))
@

\section{Influence diagnostics}\label{sec:influence}

\subsection{Observation-level influence}\label{sec:cooks-obs}

\code{cooks.distance} measures the joint Mahalanobis influence of each
observation on $(\hat\beta, \hat\sigma, \hat\theta)$, computed from the
full influence function \citep{cook1982residuals}. We illustrate on the
\code{Penicillin} data (a crossed design):

<<cooks-obs, eval=FALSE, echo=FALSE>>=
fpen <- rlmer(diameter ~ (1 | plate) + (1 | sample), Penicillin,
              method = "DASvar")
cd <- cooks.distance(fpen)
sort(cd, decreasing = TRUE)[1:5]
@
<<cooks-obs-output, echo=FALSE, results=tex, strip.white=false>>=
showChunk("cooks-obs")
@

\code{influence(fit)} returns the full per-observation influence matrix
(one row per observation, one column per parameter) that
\code{cooks.distance} condenses; \code{caseweightIF},
\code{vcov\_sandwich} and \code{implicitIF\_full} expose the underlying
substrate for advanced use.

\subsection{Cluster-level influence and robust leverage}\label{sec:cooks-cluster}

A whole group can be an outlier even when no single observation in it
is. \code{cooks.distance(fit, groups = TRUE)} (or
\code{groups = "<factor>"}) gives a per-cluster Cook's distance---the
Mahalanobis norm of the per-cluster influence function---a general
diagnostic for such groups. It is, however, not a reliable screen for
group contamination ahead of the bootstrap variance-component test of
\code{anova} (Section~\ref{sec:anova}): the robust fit absorbs a
contaminated group into a heavily downweighted random effect, so the
bounded influence saturates exactly when contamination is present. That
test's guard therefore uses the random-effect robustness weights, and
fires automatically. Cluster level requires a single or nested grouping
structure.

<<cooks-cluster>>=
cdc <- cooks.distance(fm, groups = "Subject")
sort(cdc, decreasing = TRUE)[1:5]
@

\code{hatvalues(fit)} returns the robust leverage (the self-leverage
$A_{ii}$ in the random-effect-whitened convolution; it reduces to the
classical \pkg{lme4} leverage at \code{rho = cPsi}), with a
\code{groups} argument to sum it per cluster:

<<hatvalues>>=
range(hatvalues(fm))
@

\begin{figure}[t]
\centering
<<cluster-fig, fig=TRUE, echo=FALSE, width=7, height=3.6>>=
op <- par(mar = c(4, 4, 1, 1))
dotchart(sort(cdc), xlab = "Cluster-level Cook's distance",
         ylab = "Subject", pch = 19, cex = 0.7)
par(op)
@
\caption{Per-subject (cluster-level) Cook's distances for the robust
  \code{sleepstudy} fit. Subjects at the top are the most influential
  on the joint parameter vector $(\hat\beta, \hat\sigma, \hat\theta)$.}
\label{fig:cluster}
\end{figure}

\section{A high-breakdown initial estimator: RANSAC}\label{sec:ransac}

By default the RSE estimator uses a monotone (smoothed Huber)
$\psi$-function, which is convex and largely insensitive to the starting
values. For stronger robustness against extreme outliers one can instead
choose a \emph{redescending} $\psi$-function---for example
\code{rho.e = lqqPsi}---whose influence function returns to zero so
that gross outliers are rejected outright. The bisquare
(\code{rho.e = bisquarePsi}) is the classical redescender, but it
redescends comparatively fast; the lqq (linear-quadratic-quadratic) psi
of \citet{ks2011}, the recommended redescender used by
\pkg{robustbase}'s \code{lmrob.control(setting = "KS2014")}, is generally
preferable and is pre-built as \code{lqqPsi} (the general constructor
\code{makeRobustbasePsi} additionally exposes the \code{optimal},
\code{hampel} and \code{ggw} families). Redescending $M$-estimators
are more robust, but their estimating equations are non-convex and can be
drawn to a bad local solution from a poor start. Release 3.4-5 added a
\textsc{ransac} (\textsc{ra}ndom \textsc{sa}mple \textsc{c}onsensus,
\citealp{fischler1981ransac}) high-breakdown start that makes such fits
reliable; 3.5.0 refines it with the adaptive early stopping, stratified
and nonsingular subsampling, and multi-start consensus described below.
The string form \code{init = "ransac"} plugs it straight into
\code{rlmer}:

<<ransac-init, eval=FALSE, echo=FALSE>>=
fr <- rlmer(Reaction ~ Days + (Days | Subject), sleepstudy,
            method = "DASvar", init = "ransac")
fixef(fr)
@
<<ransac-init-output, echo=FALSE, results=tex, strip.white=false>>=
showChunk("ransac-init")
@

\code{ransac\_lme4} exposes the start directly and returns the fitted
subsample, the consensus residual scale and diagnostics. The search is
adaptive: it stops early once the best $Q_n$ score plateaus, so \code{K}
is a maximum budget rather than a fixed number of iterations
(\code{adaptive = FALSE} restores the fixed loop). Subsampling is
stratified by default (\code{stratify = TRUE}): each subsample covers
every level of the finest grouping factor, keeping the per-subsample
\code{lmer} fittable on designs with many small clusters. Following the
nonsingular subsampling of \citet{ks2017nonsingular} (the
\code{setting = "KS2014"} default of \pkg{robustbase}), each draw is
full-rank \emph{by construction}: with categorical predictors a
rank-deficient subsample does not always error---\code{lmer} would
silently drop the aliased columns of a dropped factor level---so, rather
than draw-then-check, an LU decomposition with column pivoting walks a
random permutation of the observations and keeps only those that raise
the fixed-effects rank, skipping any that are collinear with those
already chosen. Whole clusters are then added until the retained data
are a fully identifiable mixed-model subsample (full-rank fixed design,
at least two levels per grouping factor, varying random slopes), which
succeeds whenever the full design is identifiable. The collinear
observations skipped by the LU are counted in \code{n\_singular}; on a
design without collinearity the draw reduces to a uniform random cluster
subsample.

<<ransac-lme4>>=
rs <- ransac_lme4(Reaction ~ Days + (Days | Subject), sleepstudy,
                  K = 20L, seed = 20260601)
rs$K_used
@

\code{rlmer\_ransac} wraps the start and the robust fit and runs three
diagnostics: a fixed-effect \emph{basin} check (warns when a
redescending \code{rho.e} left the half-weight radius
$r^\ast = s_{1/2}\,\sigma / \max_j \lVert x_j \rVert$ of the start,
where $s_{1/2}$ is the residual mean-shift at which the expected
redescender weight $m(s) = \mathrm{E}\,\psi'(Z + s)$, $Z \sim N(0,1)$, has
halved --- $s_{1/2} = 1.410$ for the bisquare, computed from the
\code{Dpsi} slot for lqq and the other redescenders --- so that the
population Hessian stays positive definite on the ball; returned by
\code{ransac\_basin\_radius}), a random-effects phony-correlation
check (warns when the fitted RE covariance is near-singular,
$|\hat\rho| \to 1$), and a refinement-singularity check
(\citealp{ks2017nonsingular}, Remark~2): a redescending $\psi$ can
zero-weight observations back into a rank-deficient design during the
robust fit, so \code{rlmer\_ransac} re-seeds and re-fits from a fresh
nonsingular start until the positive-weight fixed-effects design is full
rank, up to \code{max\_tries} (default 5), warning only if every attempt
still collapses. The multi-start consensus
\code{rlmer\_ransac(n\_starts = m)} fits from the $m$ best distinct
starts and returns the lowest-residual-scale fit whose RE covariance is
interior, recovering the good solution when the single best start falls
into the phony-correlation attractor (per-start summary in
\code{attr(fit, "consensus")}):

<<ransac-consensus, eval=FALSE, echo=FALSE>>=
fc <- rlmer_ransac(Reaction ~ Days + (Days | Subject), sleepstudy,
                   K = 20L, n_starts = 3L, seed = 20260601)
attr(fc, "consensus")
@
<<ransac-consensus-output, echo=FALSE, results=tex, strip.white=false>>=
showChunk("ransac-consensus")
@

The per-start summary reports each start's largest absolute random-effect
correlation (\code{max\_abscor}), its residual scale, whether its
random-effects covariance is interior (not a phony $|\hat\rho| \to 1$
solution), and which start was \code{chosen} (the interior fit with the
lowest residual scale).

What the RANSAC start is for here is raising the breakdown point of the
starting value, and thereby making redescending $\psi$-functions usable
on general designs: their estimating equations are non-convex, so a
low-breakdown start can be drawn to a bad local solution.

The defaults are settled and simulation-backed. Across the
contamination scenarios of
\code{inst/simulationStudy/ransacConsensus.R}, the optional multi-start
consensus (\code{n\_starts > 1}) lowers the pooled phony-correlation rate
from $\Sexpr{fmt(nums$ransacConsensus$phonySingle, 3)}$ to
$\Sexpr{fmt(nums$ransacConsensus$phonyConsensus, 3)}$ at no cost in slope
accuracy (slope RMSE $\Sexpr{fmt(nums$ransacConsensus$slopeRmseSingle)}$
and $\Sexpr{fmt(nums$ransacConsensus$slopeRmseConsensus)}$); over
$\Sexpr{nums$ransacConsensus$nPairs}$ paired fits the consensus is the
interior one in
\Sexpr{with(nums$ransacConsensus, if (singleBetter == 0) paste("all", consensusBetter) else paste(consensusBetter, "of the", consensusBetter + singleBetter))}
replicates where the two disagree. RANSAC handles a
single, crossed or nested grouping structure---subsampling is stratified
by the finest grouping factor so each subsample \code{lmer} stays
fittable on designs with many small clusters---and the start is fully
reproducible from a \code{seed} (or, for the \code{init = "ransac"}
string form, from an outer \code{set.seed}).

\section{Robustness in the design space}\label{sec:design}

\subsection{Mallows-type design weights}\label{sec:designweights}

The $\psi$-functions of the RSE bound the influence of outliers in the
\emph{response}, but a single high-leverage observation---extreme
covariates with a plausible response---can still move $\hat\beta$
arbitrarily. This is the classical gap between Huber-type and
generalised $M$-estimators. The \code{design.weights} argument adds
Mallows-type weights $\eta_i \in (0, 1]$ that downweight the score
contribution of high-leverage design points. The weights enter the
$e$-side of every estimating equation ($\beta, u, \sigma, \theta$), the
model \code{vcov} and the influence diagnostics. Passing \code{"mcd"}
derives $\eta$ from robust Mahalanobis distances of the non-constant
columns of $X$ (via \code{covMcd}; note that with many binary dummy
columns the MCD covariance estimate can be ill-conditioned); a numeric
vector may be supplied instead.

<<design-weights>>=
dat <- sleepstudy
dat$Days[1] <- 40                 # an extreme, high-leverage design point
fw <- rlmer(Reaction ~ Days + (Days | Subject), dat,
            method = "DASvar", design.weights = "mcd")
sum(getME(fw, "design.weights") < 1)
@

\code{summary} reports how many observations are downweighted by
design. On clean data the cost is small---about a
\Sexpr{fmt(nums$mallows$costVsRsePct, 1)}\% increase in the variance of
$\hat\beta$ relative to the plain RSE at the default tuning
($\gamma = 1$, cutoff $0.975$), over \Sexpr{nums$mallows$nreps}
replicates---and it vanishes for intercept-only designs (no
non-constant columns in $X$).
Under leverage contamination that drags ordinary and plain robust fits
toward a spurious slope, the Mallows fit instead stays close to the
clean-data estimate; the efficiency--robustness trade-off is quantified
in \code{inst/simulationStudy/mallowsEfficiency.R} and
\code{leverageBreakdown.R}. Active design weights are supported for a
single grouping factor; \code{rlmer} stops with an error on models with
more than one.

\subsection{Structured random-effect covariances (experimental)}\label{sec:structured}

With \pkg{lme4} $\geq$ 2.0-0, \code{rlmer} supports structured
random-effect covariances: \code{diag(f | g)} (uncorrelated slopes,
equivalent to \code{(f || g)}), \code{cs(f | g)} (compound symmetric:
one common correlation) and \code{ar1(f | g)} (autoregressive,
$\mathrm{Cor}(i, j) = \rho^{|i - j|}$), in both heterogeneous and
homogeneous (equal-variance) forms. \code{diag} is fitted directly;
\code{cs} and \code{ar1} are fitted by projecting each block's
unstructured scoring-equation update onto the structured manifold at
every iteration, so the structure is enforced exactly while the
robustness machinery is unchanged. This projection is a pragmatic
device rather than a fully established estimator, so \code{cs}/\code{ar1}
fitting is \emph{experimental}. It is validated only for a single
grouping factor (balanced repeated measures); \code{rlmer} warns when a
\code{cs}/\code{ar1} term is combined with more than one grouping factor,
where neither the fit nor its inference is validated. Inference degrades
gracefully on structured fits---the sandwich \code{vcov} is available, and
\code{summary(df = "satterthwaite")} reports the robust df where the
constrained covariance admits it (e.g.\ \code{cs}) and otherwise falls
back to the $t$-table with a note.

To show the structure being enforced we simulate a balanced
repeated-measures dataset with \pkg{robustlmm}'s own data generator
\code{generateRepeatedMeasuresDatasets}: a four-level within-subject
factor \code{visit}, replicated over $40$ subjects, with an underlying
AR(1) random-effect covariance ($\sigma = 1.5$, $\rho = 0.6$).

<<structured-data>>=
set.seed(20260623)
datasets <- generateRepeatedMeasuresDatasets(
    1, numberOfSubjects = 40, numberOfVisits = 4, numberOfReplicates = 3,
    structure = "ar1", marginalSd = 1.5, correlation = 0.6,
    trueBeta = 10, trueSigma = 0.7)
d <- datasets$generateData(1)
@

We model the factor with \code{0 + visit} (dropping the intercept) so
that \emph{each} visit gets its own random effect: the term
\code{(0 + visit | subject)} gives a four-dimensional random effect per
subject, and that is the object on which a covariance structure is
defined. An intercept-plus-slope parameterisation would instead give
only a two-dimensional random effect, on which \code{cs} and \code{ar1}
both reduce to a single correlation and so coincide with the
unstructured fit. With four dimensions the structures differ: an
unstructured fit estimates all $\binom{4}{2} = 6$ correlations, the
\code{cs} fit collapses them to one common correlation, and the
\code{ar1} fit to a single decay parameter $\rho$ with
$\mathrm{Cor}(i, j) = \rho^{|i - j|}$. \code{VarCorr} annotates the
structure it enforced:

<<structured-fits, eval=FALSE, echo=FALSE>>=
rcs  <- rlmer(y ~ 1 + cs(0 + visit | subject),  d, method = "DASvar")
rar1 <- rlmer(y ~ 1 + ar1(0 + visit | subject), d, method = "DASvar")
VarCorr(rcs)
VarCorr(rar1)
@
<<structured-fits-output, echo=FALSE, results=tex, strip.white=false>>=
showChunk("structured-fits")
@

Because the projection keeps the marginal variances and only constrains
the correlations, a structured fit resists the over-parameterisation
(and the spurious near-$\pm 1$ correlations) that an unstructured fit
can fall into on short, noisy series, while remaining fully robust to
outliers. A simulation study
(\code{inst/simulationStudy/structuredCovariances.R},
\Sexpr{nums$structured$nreps} replicates) confirms both properties. On
clean data the robust and the classical structured fits both recover the
truth. But when \Sexpr{nums$structured$contamPct}\% of the subjects
carry an outlying random-effect vector, the classical (\code{lmer})
common correlation collapses toward zero---from a true
$\Sexpr{fmt(nums$structured$cs$trueRho)}$ to
$\Sexpr{fmt(nums$structured$cs$contLmer)}$ for \code{cs} and from
$\Sexpr{fmt(nums$structured$ar1$trueRho)}$ to
$\Sexpr{fmt(nums$structured$ar1$contLmer)}$ for \code{ar1}, with its
marginal standard deviations inflating from
$\Sexpr{fmt(nums$structured$cs$trueSd)}$ to
$\Sexpr{fmt(nums$structured$cs$contSdLmer)}$---whereas the robust
structured fit holds near the truth (correlation
$\Sexpr{fmt(nums$structured$cs$contRlmer)}$ for \code{cs} and
$\Sexpr{fmt(nums$structured$ar1$contRlmer)}$ for \code{ar1}, marginal
standard deviation $\Sexpr{fmt(nums$structured$cs$contSdRlmer)}$).

A limitation of the default estimation method carries over to
structured blocks: \code{method = "DAStau"} computes its consistency
factors by Gauss--Hermite quadrature only for random-effect blocks of
dimension $\leq 2$, so for larger blocks---such as the four-dimensional
\code{cs}/\code{ar1} blocks above---the fit falls back to \code{DASvar}
with a warning (which is why the fits above request \code{DASvar}
directly). The experimental option
\code{options(robustlmm.dastau.mc = TRUE)} enables a Monte-Carlo
calibration of the same DAS-tau fixed point for such blocks; the sample
is drawn once per fit from a deterministic seed, so fits remain
reproducible. In a companion simulation on clean Gaussian data this
removes the small calibration residual that the \code{DASvar}
approximation leaves in the fitted correlation: for \code{ar1} blocks of
dimension $\Sexpr{nums$dasmc$blockDim}$ with $\Sexpr{nums$dasmc$J}$
subjects the mean difference to \code{lmer} was at most
$\Sexpr{fmt(nums$dasmc$dasvarResidual, 3)}$ in absolute value with
\code{DASvar} and $\Sexpr{fmt(nums$dasmc$dasmcResidual, 3)}$ with the
Monte-Carlo calibration
(\code{inst/simulationStudy/dasmcValidation.R}).
It is validated by simulation only---not by a finite-sample
theorem---and for blocks of dimension $\leq 2$ the quadrature path
remains in place, as the Monte-Carlo path offers no improvement there.
See \code{?"robustlmm-options"} for the option's companions and
caveats.

\section{A robust sensitivity analysis for an MMRM}\label{sec:mmrm}

In clinical trials with repeated measurements the primary analysis is
often a mixed model for repeated measures (MMRM): treatment, visit and
their interaction as fixed effects, and an unstructured covariance of
the visits within a subject. A robust fit can serve as a sensitivity
analysis next to it: if a few unusual observations or subjects drive the
result, the robust estimate moves away from the classical one, and the
robustness weights show where. This section goes through such an analysis
with the example data \code{fev\_data} of the \pkg{mmrm} package
\citep{mmrm}: \Sexpr{mm$weights$rlmer$nobs} lung-function measurements
(\code{FEV1}) of \Sexpr{mm$weights$rlmer$nsubj} subjects at four visits,
in a treatment and a placebo arm. It uses \code{augment()} and
\code{contrast\_test()}, which are new in version 3.5.1. The code is
shown but not run here, so
that the vignette does not need \pkg{mmrm}; the results were computed
with \pkg{mmrm} \Sexpr{mm$versions[["mmrm"]]} by
\code{vignettes/precompute-mmrm.R}.

\subsection{The models differ}

The MMRM specifies the marginal covariance of the four visits directly
and leaves it unstructured:
<<mmrm-fit, eval=FALSE>>=
library("mmrm")
fitMmrm <- mmrm(FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
                data = fev_data)
@
\code{rlmer} cannot fit this model. Its covariance comes from random
effects plus an independent residual error, and with one observation per
subject and visit a visit-level random effect cannot be separated from
the residual error. An unstructured block over the visits is therefore
not identifiable, and the structured alternative \code{cs(0 + AVISIT |
USUBJID)} is refused for the same reason (\code{rlmer} asks for more than
four observations per subject). What can be fitted are random-effects
models for the subject: a random intercept, which implies the same
variance at every visit and the same correlation between any two, and a
random intercept and slope in the visit number, which lets the variance
change smoothly over the visits.
<<mmrm-rlmer, eval=FALSE>>=
fitRob <- rlmer(FEV1 ~ RACE + SEX + ARMCD * AVISIT + (1 | USUBJID),
                data = fev_data)
fev_data$visit <- as.numeric(fev_data$AVISIT) - 2.5
fitRobSlope <- rlmer(FEV1 ~ RACE + SEX + ARMCD * AVISIT +
                         (visit | USUBJID),
                     data = fev_data, method = "DASvar")
@
The robust analysis is thus a check under a different, more restrictive
covariance model, not a robust version of the MMRM.

\subsection{Treatment difference at the last visit}

For the MMRM and the random-intercept fit the treatment difference at the
last visit comes from \pkg{emmeans}, with Satterthwaite degrees of
freedom in both cases:
<<mmrm-emmeans, eval=FALSE>>=
library("emmeans")
contrast(emmeans(fitMmrm, ~ ARMCD | AVISIT), "revpairwise")
contrast(emmeans(fitRob, ~ ARMCD | AVISIT), "revpairwise")
@
The same contrast can be written in terms of the coefficients and tested
with \code{contrast\_test()}, which also offers the cluster sandwich:
<<mmrm-contrast, eval=FALSE>>=
L <- c("ARMCDTRT" = 1, "ARMCDTRT:AVISITVIS4" = 1)
contrast_test(fitRob, L)
contrast_test(fitRob, L, vcov_type = "sandwich")
@

<<mmrm-table, echo=FALSE, results=tex>>=
ctab <- mm$contrasts
lab <- c(mmrm = "MMRM, unstructured covariance",
         lmer = "\\code{lmer}, random intercept",
         rlmer = "\\code{rlmer}, random intercept",
         rlmerSandwich = "\\quad with the cluster sandwich",
         slope = "\\code{rlmer}, random intercept and slope",
         slopeSandwich = "\\quad with the cluster sandwich")
cat("\\begin{table}[htb]\n\\centering\n",
    "\\begin{tabular}{lrrrr}\n\\toprule \n",
    "Fit & Estimate & SE & df & $p$ \\\\\n\\midrule \n", sep = "")
for (r in rownames(ctab))
    cat(lab[[r]], " & ", fmt(ctab[r, "estimate"]), " & ",
        fmt(ctab[r, "se"]), " & ",
        if (is.finite(ctab[r, "df"])) fmt(ctab[r, "df"], 0) else "--",
        " & ", if (ctab[r, "p"] < 0.001) "$<$0.001" else fmt(ctab[r, "p"], 3),
        " \\\\\n", sep = "")
cat("\\bottomrule \n\\end{tabular}\n",
    "\\caption[Treatment difference at the last visit]{Treatment ",
    "difference in FEV1 at the last visit. The standard errors without ",
    "the sandwich are model-based; ``--'' is the normal reference.}\n",
    "\\label{tab:mmrm}\n\\end{table}\n", sep = "")
@

Table~\ref{tab:mmrm} shows the results. The robust estimates,
$\Sexpr{fmt(mm$contrasts["rlmer", "estimate"])}$ and
$\Sexpr{fmt(mm$contrasts["slope", "estimate"])}$, are somewhat smaller
than the MMRM's $\Sexpr{fmt(mm$contrasts["mmrm", "estimate"])}$, by less
than one of its standard errors. The standard errors differ more than
the estimates. The model-based standard errors of the random-effects fits,
classical and robust, are about
$\Sexpr{fmt(mm$contrasts["lmer", "se"], 1)}$, against
$\Sexpr{fmt(mm$contrasts["mmrm", "se"], 1)}$ for the MMRM: they rely on
a common variance for all visits, while the MMRM estimates a standard
deviation of $\Sexpr{fmt(mm$variance$mmrmSd[["VIS4"]], 1)}$ at the last
visit and $\Sexpr{fmt(mm$variance$mmrmSd[["VIS3"]], 1)}$ at the third.
The cluster sandwich does not rely on the covariance model, and with it
the robust fits give standard errors of
$\Sexpr{fmt(mm$contrasts["rlmerSandwich", "se"], 1)}$ and
$\Sexpr{fmt(mm$contrasts["slopeSandwich", "se"], 1)}$, close to the
MMRM's; the $p$-value of the random-intercept fit is then
$\Sexpr{fmt(mm$contrasts["rlmerSandwich", "p"], 3)}$. When the covariance
model of the robust fit is cruder than the primary model's, as here, the
sandwich is the standard error to compare.

\subsection{Where the fits disagree}

The robustness weights come with the data from \code{augment()}:
<<mmrm-augment, eval=FALSE>>=
aug <- generics::augment(fitRob)
aug[order(aug$.wt_e)[1:6], ]
table(aug$AVISIT[aug$.wt_e < 0.8])
@
In the random-intercept fit, \Sexpr{mm$weights$rlmer$nLow} observations
have a weight below $0.5$, and of the
\Sexpr{sum(mm$weights$rlmer$belowByVisit)} observations with a weight
below $0.8$, \Sexpr{mm$weights$rlmer$belowByVisit[["VIS4"]]} are from the
last visit, which has a quarter of the observations. The downweighted
observations are thus not scattered recording errors: they sit where the
response is more variable than the model allows. The same shows in the
variance components. The robust random-intercept fit estimates a subject
standard deviation of $\Sexpr{fmt(mm$variance$rlmerSd[["subject"]], 1)}$
(\code{lmer}: $\Sexpr{fmt(mm$variance$lmerSd[["subject"]], 1)}$), a
boundary fit that says the common-correlation model does not describe
these data. With the random slope the picture is calmer:
\Sexpr{mm$weights$slope$nLow} observations below $0.5$, and
\Sexpr{mm$weights$slope$belowByVisit[["VIS4"]]} of
\Sexpr{sum(mm$weights$slope$belowByVisit)} below $0.8$ at the last
visit. Here the weights point to a shortcoming of the model, the
unequal variances of the visits, which the unstructured MMRM accommodates
by construction.

The subjects that drive a fit are found with the cluster-level Cook's
distance of Section~\ref{sec:cooks-cluster}:
<<mmrm-cooks, eval=FALSE>>=
cd <- cooks.distance(fitRobSlope, groups = TRUE)
head(sort(cd, decreasing = TRUE))
@
In the random-slope fit the largest value is
$\Sexpr{fmt(mm$cooks$slope$max, 1)}$ (subject
\Sexpr{names(mm$cooks$slope$top)[1]}) against a median of
$\Sexpr{fmt(mm$cooks$slope$median, 1)}$, and the next ones follow
closely: no single subject dominates. The subject with the smallest
random-effect weight is \Sexpr{mm$weights$slope$minWtBSubject}
($\Sexpr{fmt(mm$weights$slope$minWtB)}$); the number of subjects with a
weight below one is \Sexpr{mm$weights$slope$nSubjBelow1}.
This computation needs the influence function of the fit and takes
several minutes for a model of this size; \code{rlmer} announces it with
a message.

\subsection{What this analysis does and does not show}

For these data the sensitivity analysis gives no sign that outliers
drive the primary result: the robust estimates are within one standard
error of the MMRM's, no subject dominates, and the downweighted
observations are explained by the larger variance at the last visit. Three limits should be stated with any such analysis.
\begin{itemize}
\item The models differ. The MMRM has an unstructured marginal covariance
  of the visits; \code{rlmer} has random effects and an independent
  error. A difference between the two fits can come from the covariance
  model as well as from outliers.
\item The MMRM with an unstructured covariance over the visits cannot be
  fitted by \code{rlmer}, because that block is not identifiable next to
  the residual error.
\item Missing visits: the MMRM is valid when the dropout is missing at
  random. Whether the robust estimating equations are, has not been
  studied. With dropout that depends on the observed responses the robust
  fit should be read as descriptive.
\end{itemize}

\section{Summary}\label{sec:summary}

Release 3.5.0 turns \pkg{robustlmm} from an estimator into a more
complete robust-modelling toolbox: robust standard errors and
confidence intervals, ANOVA tables, prediction intervals and honest
degrees of freedom for inference; observation- and cluster-level
influence diagnostics; a high-breakdown \textsc{ransac} start; and
design weights and structured covariances for robustness in the design
and covariance space. For the underlying estimator and its tuning, see
the main vignette \citep{vignette}; the simulation studies that
validate these features live in \code{inst/simulationStudy}.

<<sessionInfo, echo=FALSE, results=hide>>=
## (session info intentionally not printed to keep the vignette stable)
@

\bibliography{rlmer}

\end{document}