R/rlmer.R

In robustlmm: Robust Linear Mixed Effects Models

##' Robust estimation of linear mixed effects models, for hierarchical nested
##' and non-nested, e.g., crossed, datasets.
##'
##' \describe{ \item{Overview:}{
##'
##' This function implements the Robust Scoring Equations estimator for linear
##' mixed effect models. It can be used much like the function
##' \code{\link[lme4]{lmer}} in the package \code{lme4}. The supported models
##' are the same as for \code{\link[lme4]{lmer}} (gaussian family only). The
##' robust approach used is based on the robustification of the scoring
##' equations and an application of the Design Adaptive Scale approach.
##'
##' Example analyses and theoretical details on the method are available in the
##' vignette (see \code{vignette("rlmer")}).
##'
##' Models are specified using the \code{formula} argument, using the same
##' syntax as for \code{\link[lme4]{lmer}}. Additionally, one also needs to
##' specify what robust scoring or weight functions are to be used (arguments
##' starting with \code{rho.}). By default a smoothed version of the Huber
##' function is used. Furthermore, the \code{method} argument can be used to
##' speed up computations at the expense of accuracy of the results. }
##'
##' \item{Computation methods:}{
##'
##' Currently, there are two different methods available for fitting models.
##' They only differ in how the consistency factors for the Design Adaptive
##' Scale estimates are computed. Available fitting methods for theta and
##' sigma.e: \itemize{
##'
##' \item \code{DAStau} (default): For this method, the consistency factors are
##' computed using numerical quadrature. This is slower but yields more accurate
##' results. This is the direct analogue to the DAS-estimate in robust linear
##' regression.
##'
##' \item \code{DASvar}: This method computes the consistency factors using a
##' direct approximation which is faster but less accurate. For complex models
##' with correlated random effects with more than one correlation term, this is
##' the only method available.
##'
##' }
##'
##' \code{DAStau} supports blocks of random effects of dimension at most 2;
##' for fits containing a larger block it falls back to \code{DASvar} with a
##' warning. An \emph{experimental} option,
##' \code{options(robustlmm.dastau.mc = TRUE)}, lifts this restriction via a
##' Monte-Carlo calibration of the consistency factors; it is
##' simulation-validated but not backed by a finite-sample theorem, see
##' \code{\link{robustlmm-options}}. }
##'
##' \item{Weight functions:}{
##'
##' The tuning parameters of the weight functions \dQuote{rho} can be used to
##' adjust robustness and efficiency of the resulting estimates (arguments
##' \code{rho.e}, \code{rho.b}, \code{rho.sigma.e} and \code{rho.sigma.b}).
##' Better robustness will lead to a decrease of the efficiency. With the default
##' setting, \code{setting = "RSEn"}, the tuning parameters are set to yield
##' estimates with approximately 95\% efficiency for the fixed effects. The
##' variance components are estimated with a lower efficiency but better
##' robustness properties.
##'
##' One has to use different weight functions and tuning parameters for simple
##' variance components and for such including correlation parameters. By
##' default, they are chosen appropriately to the model at hand. However, when
##' using the \code{rho.sigma.e} and \code{rho.sigma.b} arguments, it is up to
##' the user to specify the appropriate function. See
##' \code{\link{asymptoticEfficiency}} for methods to find tuning parameters
##' that yield a given asymptotic efficiency. \itemize{
##'
##' \item For simple variance components and the residual error scale use the
##' function \code{\link{psi2propII}} to change the tuning parameters. This is
##' similar to Proposal 2 in the location-scale problem (i.e., using the
##' squared robustness weights of the location estimate for the scale estimate;
##' otherwise the scale estimate is not robust).
##'
##' \item For multi-dimensional blocks of random effects modeled, e.g.,
##' a model with correlated random intercept and slope, (referred to as
##' block diagonal case below), use the \code{\link{chgDefaults}} function to
##' change the tuning parameters. The parameter estimation problem is
##' multivariate, unlike the case without correlation where the problem was
##' univariate. For the employed estimator, this amounts to switching from
##' simple scale estimates to estimating correlation matrices. Therefore
##' different weight functions have to be used. Squaring of the weights (using
##' the function \code{\link{psi2propII}}) is no longer necessary. To yield
##' estimates with the same efficiency, the tuning parameters for the
##' block diagonal are larger than for the simple case. Tables of tuning parameters
##' are given in Table 2 and 3 of the vignette (\code{vignette("rlmer")}).
##'
##' } }
##'
##' \item{Recommended tuning parameters:}{
##'
##' For a more robust estimate, use \code{setting = "RSEn"} (the default). For
##' higher efficiency, use \code{setting = "RSEa"}. The settings described in
##' the following paragraph are used when \code{setting = "RSEa"} is specified.
##'
##' For the smoothed Huber function the tuning parameters to get approximately
##' 95\% efficiency are \eqn{k=1.345}{k=1.345} for \code{rho.e} and
##' \eqn{k=2.28}{k=2.28} for \code{rho.sigma.e} (using the squared version). For
##' simple variance components, the same can be used for \code{rho.b} and
##' \code{rho.sigma.b}. For variance components including correlation
##' parameters, use \eqn{k=5.14}{k=5.14} for both \code{rho.b} and
##' \code{rho.sigma.b}. Tables of tuning parameter are given in Table 2 and 3 of
##' the vignette (\code{vignette("rlmer")}). }
##'
##' \item{Specifying (multiple) weight functions:}{
##'
##' If custom weight functions are specified using the argument \code{rho.b}
##' (\code{rho.e}) but the argument \code{rho.sigma.b} (\code{rho.sigma.e}) is
##' missing, then the squared weights are used for simple variance components
##' and the regular weights are used for variance components including
##' correlation parameters. The same tuning parameters will be used when
##' \code{setting = "RSEn"} is used. To get
##' higher efficiency either use \code{setting = "RSEa"} (and only set arguments
##' \code{rho.e} and \code{rho.b}). Or specify the tuning parameters by hand
##' using the \code{\link{psi2propII}} and \code{\link{chgDefaults}} functions.
##'
##' To specify separate weight functions \code{rho.b} and \code{rho.sigma.b} for
##' different variance components, it is possible to pass a list instead of a
##' psi_func object. The list entries correspond to the groups as shown by
##' \code{VarCorr(.)} when applied to the model fitted with \code{lmer}. A set
##' of correlated random effects count as just one group. }
##'
##' \item{\code{lmerNoFit}:}{
##'
##' The \code{lmerNoFit} function can be used to get trivial starting values.
##' This is mainly used to verify the algorithms to reproduce the fit by
##' \code{\link[lme4]{lmer}} when starting from trivial initial values. } }
##'
##' @section Prior weights: \code{weights} is passed through to
##'   \code{\link[lme4]{lmer}} and has the usual meaning,
##'   \eqn{\mathrm{Var}(e_i) = \sigma^2 / w_i}. \strong{Weighted
##'   \code{method = "DASvar"} fits made with robustlmm 3.2-0 to 3.4-5
##'   are wrong and should be refitted}: the Design Adaptive Scale
##'   counted the prior weights twice, so \eqn{\hat\sigma}, the
##'   variance components and every standard error derived from them
##'   were biased by an amount that depended on the scale the weights
##'   happened to be supplied on. The default \code{method = "DAStau"}
##'   was affected only through its starting value, and in practice
##'   only at \code{rho.e = rho.b = cPsi}. Unweighted fits are
##'   unaffected in every version. See the 3.5.0 entry of
##'   \code{news(package = "robustlmm")} for the magnitudes and the
##'   measured exposure of each method.
##' @title Robust Scoring Equations Estimator for Linear Mixed Models
##' @param formula a two-sided linear formula object describing the
##'   fixed-effects part of the model, with the response on the left of a
##'   \code{~} operator and the terms, separated by \code{+} operators, on the
##'   right.  The vertical bar character \code{"|"} separates an expression for
##'   a model matrix and a grouping factor.
##' @param data an optional data frame containing the variables named in
##'   \code{formula}.  By default the variables are taken from the environment
##'   from which \code{lmer} is called.
##' @param ... Additional parameters passed to lmer to find the initial
##'   estimates. See \code{\link[lme4]{lmer}}.
##' @param method method to be used for estimation of theta and sigma, see
##'   Details.
##' @param setting a string specifying suggested choices for the arguments
##'   \code{rho.e}, \code{rho.sigma.e}, \code{rho.b} and \code{rho.sigma.b}.
##'   Use \code{"RSEn"} (the default) or \code{"RSEa"}. Both use
##'   \code{\link{smoothPsi}} for all the \dQuote{rho} arguments. For
##'   \code{rho.sigma.e}, squared robustness weights are used (see
##'   \code{\link{psi2propII}}). \code{"RSEn"} uses the same tuning parameter as
##'   for \code{rho.e}, which leads to higher robustness but lower efficiency.
##'   \code{"RSEa"} adjusts the tuning parameter for higher asymptotic efficiency
##'   which results in lower robustness (\code{k = 2.28} for default \code{rho.e}).
##'   For diagonal random effects covariance matrices, \code{rho.sigma.b} is
##'   treated exactly as \code{rho.sigma.e}. For block diagonal random effects
##'   covariance matrices (with correlation terms), regular robustness weights
##'   are used for \code{rho.sigma.b}, not squared ones, as they're not needed.
##'   But the tuning parameters are adjusted for both \code{rho.b} and
##'   \code{rho.sigma.b} according to the dimensions of the blocks (for both
##'   \code{"RSEn"} or \code{"RSEa"}). For a block of dimension 2 (e.g.,
##'   correlated random intercept and slope) \code{k = 5.14} is used.
##' @param rho.e object of class psi_func, specifying the functions to use for
##'   the huberization of the residuals.
##' @param rho.b object of class psi_func or list of such objects (see Details),
##'   specifying the functions to use for the huberization of the random
##'   effects.
##' @param rho.sigma.e object of class psi_func, specifying the weight functions
##'   to use for the huberization of the residuals when estimating the variance
##'   components, use the \code{\link{psi2propII}} function to specify squared
##'   weights and custom tuning parameters.
##' @param rho.sigma.b (optional) object of class psi_func or list of such
##'   objects, specifying the weight functions to use for the huberization of
##'   the random effects when estimating the variance components (see Details).
##'   Use \code{\link{psi2propII}} to specify squared weights and custom tuning
##'   parameters or \code{\link{chgDefaults}} for regular weights for variance
##'   components including correlation parameters.
##' @param rel.tol relative tolerance used as criteria in the fitting process.
##' @param max.iter maximum number of iterations allowed.
##' @param verbose verbosity of output. Ranges from 0 (none) to 3 (a lot of
##'   output)
##' @param doFit logical scalar. When \code{doFit = FALSE} the model is not fit
##'   but instead a structure with the model matrices for the random-effects
##'   terms is returned (used to speed up tests). When \code{doFit = TRUE}, the
##'   default, the model is fit immediately.
##' @param init optional lmerMod- or rlmerMod-object to use for starting values,
##'   a list with elements \sQuote{fixef}, \sQuote{u}, \sQuote{sigma},
##'   \sQuote{theta}, the string \code{"ransac"}, or a function producing an
##'   lmerMod object.
##'
##'   \strong{An lmerMod or rlmerMod passed here supplies the model, not only
##'   the starting values}: the working object is built from it, so
##'   \code{formula} and \code{data} no longer have any effect. Pass such an
##'   object only when it was fitted on the same data; a warning is issued if
##'   its number of observations does not match \code{nrow(data)}. To start
##'   from values obtained elsewhere while keeping the model given by
##'   \code{formula} and \code{data}, use the list form. Its components are
##'   optional --- whatever is omitted keeps the value from the internal
##'   \code{lmer} initialisation, which is the right choice for \sQuote{u}
##'   when the values come from a fit on different observations.
##'
##'   When \code{init = "ransac"}, the high-breakdown RANSAC
##'   start is obtained by calling \code{\link{ransac_lme4}(formula, data)} with
##'   its default \code{K = 200} subsamples and \code{sub_frac = 0.5}. Its
##'   purpose is to raise the breakdown point of the starting value, which is
##'   what makes redescending psi-functions usable: their estimating equations
##'   are non-convex, so a low-breakdown start can be drawn to a bad local
##'   solution. \code{\link{rlmer}} takes
##'   \sQuote{fixef}, \sQuote{sigma} and \sQuote{theta} from that fit and fits
##'   the model on the full \code{data}. The string form is random
##'   (not reproducible without an outer \code{set.seed}); for fine control of
##'   \code{K}, \code{sub_frac}, or \code{seed}, use \code{\link{rlmer_ransac}}.
##' @param size_obr logical scalar; if \code{TRUE} (default \code{FALSE}),
##'   the size-controlling weight \code{w_delta} in the block-diagonal
##'   variance-components scoring equation is replaced by the
##'   Hampel-OBR (Stahel 1987 / Hampel et al. 1986) form
##'   \eqn{w_\tau(d) = \min(1, b_\tau / |d - s\kappa - a|)}, with
##'   \eqn{b_\tau} taken from the tuning constant of
##'   \code{rho.sigma.b} and \eqn{a} determined by Fisher consistency
##'   under \eqn{\chi^2_s}. The default finite-difference form
##'   \eqn{w_\delta(d) = (\psi(d) - \psi(d-s\kappa))/s} is not
##'   Hampel-OBR-optimal at the central model; this option recovers
##'   asymptotic Hampel-OBR efficiency for the variance-component
##'   magnitudes (a typically modest 1-2 percentage-point gain at
##'   matched gross-error sensitivity). Has no effect for diagonal
##'   \eqn{V_b} (block size 1).
##' @param design.weights Mallows-type design weights for robustness to
##'   high-leverage design points: \code{NULL} (default, exact
##'   current behaviour), a numeric vector of deterministic weights
##'   \eqn{\eta_i \in (0, 1]} of length \code{n}, or \code{"mcd"} to
##'   compute \eqn{\eta_i = \min(1, \chi^2_{p^*}(0.975) / d_i^2)^{1/2}}
##'   from robust squared Mahalanobis distances \eqn{d_i^2} of the
##'   non-constant fixed-effects covariates
##'   (\code{\link[robustbase]{covMcd}}); the \code{"mcd"} tuning
##'   (severity \eqn{\gamma = 1}, cutoff \eqn{c = 0.975}) is the
##'   simulation-backed default (see \code{vignette("rlmer")}). The
##'   weights multiply the
##'   e-side score contribution of each observation throughout the
##'   estimator -- the \eqn{\beta}, \eqn{u}, \eqn{\sigma} and
##'   \eqn{\theta} equations, the model \code{vcov}, the Satterthwaite
##'   degrees of freedom and the influence
##'   diagnostics -- bounding the influence of high-leverage design
##'   points. For designs without continuous covariates \code{"mcd"}
##'   reduces to \eqn{\eta \equiv 1} (the unmodified estimator). Active
##'   design weights are supported for a \emph{single grouping factor}
##'   only; \code{rlmer} stops with an error if the model has more than
##'   one. Leverage robustness costs a little efficiency at the clean model
##'   (about 1\% at the default tuning; see \code{vignette("rlmer")}).
##' @return object of class rlmerMod.
##' @seealso \code{\link[lme4]{lmer}}, \code{vignette("rlmer")}
##' @author Manuel Koller, with thanks to Vanda Lourenço for improvements.
##' @keywords models
##' @examples
##' ## dropping of VC
##' system.time(print(rlmer(Yield ~ (1|Batch), Dyestuff2, method="DASvar")))
##'
##' \dontrun{
##'   ## Default method "DAStau"
##'   system.time(rfm.DAStau <- rlmer(Yield ~ (1|Batch), Dyestuff))
##'   summary(rfm.DAStau)
##'   ## DASvar method (faster, less accurate)
##'   system.time(rfm.DASvar <- rlmer(Yield ~ (1|Batch), Dyestuff,
##'                                   method="DASvar"))
##'   ## compare the two
##'   compare(rfm.DAStau, rfm.DASvar)
##'
##'   ## Fit variance components with higher efficiency
##'   ## psi2propII yields squared weights to get robust estimates
##'   ## this is the same as using rlmer's argument `setting = "RSEa"`
##'   rlmer(diameter ~ 1 + (1|plate) + (1|sample), Penicillin,
##'         rho.sigma.e = psi2propII(smoothPsi, k = 2.28),
##'         rho.sigma.b = psi2propII(smoothPsi, k = 2.28))
##'
##'   ## use chgDefaults for variance components including
##'   ## correlation terms (regular, non squared weights suffice)
##'   ## this is the same as using rlmer's argument `setting = "RSEa"`
##'   rlmer(Reaction ~ Days + (Days|Subject), sleepstudy,
##'         rho.sigma.e = psi2propII(smoothPsi, k = 2.28),
##'         rho.b = chgDefaults(smoothPsi, k = 5.14, s=10),
##'         rho.sigma.b = chgDefaults(smoothPsi, k = 5.14, s=10))
##' }
##'
##' @importFrom lme4 lmer
##' @importFrom stats getCall
##' @export
rlmer <- function(formula, data, ..., method = c("DAStau", "DASvar"),
                  setting, rho.e, rho.b, rho.sigma.e, rho.sigma.b,
                  rel.tol = 1e-8, max.iter = 40 * (r + 1)^2, verbose = 0,
                  doFit = TRUE, init, size_obr = FALSE,
                  design.weights = NULL)
{
    lcall <- match.call()
    pf <- parent.frame()
    method <- match.arg(method)
    lobj <- .rlmerInit(lcall, pf, formula, data, method, rho.e, rho.b, rho.sigma.e,
                       rho.sigma.b, rel.tol, max.iter, verbose, init, setting, ...)$obj
    lobj@resp$eta <- .processDesignWeights(design.weights, lobj)
    lobj@pp <- as(lobj@pp, "rlmerPredD_DAS")
    lobj@pp$method <- lobj@method
    lobj@pp$initRho(lobj)
    lobj@pp$initMatrices(lobj)

    ## required for max.iter:
    r <- len(lobj, "theta")

    return(.rlmer(lobj, rel.tol, max.iter, verbose, doFit,
                  size_obr = isTRUE(size_obr)))
}

