--- title: "Convergence, Boundaries, and Singularity" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Convergence, Boundaries, and Singularity} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE ) ``` ```{r load-package} library(mixeff) ``` Every fit returned by `lmm()` or `glmm()` carries a labeled convergence status, and the label is a claim with a precise scope. This vignette is the reference for that vocabulary: what `converged_interior` does and does not establish, how boundary and reduced-rank outcomes are recorded, how singularity is actually decided, how to check an optimum with `verify_convergence()`, and what happens when a requested GLMM estimator does not certify. The other vignettes cross-reference these semantics rather than restating them. ## What `converged_interior` establishes — and what it does not Start with a fit that goes well: a small longitudinal study with twelve clinics, a week effect, and a treatment contrast. ```{r clinic-data, include = FALSE} set.seed(42) n_clinic <- 12L n_per_clinic <- 6L clinic <- factor(rep(seq_len(n_clinic), each = n_per_clinic)) week <- rep(0:(n_per_clinic - 1L), n_clinic) treatment <- factor( rep(rep(c("usual", "coached"), each = n_per_clinic), length.out = length(week)), levels = c("usual", "coached") ) clinic_shift <- rnorm(n_clinic, sd = 0.45) score <- 7.2 - 0.28 * week - 0.55 * (treatment == "coached") + clinic_shift[as.integer(clinic)] + rnorm(length(week), sd = 0.35) clinic_visits <- data.frame(score, week, treatment, clinic) ``` ```{r interior-fit} fit <- lmm( score ~ week + treatment + (1 | clinic), clinic_visits, control = mm_control(verbose = -1) ) fit_status(fit) ``` `converged_interior` means the optimizer returned a finite interior solution that passed mixeff's declared numerical checks, with every variance component away from zero. It does not establish global optimality, correct model specification, or that the data carry enough information — no numerical status can — but it is the status every other block of a `summary()` assumes. The status changes when something needs your attention. The vocabulary has three values you will meet in practice: - `converged_interior` — a finite interior optimum, as above. - `converged_reduced_rank` — a genuine optimum at which the random-effect covariance has lower rank than the model requested; the next section walks through one. - `not_optimized` — the conservative label for a run that stopped in a flat or boundary region *without* certifying a stationary point. Estimates are present but the engine does not call them an optimum, and `summary()` repeats that state in a `Notes:` line rather than relying on the header alone. Optimization is silent and runs within one native call under bounded evaluation budgets; the `lmm()`/`glmm()` help pages and `?mm_control` document the budgets and how to cap them. What the fit ran is recorded on the object. `optimizer_certificate()` reports the optimizer, the stopping condition, and the certificate-time facts the status label summarizes: ```{r interior-certificate} optimizer_certificate(fit) ``` ```{r interior-checks, include = FALSE} stopifnot(inherits(fit, "mm_lmm")) stopifnot(identical(fit_status(fit), "converged_interior")) stopifnot(isFALSE(is_singular(fit))) stopifnot(nrow(optimizer_certificate(fit)$table) > 0L) ``` ## Boundary and reduced-rank fits are labeled facts Now a fit that does not go well — deliberately. Eighteen subjects, ten daily reaction-time measurements each, with subject intercepts and slopes that are nearly perfectly correlated by construction. This design produces a singular fit in `lme4` and in `mixeff` alike. ```{r singular-data, include = FALSE} set.seed(3) n_subj <- 18L days <- 0:9 b0 <- rnorm(n_subj, sd = 30) b1 <- 0.5 * b0 / 30 * 10 + rnorm(n_subj, sd = 0.5) sleep_like <- do.call(rbind, lapply(seq_len(n_subj), function(i) { data.frame( subj = factor(i), days = days, rt = 250 + b0[i] + (10 + b1[i]) * days + rnorm(length(days), sd = 20) ) })) ``` `lme4::lmer()` fits this model and reports the situation in one short line: `boundary (singular) fit`. The fact is correct. What it leaves implicit is *which* variance component reached the boundary, what the effective rank of the random-effect covariance is, and which downstream inference methods the boundary takes off the table. ```{r lme4-side, eval = requireNamespace("lme4", quietly = TRUE)} m <- suppressMessages(lme4::lmer( rt ~ days + (1 + days | subj), data = sleep_like )) m lme4::isSingular(m) ``` `mixeff` reports the same fact, and then unpacks it. `fit_status()` names the convergence outcome, `is_singular()` answers the yes/no question, `diagnostics()` returns stable codes with severity and stage, and `changes()` shows the requested-to-effective-to-fitted arc. ```{r singular-fit} sfit <- lmm( rt ~ days + (1 + days | subj), sleep_like, control = mm_control(verbose = -1) ) fit_status(sfit) is_singular(sfit) diagnostics(sfit)$table[, c("code", "severity", "stage", "message")] changes(sfit) ``` The requested random-effect covariance for `(1 + days | subj)` has rank 2: one dimension for the subject intercept and one for the subject slope. The fitted covariance has effective rank 1 — the model asked for two random-effect dimensions, but the data support only one effective dimension. That is recorded as the `covariance_reduced` diagnostic and the `converged_reduced_rank` status, and the printed object carries the tag in its header every time it prints — where `lme4` warns once at fit time and then falls silent. ```{r singular-print} print(sfit) ``` ```{r singular-checks, include = FALSE} stopifnot(identical(fit_status(sfit), "converged_reduced_rank")) stopifnot(isTRUE(is_singular(sfit))) diag_codes <- diagnostics(sfit)$table$code stopifnot("covariance_reduced" %in% diag_codes) ch <- changes(sfit)$table stopifnot("design_time_reduction" %in% ch$stage || "covariance_transition" %in% ch$stage) sfit_print <- capture.output(print(sfit)) stopifnot(any(grepl("reduced_rank|boundary", sfit_print))) ``` A reduced-rank covariance does not automatically invalidate the model. It is a labeled fact about the fit, recorded in the object rather than printed once and lost — and it changes which inference routes the package will certify (one sentence on that at the end of this vignette). ## How singularity is decided: rank, not θ entries To read the boundary structure itself, look at the model's internal coordinates. The random-effect covariance Σ is not optimized directly. The engine searches over a vector θ that fills a lower-triangular factor Λ, with Σ = σ²ΛΛᵀ; θ is by construction exactly the set of free entries of Λ. This is the parameterization Bates et al. (2015) describe for `lme4`; `MixedModels.jl` and the engine behind `mixeff` inherit it. `parameterization()` shows the mapping for the singular fit: ```{r singular-parameterization} parameterization(sfit)$table[, c("theta_name", "theta_value", "lambda_value")] ``` The *diagonal* entry for the `days` slope is pinned at numerical zero — that is the boundary event behind the `converged_reduced_rank` status. It is tempting to generalize this to "a zero θ entry means the covariance is reduced-rank". That reading is wrong. A zero **off-diagonal** entry of Λ only says the corresponding covariance is zero; the covariance matrix can keep full rank. Singularity is a property of the covariance matrix's rank — equivalently its eigenstructure — and in this parameterization the rank collapses exactly when a **diagonal** entry of Λ is zero. A two-by-two example makes the difference concrete: ```{r lambda-rank-demo} l11 <- 1.2 l21 <- 0.4 l22 <- 0.9 # Zero the OFF-DIAGONAL entry (l21 = 0): zero covariance, full rank. L_offdiag_zero <- matrix(c(l11, 0, 0, l22), nrow = 2) Sigma_offdiag_zero <- L_offdiag_zero %*% t(L_offdiag_zero) Sigma_offdiag_zero eigen(Sigma_offdiag_zero, symmetric = TRUE)$values # Zero a DIAGONAL entry (l22 = 0): the rank collapses to 1. L_diag_zero <- matrix(c(l11, l21, 0, 0), nrow = 2) Sigma_diag_zero <- L_diag_zero %*% t(L_diag_zero) Sigma_diag_zero eigen(Sigma_diag_zero, symmetric = TRUE)$values ``` With the off-diagonal at zero, both eigenvalues are positive: the intercept and slope are uncorrelated, and each still has its own variance — an ordinary full-rank fit (the model `(1 + days || subj)` *asks* for exactly this shape). With a diagonal entry at zero, one eigenvalue is zero and the two random effects collapse onto a single dimension: that is singularity. ```{r lambda-rank-checks, include = FALSE} stopifnot(identical(qr(Sigma_offdiag_zero)$rank, 2L)) stopifnot(all(eigen(Sigma_offdiag_zero, symmetric = TRUE)$values > 0)) stopifnot(identical(qr(Sigma_diag_zero)$rank, 1L)) ev_diag <- eigen(Sigma_diag_zero, symmetric = TRUE)$values stopifnot(abs(ev_diag[[2L]]) < 1e-12) stopifnot(identical(Sigma_offdiag_zero[1L, 2L], 0)) # The fitted singular case: the slope's DIAGONAL theta entry is at zero. pt <- parameterization(sfit)$table stopifnot(pt$theta_value[pt$theta_name == "theta[0:days,days]"] < 1e-6) ``` Consistent with this, `is_singular()` does not scan θ for zeros. It reads the covariance status the engine recorded on the fit — the fit status plus the per-block effective-covariance status — so its answer is the engine's rank verdict, not an R-side reconstruction. ## Verifying an optimum: `verify_convergence()` A convergence status describes the run that produced the fit. If you want independent evidence, `verify_convergence()` re-runs the fit under the engine's bounded verification workflow: a restart from the fitted optimum, one or more jittered restarts, and (LMM-only, opt-in) an alternate-optimizer consensus pass. The engine judges agreement against explicit tolerances on the objective, θ, and the fixed effects, and owns the verdict; R only formats it. ```{r verify-lmm} vc <- verify_convergence(fit) vc ``` Each run either agrees with the reference optimum within the printed tolerances or it does not; the verdict summarizes the runs (`restart_agrees`, `optimizer_consensus`, `fragile`, `unstable`, or `not_run`). A verdict of `restart_agrees` is evidence of a *stable* optimum — the bounded re-runs came back to the same place — not proof of a global one. The `consensus` pass defaults to `FALSE` because this vendored build compiles without the optional `nlopt` backend; for some models the engine's alternate optimizer choice is an nlopt optimizer, and its absence would read as a spurious `fragile` verdict that reflects the build, not the fit. Verification and rank are different questions. Running the verifier on the singular fit from above also comes back `restart_agrees`: a reduced-rank optimum can be a perfectly stable optimum. ```{r verify-singular} vs <- verify_convergence(sfit) vs$status ``` The GLMM method uses the same workflow with wider default tolerances — objective `1e-4` and beta `1e-3`, against the LMM defaults of `1e-5` and `1e-4` — because every GLMM objective evaluation carries inner-PIRLS noise and the joint path's derivative-free beta search bounds beta reproducibility at its stopping tolerance. The consensus pass is refused for GLMM fits (the GLMM fit drivers own optimizer selection). ```{r glmm-data, include = FALSE} set.seed(7) g <- factor(rep(1:12, each = 12)) x <- rnorm(144) eta <- -0.3 + 0.8 * x + rep(rnorm(12, sd = 0.5), each = 12) trial <- data.frame(y = rbinom(144, 1, plogis(eta)), x = x, g = g) ``` ```{r verify-glmm} gfit <- glmm( y ~ x + (1 | g), trial, family = binomial(), control = mm_control(verbose = -1) ) gv <- verify_convergence(gfit) gv ``` ```{r verify-checks, include = FALSE} stopifnot(identical(vc$status, "restart_agrees")) stopifnot(all(vc$table$agrees)) stopifnot(identical(vc$tolerances$objective, 1e-5)) stopifnot(identical(vc$tolerances$beta, 1e-4)) stopifnot(identical(vs$status, "restart_agrees")) stopifnot(identical(gv$status, "restart_agrees")) stopifnot(all(gv$table$agrees)) stopifnot(identical(gv$tolerances$objective, 1e-4)) stopifnot(identical(gv$tolerances$beta, 1e-3)) ``` ## When the requested estimator does not certify (GLMM) `glmm(method = "joint_laplace")` requests the joint-Laplace estimator (the route that certifies Wald inference; `vignette("glmm")` is the operational quickstart). On some designs the joint route does not certify an optimum. The engine then returns a *labelled* fallback fit from the fast profiled-PIRLS estimator instead of failing or — worse — silently presenting the fallback as the thing you asked for. The substitution crosses the boundary as a typed record, in keeping with the package's no-silent-surgery contract: - **The typed record.** `optimizer_certificate(fit)$raw$estimator_substitution` carries the requested method, the effective (fallback) method, the requested route's fit status and return code, and the engine's reason. It is the machine-readable contract surface; scripts should gate on it, not on message text. - **The fit-time notice.** A condition of class `mm_estimator_substitution_notice` announces the substitution when the fit returns (suppressible with `mm_control(verbose = -1)`). - **The ungated summary note.** `summary()` prints an `estimator substitution:` note naming both estimators regardless of verbosity settings, and independent of the fallback's own `converged_*` status — a clean-looking status line must not hide the fact that a different estimator than requested produced the numbers. - **The checking verb.** `verify_convergence()` on a substituted fit re-requests the joint route; when the verification refit falls back the same way — the expected case — every run verifies the profiled objective the fitted numbers actually came from. On an ordinary fit the record is absent, and its absence is meaningful (no substitution happened): ```{r substitution-surface} optimizer_certificate(gfit)$raw$estimator_substitution ``` ```{r substitution-checks, include = FALSE} stopifnot(is.null(optimizer_certificate(gfit)$raw$estimator_substitution)) ``` A real example, from reproducing a published `glmer()` analysis of OSF study data (an error-rate model with an open-practice × year interaction): requesting `method = "joint_laplace"` with the raw, uncentered `Year` predictor hands the engine an ill-scaled design, the joint route does not certify, and the result comes back as a labelled fast-PIRLS fallback — fit-time notice, summary note, and typed record all present, with Wald inference withheld. Centering the year (`Year - 2015`) makes the joint route certify and match `glmer()` on the same data. mixeff flags this situation at fit time: a scaling advisory names the ill-scaled columns, and centering or rescaling the predictor — as the year centering does here — is the standard remedy. ## What a boundary refuses A boundary fit also changes which inference routes the package will certify: Satterthwaite and Kenward–Roger degrees of freedom, for example, are refused on singular fits by mixeff's validated-support policy, with stable reason codes and a bootstrap alternative offered under the same labeling rules — the route map, the refusal codes, and the alternatives are documented in `vignette("inference")`. ## References - Bates, D., Mächler, M., Bolker, B., and Walker, S. (2015). Fitting linear mixed-effects models using lme4. *Journal of Statistical Software*, 67(1), 1–48.