## WS11 step 1: validate / compute the Mallows design weights eta.
## Returns a length-n numeric in (0, 1]. "mcd" follows
## PLAN-WS11-mallows.md section 2.2 with gamma = 1, cutoff 0.975:
## eta_i = min(1, chi^2_{p*}(0.975) / d_i^2)^{1/2}, d_i^2 the robust
## squared Mahalanobis distance of the non-constant columns of X.
.processDesignWeights <- function(design.weights, lobj) {
    n <- nrow(lobj@pp$X)
    if (is.null(design.weights))
        return(rep.int(1, n))
    if (is.character(design.weights)) {
        design.weights <- match.arg(design.weights, "mcd")
        X <- as.matrix(lobj@pp$X)
        nonconst <- apply(X, 2L, function(col) length(unique(col)) > 1L)
        Xs <- X[, nonconst, drop = FALSE]
        pstar <- ncol(Xs)
        if (pstar == 0L) {
            ## no continuous design information: exact no-op
            eta <- rep.int(1, n)
        } else if (n < 2L * (pstar + 1L)) {
            warning("design.weights = \"mcd\": n = ", n,
                    " too small for covMcd with ", pstar,
                    " non-constant column(s); using eta = 1.")
            eta <- rep.int(1, n)
        } else {
            mcd <- robustbase::covMcd(Xs)
            d2  <- stats::mahalanobis(Xs, mcd$center, mcd$cov)
            eta <- pmin(1, sqrt(stats::qchisq(0.975, df = pstar) /
                                pmax(d2, .Machine$double.eps)))
        }
    } else {
        eta <- as.numeric(design.weights)
        if (length(eta) != n)
            stop("design.weights must have length ", n,
                 " (one weight per observation)")
        if (any(!is.finite(eta)) || any(eta <= 0) || any(eta > 1))
            stop("design.weights must be finite and in (0, 1]")
    }
    ## Mallows design weights are supported for a single grouping factor
    ## only: the e-side weighting flows into the beta/u/sigma/theta
    ## equations, the model vcov, the Satterthwaite df and the influence
    ## diagnostics, all validated for one grouping factor (WS11; see
    ## vignette("rlmer")). Fail loudly outside that scope -- when weights
    ## are actually active -- rather than return a silently unvalidated
    ## fit. (A no-op request, eta == 1 throughout, is left alone.)
    if (any(eta < 1)) {
        ngf <- length(lobj@flist)
        if (ngf > 1L)
            stop("design.weights (Mallows eta) currently supports a single ",
                 "grouping factor, but this model has ", ngf, " (",
                 paste(names(lobj@flist), collapse = ", "),
                 "). Fit without design.weights, or reduce the model to one ",
                 "grouping factor.", call. = FALSE)
    }
    eta
}

## Covariance structures (lme4 >= 2.0-0) that robustlmm can fit.
## us/diag are estimated directly; cs/ar1 are estimated by projecting each
## block's unstructured covariance update onto the structured manifold
## (see .structuredCovInfo() and rlmer.fit.DAS.nondiag()). Both the
## heterogeneous and homogeneous (equal-variance) forms are supported.
.supportedCovClasses <- c("Covariance.us", "Covariance.diag",
                          "Covariance.cs", "Covariance.ar1")

## Covariance structures handled by the projection step.
.projectedCovClasses <- c("Covariance.cs", "Covariance.ar1")

## Check for structured covariance matrices (lme4 >= 2.0-0)
.checkStructuredCovariance <- function(object) {
    if (!lme4::anyStructured(object)) {
        return(invisible(NULL))  # No structured covariances
    }

    # Check each random effect's covariance structure
    reCovs <- attr(object, "reCovs")
    if (is.null(reCovs)) {
        return(invisible(NULL))  # No reCovs attribute (shouldn't happen)
    }

    for (i in seq_along(reCovs)) {
        covClass <- class(reCovs[[i]])[1]
        if (!covClass %in% .supportedCovClasses) {
            structName <- sub("Covariance\\.", "", covClass)
            stop("robustlmm does not yet support '", structName,
                 "()' covariance structure (from lme4 >= 2.0-0). ",
                 "Supported: (f | g), (f || g), diag(f | g), ",
                 "cs(f | g), ar1(f | g).",
                 call. = FALSE)
        }
    }

    ## cs / ar1 are fitted by an iterative projection onto the structured
    ## manifold, validated only for a single grouping factor (balanced
    ## repeated measures); the fit AND its inference are untested with more
    ## than one grouping factor (crossed/nested). Warn rather than silently
    ## return an unvalidated fit. diag is fitted directly and is exempt.
    isProj <- vapply(reCovs, function(rc)
        class(rc)[1L] %in% .projectedCovClasses, logical(1))
    if (any(isProj) && length(object@flist) > 1L)
        warning("Structured cs()/ar1() covariances are experimental and ",
                "validated only for a single grouping factor; this model has ",
                length(object@flist), " (",
                paste(names(object@flist), collapse = ", "), "). The ",
                "projected fit and its inference are not validated for more ",
                "than one grouping factor; interpret the results with care.",
                call. = FALSE)

    ## A projected cs()/ar1() block of random-effect dimension nc is only
    ## identifiable when every cluster carries MORE than nc observations: with
    ## at most nc observations per cluster there is no within-cluster residual
    ## degree of freedom to separate the structured random effect from the
    ## residual, and the robust DAS correction matrices go singular (the fit
    ## then dies deep inside the projection with a cryptic error). Detect this
    ## "saturated" case from the model frame and stop with a clear message.
    ## reCovs[[i]] aligns positionally with object@cnms / the random-effect
    ## terms, so names(object@cnms)[i] is the grouping factor for term i.
    gfNames <- names(object@cnms)
    for (i in seq_along(reCovs)) {
        if (!isProj[i])
            next
        nc <- reCovs[[i]]@nc
        gf <- gfNames[i]
        clusterSizes <- tabulate(object@flist[[gf]])
        minObs <- if (length(clusterSizes)) min(clusterSizes) else 0L
        if (minObs <= nc) {
            structName <- sub("Covariance\\.", "", class(reCovs[[i]])[1L])
            stop("Structured ", structName, "() covariance with random-effect ",
                 "dimension nc=", nc, " requires more than nc observations per ",
                 "cluster of '", gf, "' (the smallest cluster has ", minObs,
                 "); the structured fit is not identifiable without ",
                 "within-cluster replication. Add replication or use a smaller ",
                 "structured dimension.", call. = FALSE)
        }
    }
    invisible(NULL)
}

## Starting values taken from a fit on a DIFFERENT set of observations
## (a RANSAC subsample), in the form `.rlmerInit()`'s list branch wants.
##
## `fixef`, `sigma` and `theta` have lengths fixed by the model, not by
## the number of observations, so they carry over.  `u` does not -- it
## has one entry per random-effect coefficient of the fit it came from --
## and it is deliberately omitted rather than zeroed: see the comment in
## `.rlmerInit()`.  The high-breakdown start is for beta anyway.
.ransacInitList <- function(fit) {
    list(fixef = unname(fixef(fit)),
         sigma = sigma(fit),
         theta = unname(getME(fit, "theta")))
}

.rlmerInit <- function(lcall, pf, formula, data, method, rho.e, rho.b, rho.sigma.e,
                       rho.sigma.b, rel.tol, max.iter, verbose, init,
                       setting = c("RSEn", "RSEa"), ...) {
    ## Resolve init = "ransac" to a PARAMETER LIST before the branches
    ## below.  It must not resolve to the merMod that ransac_lme4()
    ## returns: an merMod passed as `init` supplies the model as well as
    ## the starting values -- `lobj <- as(init, "rlmerMod")` further down
    ## takes the model frame from it and `formula`/`data` are ignored
    ## from that point on -- and the RANSAC fit lives on the winning
    ## SUBSAMPLE.  Resolving to the fit object therefore silently
    ## estimated the model on roughly half the data (robustlmm 3.4-5 and
    ## 3.5.0: 510 of 1000 observations, 27 of 50 clusters, on the
    ## two-dimensional-block design of the regression test).
    if (!missing(init) && !is.null(init) && identical(init, "ransac")) {
        rr <- ransac_lme4(formula, data)
        if (is.null(rr$fit))
            stop("init = \"ransac\": no successful lmer fit across the RANSAC subsamples")
        init <- .ransacInitList(rr$fit)
    }
    if (missing(init) || is.null(init) || is.list(init)) {
        lcall2 <- lcall
        lcall2[setdiff(names(formals(rlmer)), names(formals(lmer)))] <- NULL
        lcall2$doFit <- NULL
        lcall2$REML <- TRUE
        lcall2[[1]] <- substitute(lme4::lmer)
        linit <- eval(lcall2, pf)
        .checkStructuredCovariance(linit)
        if (!missing(init) && is.list(init)) {
            ## Components may be omitted; whatever is absent keeps the
            ## value the full-data lmer initialisation produced.  When
            ## `u` is omitted but `theta` is given -- the RANSAC case,
            ## where `u` has the wrong length and cannot be carried over
            ## -- the effects are re-solved instead, so the starting
            ## point is self-consistent rather than a mixture of
            ## supplied parameters and an lmer BLUP fitted under
            ## different ones.  Note what that solve is: fitEffects()
            ## is the robust IRLS step for (beta, u) JOINTLY, run to
            ## convergence, so the supplied `fixef` is a warm start that
            ## gets re-solved, not an anchor that is held.  Only sigma
            ## and theta are held at the supplied values.
            ##
            ## Requiring the names to be exactly right is not pedantry:
            ## since the components became optional, `list(fixed = ...)`
            ## or `list()` would otherwise fall through silently to the
            ## plain lmer start and return a fit that looks fine.  It
            ## also removes the inconsistency of `$` partial matching,
            ## under which `fixeff` resolved to `fixef` but `fixed` did
            ## not.  Extraction below uses [[ ]], which matches exactly.
            valid <- c("fixef", "u", "sigma", "theta")
            nms <- names(init)
            if (!length(init) || is.null(nms) || !all(nzchar(nms)) ||
                !all(nms %in% valid) || anyDuplicated(nms))
                stop("'init' list components must be named ",
                     "'fixef', 'u', 'sigma' and/or 'theta'", call. = FALSE)
            ## check sanity of whatever was supplied
            if (!is.null(init[["fixef"]]))
                stopifnot(length(init[["fixef"]]) == length(fixef(linit)))
            if (!is.null(init[["u"]]))
                stopifnot(length(init[["u"]]) == length(getME(linit, "u")))
            if (!is.null(init[["sigma"]]))
                stopifnot(length(init[["sigma"]]) == length(sigma(linit)))
            if (!is.null(init[["theta"]]))
                stopifnot(length(init[["theta"]]) == length(getME(linit, "theta")))
            ## convert object to rlmerMod
            linit <- as(linit, "rlmerMod")
            ## set the supplied initial parameters, but do not fit yet
            ## Order is load-bearing, not stylistic: the effects solve
            ## triggered below reads beta and b from the object at call
            ## time and its robustness weights depend on sigma, so beta
            ## and sigma must already be set.  Reordering would start
            ## the solve from stale values and fail silently.
            if (!is.null(init[["fixef"]])) setFixef(linit, unname(init[["fixef"]]))
            if (!is.null(init[["sigma"]])) setSigma(linit, init[["sigma"]])
            if (!is.null(init[["theta"]])) {
                ## fit.effects = TRUE re-solves (beta, u) at the sigma
                ## and theta just set; update.sigma = FALSE keeps the
                ## supplied sigma.  Only do this when `u` was not
                ## supplied, so a fully specified init is used verbatim.
                refitU <- is.null(init[["u"]])
                setTheta(linit, init[["theta"]], fit.effects = refitU,
                         update.sigma = FALSE)
            }
            if (!is.null(init[["u"]])) setU(linit, init[["u"]])
        }
        init <- linit
    } else if (is.function(init)) {
        init <- do.call(init,list(formula=formula, data=data, REML=TRUE, ...))
    } else if (!is(init, "merMod") && !is(init, "rlmerMod")) {
        stop("Unsuitable init object, aborting.",
                "Expecting no, list (see ?rlmer), rlmerMod or merMod object")
    }
    ## Check for structured covariance matrices when init is a merMod
    ## (the check for internally-fitted models is above)
    if (is(init, "merMod") && !is(init, "rlmerMod")) {
        .checkStructuredCovariance(init)
    }
    ## An merMod passed as `init` becomes the model: the conversion below
    ## takes the model frame from it, so `formula` and `data` no longer
    ## have any effect.  That is intended for a starting fit on the same
    ## data, and silently wrong for one fitted on a subset.  Warn rather
    ## than stop, because dropping rows with NAs is a legitimate reason
    ## for the counts to differ; pass starting VALUES through the list
    ## form of `init` if the model should come from `formula`/`data`.
    if ((is(init, "merMod") || is(init, "rlmerMod")) &&
        !missing(data) && is.data.frame(data)) {
        ninit <- tryCatch(nobs(init), error = function(e) NA_integer_)
        ## Compare against the rows the model would actually use, not
        ## nrow(data): rows dropped for NAs are a legitimate reason for
        ## the counts to differ and must not raise a false alarm.
        ##
        ## Counted from the model variables rather than via
        ## model.frame(lme4::subbars(formula)).  `subbars` is exported by
        ## every lme4 this package supports, but it now lives in
        ## `reformulas` and lme4 emits a deprecation warning on every
        ## call -- one that tells the user to ask the upstream
        ## maintainer to stop.  A heuristic behind a warning is not worth
        ## a warning of its own.
        ndata <- tryCatch({
            vars <- intersect(all.vars(formula), names(data))
            if (!length(vars)) NA_integer_
            else sum(stats::complete.cases(data[, vars, drop = FALSE]))
        }, error = function(e) NA_integer_)
        if (!is.na(ninit) && !is.na(ndata) && ninit != ndata)
            warning(gettextf(
                paste0("'init' was fitted on %d observations but 'formula' ",
                       "and 'data' give %d; the fit will use the model from ",
                       "'init', not from 'formula' and 'data'. If that is ",
                       "not what you want (e.g. 'init' is a subsample fit), ",
                       "pass starting values as a list, see ?rlmer."),
                ninit, ndata), call. = FALSE)
    }
    lobj <- as(init, "rlmerMod")
    lobj@call <- lcall

    if (any(lobj@resp$weights == 0))
        stop("Observations with zero weights are not allowed.")
    if (any(lobj@resp$weights < 0))
        stop("Observations with negative weights are not allowed.")
    tooSmall <- lobj@resp$weights < 0.1
    if (any(tooSmall)) {
        obsString <- createObservationsString(tooSmall)
        warn(paste0("Detected very small weights (", obsString, "). ",
                    "You may encounter issues due to numerical instability."))
    }

    if (!missing(setting)) {
        if (!missing(rho.e) || !missing(rho.sigma.e) ||
            !missing(rho.b) || !missing(rho.sigma.b)) {
            overridden <- c()
            if (!missing(rho.e)) {
                overridden <- c(overridden, "'rho.e'")
            }
            if (!missing(rho.sigma.e)) {
                overridden <- c(overridden, "'rho.sigma.e'")
            }
            if (!missing(rho.b)) {
                overridden <- c(overridden, "'rho.b'")
            }
            if (!missing(rho.sigma.b)) {
                overridden <- c(overridden, "'rho.sigma.b'")
            }
            if (length(overridden) > 1) {
                args <- "Arguments "
            } else {
                args <- "Argument "
            }
            overridden <- paste(overridden, collapse = ", ")
            setting <- match.arg(setting)
            warning("Argument 'setting' specified together with ",
                    overridden, ". ", args, overridden,
                    " will override defaults of \"", setting, "\".")
        }
        setting <- match.arg(setting)
    } else {
        setting <- "RSEn"
    }
    if (missing(rho.e)) {
        rho.e <- smoothPsi
    }
    adjustForEfficiency <- setting == "RSEa"
    if (missing(rho.sigma.e)) {
        rho.sigma.e <- psi2propII(rho.e, adjust = adjustForEfficiency)
    }
    if (missing(rho.b)) {
        rho.b <- lapply(lobj@dim, getDefaultRhoB, rho = rho.e)
    } else if (!is.list(rho.b)) {
        rho.b <- rep.int(list(rho.b), length(lobj@dim))
        ## TODO warn if low asymptotic efficiency?
    }
    if (missing(rho.sigma.b)) {
        rho.sigma.b <- list()
        for (bt in seq_along(lobj@dim)) {
            if (lobj@dim[bt] == 1) {
                rho.sigma.b[[bt]] <- psi2propII(rho.b[[bt]], adjust = adjustForEfficiency)
            } else {
                rho.sigma.b[[bt]] <- rho.b[[bt]]
            }
        }
    } else if (!is.list(rho.sigma.b)) {
        rho.sigma.b <- rep.int(list(rho.sigma.b), length(lobj@dim))
        ## TODO warn if low asymptotic efficiency?
    }
    ## set arguments only relevant to rlmerMod
    lobj@rho.b <- rho.b
    lobj@rho.sigma.b <- rho.sigma.b
    if (!isTRUE(chk <- validObject(lobj))) stop(chk)
    lobj@rho.e <- rho.e
    lobj@rho.sigma.e <- rho.sigma.e
    if (method == "DAStau" & any(sapply(lobj@idx, nrow) > 2)) {
        if (.dasMcEnabled()) {
            ## experimental Monte-Carlo DAS-tau calibration: it handles
            ## any block dimension, so do NOT fall back to DASvar (see
            ## calcTau.nondiag / .buildDasMcSamples)
        } else {
            warning("Method 'DAStau' does not support blocks of size larger than 2. ",
                    "Falling back to method 'DASvar'.")
            method <- "DASvar"
        }
    }
    lobj@method <- method
    return(list(obj = lobj, init = init))
}

getDefaultRhoB <- function(dimension, rho) {
    if (dimension == 1) {
        return(rho)
    }
    if (isDefaultHuberOrSmoothPsi(rho) && dimension < 8) {
        k <- switch(dimension - 1,
                    5.14, 5.55, 5.91, 6.25, 6.55, 6.84)
    } else {
        k <-
            findTuningParameter(asymptoticEfficiency(rho, "location"),
                                rho,
                                "tau",
                                dimension)
    }
    return(chgDefaults(rho, k = k))
}

.rlmer <- function(lobj, rel.tol, max.iter, verbose, doFit,
                   size_obr = FALSE) {
    if (!doFit) {
        lobj@pp$updateMatrices()
        return(updateWeights(lobj))
    }

    ## do not start with theta == 0
    if (any(theta(lobj)[lobj@lower == 0] == 0)) {
        if (verbose > 0)
            cat("Setting variance components from 0 to 1\n")
        theta0 <- theta(lobj)
        theta0[lobj@lower == 0 & theta0 == 0] <- 1
        setTheta(lobj, theta0, fit.effects = TRUE, update.sigma = TRUE)
    } else {
        ## set theta at least once
        setTheta(lobj, theta(lobj), fit.effects = FALSE)
    }

    if (verbose > 0) {
        cat("\nrlmer starting values:\n")
        cat("sigma, theta: ", .sigma(lobj), ", ", theta(lobj), "\n")
        cat("coef: ", .fixef(lobj), "\n")
        if (verbose > 1)
            cat("b.s: ", b.s(lobj), "\n")
    }

    curWarnings <- list()
    lobj <- withCallingHandlers(
        .rlmer.fit(lobj, rel.tol, max.iter, verbose, size_obr = size_obr),
        warning = function(w) {
            curWarnings <<- append(curWarnings, list(conditionMessage(w)))
            invokeRestart("muffleWarning")
        })
    lobj@optinfo$warnings <- curWarnings
    lobj@optinfo$size_obr <- size_obr

    if (verbose > 0) {
        cat("sigma, theta: ", .sigma(lobj), ", ", theta(lobj), "\n")
        cat("coef: ", .fixef(lobj), "\n")
        if (verbose > 1)
            cat("b.s: ", b.s(lobj), "\n")
    }

    return(updateWeights(lobj))
}

.rlmer.fit <- function(lobj, rel.tol, max.iter, verbose,
                       size_obr = FALSE) {
    ## do fit: non diagonal case differently
    if (!isDiagonal(.U_b(lobj))) {
        lobj <- rlmer.fit.DAS.nondiag(lobj, verbose, max.iter, rel.tol,
                                      size_obr = size_obr)
    } else {
        lobj <- rlmer.fit.DAS(lobj, verbose, max.iter, rel.tol)
    }
    return(lobj)
}

## Length of the theta (relative Cholesky) vector for a Covariance
## object. Reimplements lme4's getThetaLength() from the documented
## structure (class, nc, hom) so robustlmm does not depend on the
## unexported lme4 generic (lme4 < 2.0-2 does not export it).
.covThetaLength <- function(reCov) {
    nc <- reCov@nc
    cls <- class(reCov)[1L]
    if (cls == "Covariance.diag") {
        nc
    } else if (cls %in% .projectedCovClasses && isTRUE(reCov@hom)) {
        2L * nc - (nc > 0L)          # homogeneous cs / ar1
    } else {
        nc * (nc + 1L) / 2L          # us, or heterogeneous cs / ar1
    }
}

## Map structured parameters par = (sigma..., rho) to lme4's relative
## Cholesky theta. This reimplements lme4's (unexported) getTheta() for
## the cs and ar1 families -- the formulas are the analytic Cholesky of
## the structured correlation matrix, identical across lme4 versions, so
## the result is byte-compatible with what VarCorr() reconstructs from
## the Covariance object's par slot (which we set via initialize()).
.covGetTheta <- function(reCov, par) {
    nc <- reCov@nc
    if (nc <= 1L) {
        return(par)
    }
    hom <- isTRUE(reCov@hom)
    rho <- par[length(par)]
    rho2 <- rho * rho
    if (class(reCov)[1L] == "Covariance.ar1") {
        v1 <- rho^(0L:(nc - 1L))
        v2 <- rho^(0L:(nc - 2L)) * sqrt(1 - rho2)
        if (hom) {
            par[1L] * c(v1, v2)
        } else {
            par[sequence.default(from = 1L:nc, nvec = nc:1L)] *
                c(v1, v2[sequence.default(from = 1L, nvec = (nc - 1L):1L)])
        }
    } else {                         # cs
        a <- double(nc)
        a. <- 0
        for (j in 1L:nc) {
            a[j] <- a.
            a. <- a. + (1 - rho * a.)^2 / (1 - rho2 * a.)
        }
        v. <- sqrt(1 - rho2 * a)
        v <- rbind(v., (rho - rho2 * a) / v.)
        if (hom) {
            par[1L] * v[1L:(2L * nc - 1L)]
        } else {
            par[sequence.default(from = 1L:nc, nvec = nc:1L)] *
                rep(v, rbind(1L, (nc - 1L):0L))
        }
    }
}

## Structured parameters par = (sigma..., rho) of the nearest cs / ar1
## covariance to a full nc x nc covariance Sigma: marginal variances are
## kept (per-component for heterogeneous, pooled for homogeneous) and the
## single correlation parameter rho is estimated from the sample
## correlation matrix Cor = D^{-1} Sigma D^{-1} (D = diag of marginal sds).
##
## cs: rho is the average off-diagonal correlation -- which is also the
## constrained estimator for the (linear) compound-symmetric structure
## (Anderson 1973), so the projected fixed point is the constrained
## estimating-equation solution.
##
## ar1: rho is the conditional maximiser of the working-Gaussian objective
## over the AR(1) manifold given the marginal variances -- i.e. the
## minimiser of the Stein/entropy loss tr(R(rho)^{-1} Cor) + (nc-1)
## log(1 - rho^2), where R(rho)_{ij} = rho^|i-j|. This is the ECME / CM
## step (Meng & Rubin 1993; Liu & Rubin 1994): it makes the projected
## fixed point a stationary point of the constrained Gaussian objective,
## unlike the earlier ad-hoc log-correlation least-squares fit (which is
## not such a maximiser and cannot represent rho < 0). The 1-D problem is
## solved by optimize() on (-1, 1).
##
## Returns NULL on a degenerate covariance (non-positive marginal variance).
.covParFromCov <- function(reCov, Sigma) {
    dvar <- diag(Sigma)
    if (any(dvar <= 0)) {
        return(NULL)
    }
    sdv <- sqrt(dvar)
    Cor <- Sigma / outer(sdv, sdv)
    if (class(reCov)[1L] == "Covariance.ar1") {
        nc  <- nrow(Sigma)
        lag <- abs(outer(seq_len(nc), seq_len(nc), "-"))
        ## negative working-Gaussian log-likelihood profiled over rho
        ## (Stein loss up to the rho-free constant), to be minimised:
        steinLoss <- function(rho) {
            R   <- rho ^ lag
            sol <- tryCatch(sum(diag(solve(R, Cor))),
                            error = function(e) NA_real_)
            if (!is.finite(sol)) return(.Machine$double.xmax)
            sol + (nc - 1) * log1p(-rho * rho)
        }
        opt <- stats::optimize(steinLoss, interval = c(-0.999, 0.999))
        rho <- opt$minimum
    } else {
        rho <- mean(Cor[lower.tri(Cor)])
    }
    sigma <- if (isTRUE(reCov@hom)) sqrt(mean(dvar)) else sdv
    c(sigma, rho)
}

## Lower bounds for a full nc x nc Cholesky factor in lme4's column-major
## lower-triangular theta order: diagonal entries >= 0, off-diagonal
## entries unconstrained. Matches getME(., "lower") for unstructured (us)
## terms and is the correct theta-space bound for cs/ar1 as well (lme4's
## getME(., "lower") returns the par-space bounds for those, which do not
## align with robustlmm's Cholesky-space theta).
.usThetaLower <- function(nc) {
    unlist(lapply(seq_len(nc), function(j) c(0, rep(-Inf, nc - j))))
}

## Build the theta-space (Cholesky) lower bounds for a converted model.
## For old lme4 (no reCovs) the lme4 bounds are already theta-space and
## are returned unchanged.
.structuredThetaLower <- function(from) {
    reCovs <- attr(from, "reCovs")
    deflt <- getME(from, "lower")
    if (is.null(reCovs)) {
        return(deflt)
    }
    ## Only cs/ar1 need a rebuild (their getME() bounds are par-space).
    ## For us/diag the lme4 bounds are already correct theta-space bounds
    ## and carry the theta names downstream code expects, so return them
    ## unchanged -- this keeps non-structured fits byte-identical.
    isProjected <- vapply(reCovs,
                          function(rc) class(rc)[1L] %in% .projectedCovClasses,
                          logical(1))
    if (!any(isProjected)) {
        return(deflt)
    }
    out <- numeric(0)
    for (rc in reCovs) {
        cls <- class(rc)[1L]
        if (cls == "Covariance.diag") {
            ## theta == diagonal variances, all >= 0
            out <- c(out, rep(0, .covThetaLength(rc)))
        } else if (cls %in% .projectedCovClasses) {
            ## cs / ar1. theta is a Cholesky-derived vector; the par-space
            ## bounds from getME() do not align with it. Constrain the
            ## diagonal (Cholesky) entries to be non-negative to keep the
            ## block covariance positive-definite during iteration, leaving
            ## off-diagonal entries free. Heterogeneous structures store a
            ## full lower triangle; homogeneous ones a reduced vector whose
            ## leading entry is the (non-negative) scale.
            nt <- .covThetaLength(rc)
            nc <- rc@nc
            if (nt == nc * (nc + 1L) / 2L) {
                out <- c(out, .usThetaLower(nc))
            } else {
                out <- c(out, c(0, rep(-Inf, nt - 1L)))
            }
        } else {
            ## us: full lower-triangular Cholesky
            out <- c(out, .usThetaLower(rc@nc))
        }
    }
    out
}

## Map robustlmm fitting blocks to their lme4 Covariance object, for the
## structured families (cs, ar1) that need a projection step.
##
## Returns a list aligned with lobj@blocks: entry [[type]] is the
## Covariance object when block `type` belongs to a cs/ar1 random effect,
## and NULL otherwise (us/diag need no projection -- their natural
## parametrisation already coincides with the unstructured Cholesky).
##
## The mapping is by theta index: each Covariance owns a contiguous range
## of the global theta vector (in cnms order), and each block's theta
## indices fall entirely inside exactly one such range. This is robust to
## diag() terms, which expand into several 1x1 blocks.
.structuredCovInfo <- function(lobj) {
    info <- vector("list", length(lobj@blocks))
    reCovs <- attr(lobj, "reCovs")
    if (is.null(reCovs)) {
        return(info)
    }
    ## contiguous theta range owned by each Covariance object
    start <- 1L
    ranges <- vector("list", length(reCovs))
    for (i in seq_along(reCovs)) {
        n <- .covThetaLength(reCovs[[i]])
        ranges[[i]] <- start:(start + n - 1L)
        start <- start + n
    }
    for (type in seq_along(lobj@blocks)) {
        Ub <- lobj@blocks[[type]]
        lind <- Ub[Ub != 0]
        for (i in seq_along(reCovs)) {
            if (all(lind %in% ranges[[i]])) {
                if (class(reCovs[[i]])[1L] %in% .projectedCovClasses) {
                    info[[type]] <- reCovs[[i]]
                }
                break
            }
        }
    }
    info
}

## Project an unstructured block covariance Sigma (= relative covariance
## factor times its transpose) onto the structured manifold (cs / ar1)
## defined by `reCov`, returning the structured theta vector.
##
## The block's scoring-equation update produces a full nc x nc Sigma. We
## retract it onto the structured family by keeping the marginal variances
## and estimating the structure's single correlation parameter rho from
## the sample correlations:
##   cs : rho is the average off-diagonal correlation;
##   ar1: rho is the Stein-loss / ECME profile-rho minimiser over the
##        AR(1) manifold given the marginal variances (see .covParFromCov;
##        this replaced the earlier log-correlation least-squares fit,
##        which could not represent rho < 0).
## Variances become per-component (heterogeneous, sigma_i = sqrt(Sigma_ii))
## or pooled (homogeneous, sigma = sqrt(mean diag)) to match the family's
## par layout (sigma..., rho). getTheta() then rebuilds the structured
## theta, so VarCorr() stays consistent. This is exact (identity) for a
## Sigma already on the manifold.
##
## We deliberately do NOT use lme4's setTheta() retraction: it reads rho
## from a single Cholesky column and distorts the marginal variances off
## the manifold.
##
## Returns NULL on the rare degenerate case (a non-positive marginal
## variance, or a theta-length mismatch) so the caller can keep the
## previous block estimate and continue iterating.
.projectCovToStructure <- function(reCov, Sigma, ntheta) {
    par <- .covParFromCov(reCov, Sigma)
    if (is.null(par)) {
        return(NULL)
    }
    proj <- .covGetTheta(reCov, par)
    if (length(proj) != ntheta || anyNA(proj)) {
        return(NULL)
    }
    proj
}

## DAS method
rlmer.fit.DAS.nondiag <- function(lobj, verbose, max.iter, rel.tol, method=lobj@method,
                                  checkFalseConvergence = TRUE,
                                  size_obr = FALSE) {
    if (!.isREML(lobj))
        stop("can only do REML when using averaged DAS-estimate for sigma")

    ## per-block projection onto cs/ar1 manifolds (NULL for us/diag)
    structInfo <- .structuredCovInfo(lobj)

    ## Prepare for DAStau
    if (method == "DAStau") {
        ## 4d int
        ## vectorize it!
        ghZ <- as.matrix(expand.grid(lobj@pp$ghz, lobj@pp$ghz, lobj@pp$ghz, lobj@pp$ghz))
        ghZ12 <- ghZ[, 1:2]
        ghZ34 <- ghZ[, 3:4]
        ghw <- apply(as.matrix(expand.grid(lobj@pp$ghw, lobj@pp$ghw, lobj@pp$ghw, lobj@pp$ghw)), 1, prod)
        ## experimental Monte-Carlo DAS-tau calibration: draw the
        ## common-random-numbers sample ONCE per fit; NULL if disabled
        mcSamples <- .buildDasMcSamples(lobj)
        if (!is.null(mcSamples)) {
            mcTypes <- which(!vapply(mcSamples, is.null, logical(1)))
            message("Using EXPERIMENTAL Monte-Carlo DAS-tau calibration ",
                    "(robustlmm.dastau.mc) for block type(s) ",
                    paste(mcTypes, collapse = ", "), " [dim ",
                    paste(lobj@dim[mcTypes], collapse = ", "), "] with N = ",
                    attr(mcSamples, "nsim"), " draws.")
        }
    } else {
        ghZ <- ghw <- c()
        mcSamples <- NULL
    }

    ## fit model using EM algorithm
    conv <- FALSE
    convBlks <- rep(FALSE, length(lobj@blocks))
    iter <- 0
    rel.tol <- sqrt(rel.tol)
    ## compute kappa
    kappas <- .kappa_b(lobj)
    ## zero pattern for T matrix
    nzT <- as.matrix(crossprod(.bdiag(lobj@blocks[lobj@ind])) == 0)
    expectedEntriesLength <- prod(dim(nzT)) - sum(nzT)
    q <- lobj@pp$q
    ## false convergence indicator
    fc <- rep(FALSE, length(lobj@blocks))
    if (verbose > 0) {
        theta0 <- theta(lobj)
        if (verbose > 1) {
            coef0 <- .fixef(lobj)
            b.s0 <- b.s(lobj)
            sigma0 <- .sigma(lobj)
        }
    }
    ## iterate
    while(!conv && (iter <- iter + 1) < max.iter) {
        if (verbose > 0) cat("---- Iteration", iter, " ----\n")
        thetatilde <- theta(lobj)
        sigma <- .sigma(lobj)

        ## get expected value of cov(\tvbs)
        q <- len(lobj, "b")
        T <- switch(method,
                    DASvar=lobj@pp$Tb(),
                    DAStau=calcTau.nondiag(lobj, ghZ12, ghZ34, ghw, .S(lobj), kappas, max.iter,
                                           rel.tol = rel.tol, verbose = verbose,
                                           mcSamples = mcSamples),
                    stop("Non-diagonal case only implemented for DASvar and DAStau"))
        T <- as(T, "CsparseMatrix")
        ## compute robustness weights and add to t and bs
        T[nzT] <- 0
        ## apply chol to non-zero part only
        idx <- !.zeroB(lobj)
        ## stop if all are zero
        if (!any(idx)) break
        Tidx <- T[idx, idx]
        if (any(diag(Tidx) == 0.0)) {
            ## partially dropped block: can't handle yet.
            idxZeroes <- which(idx)[which(diag(Tidx) == 0.0)]
            blockAffected <- which(sapply(lobj@idx, function(cand) any(cand == idxZeroes)))
            stop("Covariance matrix for random effects block ", blockAffected,
                 " with grouping factor ", names(lobj@cnms)[blockAffected],
                 " is singular. Please simplify the model and run again.")
        }
        L <- t(chol(Tidx))
        T.bs <- numeric(q) ## set the others to zero
        T.bs[idx] <- forwardsolve(L, b.s(lobj)[idx])
        ## compute weights
        db <- .dk(lobj, sigma, FALSE, T.bs)[lobj@k]
        wbsEta <- wbsDelta <- numeric(q)
        for (type in seq_along(lobj@blocks)) {
            s <- lobj@dim[type]
            lidx <- as.vector(lobj@idx[[type]])
            if (s > 1) {
                ## for eta, we would actually need a little smaller
                ## tuning constants than for delta to get the same efficiency
                wbsEta[lidx] <- lobj@rho.sigma.b[[type]]@wgt(db[lidx])
                if (size_obr) {
                    ## Hampel-OBR size weight (Stahel 1987 / Hampel et al
                    ## 1986): w_tau(d) = min(1, b_tau / |d - s*kappa - a|).
                    ## b_tau comes from the rho.sigma.b tuning constant; a
                    ## is solved for Fisher consistency under chi^2_s.
                    b_tau <- .obrBTau(lobj@rho.sigma.b[[type]])
                    a_obr <- .stahelObrA(s, kappas[type], b_tau)
                    wbsDelta[lidx] <- pmin(
                        1, b_tau / pmax(abs(db[lidx] - s * kappas[type] - a_obr),
                                        .Machine$double.eps))
                } else {
                    wbsDelta[lidx] <- (lobj@rho.sigma.b[[type]]@psi(db[lidx]) -
                                           lobj@rho.sigma.b[[type]]@psi(db[lidx] - s*kappas[type]))/s
                }
            } else {
                lw <- lobj@rho.sigma.b[[type]]@wgt(db[lidx])
                wbsEta[lidx] <- lw
                wbsDelta[lidx] <- lw*kappas[type] ## adding kappa to wbsDelta in 1d case
            }
        }
        WbDelta <- Diagonal(x=wbsDelta)
        T <- as(WbDelta %*% T, "CsparseMatrix")
        bs <- sqrt(wbsEta) * b.s(lobj)
        Tstart <- 1

        ## cycle block types
        for(type in seq_along(lobj@blocks)) {
            bidx <- lobj@idx[[type]]
            s <- nrow(bidx)
            K <- ncol(bidx)
            Tend <- Tstart + s * s * K - 1
            if (convBlks[type]) {
                Tstart <- Tend + 1
                next
            }
            if (verbose > 5) {
                cat("Tau for blocktype ", type, ":", as.vector(T[bidx[,1],bidx[,1]]), "\n")
            }
            ## catch dropped vc
            if (all(abs(bs[bidx]) < 1e-7)) {
                if (verbose > 1)
                    cat("Block", type, "dropped (all = 0), stopping iterations.\n")
                Ubtilde <- lobj@blocks[[type]]
                pat <- Ubtilde != 0
                Lind <- Ubtilde[pat]
                thetatilde[Lind] <- 0
                convBlks[type] <- TRUE
                Tstart <- Tend + 1
                next
            }
            ## right hand side
            if (length(T@x) == expectedEntriesLength) {
                rhs <- matrix(rowSums(matrix(T@x[Tstart:Tend], s * s)), s)
                if (isTRUE(getOption("robustlmm.check_rhs_optimisation"))) {
                    check <- matrix(0, s, s)
                    for (k in 1:K) {
                        check <- check + as.matrix(T[bidx[,k],bidx[,k]])
                    }
                    stopifnot(all.equal(rhs, check))
                }
            } else {
                rhs <- matrix(0, s, s)
                for (k in 1:K) {
                    rhs <- rhs + as.matrix(T[bidx[,k],bidx[,k]])
                }
            }
            rhs <- rhs / K
            Tstart <- Tend + 1
            ## left hand side
            lbs <- matrix(bs[bidx], ncol(bidx), nrow(bidx), byrow=TRUE)
            ## add left hand side
            lhs <- crossprod(lbs / sigma) / K
            if (verbose > 2) {
                cat("LHS:", as.vector(lhs), "\n")
                cat("RHS:", as.vector(rhs), "\n")
                cat("sum(abs(LHS - RHS)):", sum(abs(lhs - rhs)), "\n")
            }
            diff <- abs(rhs - lhs)
            if (all(diff < rel.tol * max(diff, rel.tol))) {
                if (verbose > 1)
                    cat("Estimating equations satisfied for block", type,
                        ", stopping iterations.\n")
                convBlks[type] <- TRUE
                next
            }
            deltaT <- backsolve(lchol(rhs), lchol(lhs))
            if (verbose > 1) cat("deltaT:", c(deltaT), "\n")
            ## get old parameter estimates for this block
            Ubtilde <- lobj@blocks[[type]]
            pat <- Ubtilde != 0
            Lind <- Ubtilde[pat]
            diagLind <- diag(Ubtilde)
            Ubtilde[pat] <- thetatilde[Lind]
            ## new (unstructured) relative covariance factor for this block
            Unew <- tcrossprod(Ubtilde, deltaT)
            if (is.null(structInfo[[type]])) {
                ## unstructured / diagonal: update theta directly
                thetatilde[Lind] <- Unew[pat]
            } else {
                ## structured covariance (cs / ar1): project the block's
                ## covariance onto the structured manifold before the
                ## convergence check, so the algorithm converges to the
                ## structured fixed point. We project the covariance (not
                ## the Cholesky theta) because homogeneous structures reuse
                ## theta entries across Cholesky positions, so theta is not
                ## a free lower triangle. The structured theta is written
                ## to the block's contiguous theta range (lme4's canonical
                ## order). On a degenerate projection keep the prior block.
                rng <- min(Lind):max(Lind)
                proj <- .projectCovToStructure(structInfo[[type]],
                                               tcrossprod(Unew), length(rng))
                if (!is.null(proj)) {
                    thetatilde[rng] <- proj
                }
            }
            ## FIXME: check boundary conditions?
            ## check if varcomp is dropped
            if (all(thetatilde[diagLind] < 1e-7)) {
                thetatilde[Lind] <- 0
                convBlks[type] <- TRUE
                next
            }
            ## check if this block is converged
            diff <- abs(thetatilde[Lind] - theta(lobj)[Lind])
            if (verbose > 3)
                cat("criterion:", sum(diff), ">=",
                    rel.tol * max(diff, rel.tol), ":",
                    sum(diff) < rel.tol * max(diff, rel.tol), "\n")
            if (sum(diff) < rel.tol * max(diff, rel.tol)) {
                convBlks[type] <- TRUE
                ## check if estimating equations are satisfied. Skipped for
                ## structured (cs/ar1) blocks: the test compares the full
                ## unstructured score (lhs vs rhs), which cannot vanish
                ## within the constrained family, so it would always fire a
                ## spurious false-convergence warning.
                if (checkFalseConvergence && is.null(structInfo[[type]])) {
                    if (verbose > 3)
                        cat("checking estimating equations:", sum(abs(lhs - rhs)),
                            ">", sqrt(rel.tol), ":", sum(abs(lhs - rhs)) > sqrt(rel.tol), "\n")
                    if (sum(abs(lhs - rhs)) > sqrt(rel.tol))
                        fc[type] <- TRUE
                }
                next
            }
        }
        ## set theta
        setTheta(lobj, thetatilde, fit.effects = TRUE,
                 update.sigma = FALSE)
        ## update sigma without refitting effects
        updateSigma(lobj, fit.effects = FALSE)
        if (verbose > 0) {
            cat("delta theta:", format(theta0 - thetatilde, nsmall=20, scientific=FALSE),
                "\n")
            theta0 <- thetatilde
            if (verbose > 1) {
                cat(sprintf("delta coef:  %.12f\n", sum(abs(coef0 - .fixef(lobj)))))
                cat(sprintf("delta u:     %.12f\n", sum(abs(b.s0 - b.s(lobj)))))
                cat(sprintf("delta sigma: %.12f\n", abs(sigma0 - .sigma(lobj))))
                coef0 <- .fixef(lobj)
                b.s0 <- b.s(lobj)
                sigma0 <- .sigma(lobj)
                if (verbose > 2) {
                    cat("theta:", format(thetatilde, nsmall=20, scientific=FALSE), "\n")
                    cat("coef:   ", .fixef(lobj),"\n")
                    cat("b.s:    ", b.s(lobj), "\n")
                    cat("sigmae: ", .sigma(lobj), "\n")
                }
            }
        }
        if (all(convBlks)) conv <- TRUE
    }

    optinfo <- list(optimizer = "rlmer.fit.DAS.nondiag",
                    conv = list(opt = 0),
                    feval = iter,
                    warnings = list(),
                    val = diff)

    if (iter == max.iter) {
        warning(wt <- "iterations did not converge, returning unconverged estimate.")
        optinfo$warnings <- list(wt)
        optinfo$conv$opt <- 1
    }
    if (any(fc)) {
        warning(wt <- "algorithm converged, but estimating equations are not satisfied.")
        optinfo$warnings <- c(optinfo$warnings, list(wt))
        optinfo$conv$opt <- 2
    }
    lobj@optinfo <- optinfo
    return(lobj)
}

## DAS method
rlmer.fit.DAS <- function(lobj, verbose, max.iter, rel.tol) {
    if (!.isREML(lobj))
        stop("can only do REML when using averaged DAS-estimate for sigma")

    ## fit
    converged <- FALSE
    theta0 <- theta(lobj)
    if (verbose > 1) {
        coef0 <- .fixef(lobj)
        b.s0 <- b.s(lobj)
        sigma0 <- .sigma(lobj)
    }
    iter <- 0
    while (!converged && iter < max.iter) {
        iter <- iter + 1
        if (verbose > 0) cat("Iteration", iter, "\n")
        ## fit theta
        updateThetaTau(lobj, max.iter, rel.tol/10, verbose)
        theta1 <- theta(lobj)

        if (verbose > 0) {
            cat(sprintf("delta theta: %.12f\n", sum(abs(theta0 - theta1))))
            if (verbose > 1) {
                cat(sprintf("delta coef:  %.12f\n", sum(abs(coef0 - .fixef(lobj)))))
                cat(sprintf("delta u:     %.12f\n", sum(abs(b.s0 - b.s(lobj)))))
                cat(sprintf("delta sigma: %.12f\n", abs(sigma0 - .sigma(lobj))))
                coef0 <- .fixef(lobj)
                b.s0 <- b.s(lobj)
                sigma0 <- .sigma(lobj)
                if (verbose > 2) {
                    cat("theta:  ", theta(lobj),"\n")
                    cat("coef:   ", .fixef(lobj),"\n")
                    cat("b.s:    ", b.s(lobj), "\n")
                    cat("sigmae: ", .sigma(lobj), "\n")
                }
            }
        }

        ## all zero or change smaller than relative tolerance
        ## all zero: we can't get out of this anyway, so we have to stop.
        converged <- all(theta1 == 0) || sum(abs(theta0 - theta1)) <
            200*rel.tol*sum(abs(theta0))
        if (verbose > 1)
            cat(sprintf("Criterion: %.12f, %.12f", sum(abs(theta0 - theta1)),
                        sum(abs(theta0 - theta1)) / rel.tol / sum(abs(theta0)) / 200), "\n")
        theta0 <- theta1
    }

    optinfo <- list(optimizer = "rlmer.fit.DAS",
                    conv = list(opt = 0),
                    feval = iter,
                    warnings = list(),
                    val = diff)

    if (iter == max.iter) {
        warning(wt <- "iterations did not converge, returning unconverged estimate.")
        optinfo$warnings <- list(wt)
        optinfo$conv$opt <- 1
    }
    lobj@optinfo <- optinfo
    return(lobj)
}