Struct Fit 

#[non_exhaustive]
pub struct Fit {

    pub beta: Vec<f64>,
    pub se: Vec<f64>,
    pub vcov: Vec<Vec<f64>>,
    pub tau2: Vec<f64>,
    pub dispersion: f64,
    pub diagnostics: Diagnostics,
    pub varcorr: Vec<Vec<f64>>,
    pub stddev_se: Vec<f64>,
    pub n_eval: usize,
    pub deviance: f64,
    pub loglik: f64,
    pub df: usize,
    pub reml: bool,
    pub fitted: Vec<f64>,
    pub ranef: Vec<f64>,
    pub ranef_levels: Vec<usize>,
}

Result of fit. Fixed-effect estimates cover all p predictors; SE and tau2 have the ranges below. Non-target SE slots are NaN.

#[non_exhaustive]: construct one only through fit_cold/fit_warm, and match on it with a trailing ... Everything the fit reports about ITSELF — as opposed to about the data — lives on Fit::diagnostics, which is #[non_exhaustive] for the same reason: a seventh diagnostic is then an additive change instead of a major version.

Fields (Non-exhaustive)

Non-exhaustive structs could have additional fields added in future. Therefore, non-exhaustive structs cannot be constructed in external crates using the traditional Struct { .. } syntax; cannot be matched against without a wildcard ..; and struct update syntax will not work.

beta: Vec<f64>

Fixed-effect estimates, length p.

se: Vec<f64>

Standard errors: se[j] = sqrt(Var(β̂_j)) for target predictors, NaN for non-targets. Length p.

vcov: Vec<Vec<f64>>

Fixed-effect covariance Cov(β̂) — a full symmetric p×p matrix (vcov[i][j]), NOT vech-packed: unlike varcorr, which packs one variable-sized block per grouping, this is a single square block and every consumer (R’s vcov()/confint, multcomp::glht, any hand-built Wald contrast) indexes it directly. se is its diagonal: se[j] == vcov[j][j].sqrt() wherever both are finite.

Finite exactly where se is — the off-diagonal vcov[i][j] is finite iff se[i] and se[j] both are. So a FitOptions::target_indices subset leaves everything outside the target block NaN (there is no covariance to report for a coefficient whose variance was never computed), and a non-converged fit is all-NaN.

Sources, by path: OLS/GLM/LMM invert the same Cholesky factor se’s forward solve already walks; GLMM WaldSe::Hessian takes the β block of the joint (θ,β) Hessian covariance (exact on every shape derivative::supports_shape accepts, finite-difference elsewhere), and WaldSe::Rx the p×p Schur inverse — both already formed in full and previously discarded down to a diagonal. Unlike stddev_se, this is populated on the Hessian’s RX fallback too (that fallback inverts a full p×p covariance; only a double failure, where the Hessian AND the fallback both fail, NaN-fills — as a non-converged fit).

tau2: Vec<f64>

Per-element Cholesky-scaled values theta[k]^2 * sigma_sq. These equal the random-effect variance components only for diagonal/scalar RE components (q=1 / scalar-extra — the currently reachable case); slope (q≥2) models are not yet validated through this field. Empty for OLS.

dispersion: f64

Estimated dispersion: φ for Gamma (Pearson moment estimator), the estimated shape θ for negative-binomial, the residual variance σ̂² for Gaussian — RSS/(n−p) for OLS (raw-row df, matching R summary.lm‘s sigma²; the same sigma_sq that scales se/vcov) and the REML pwrss/(n−p) for LMM (matching lme4 sigma()²; oracle: validation/goldens/sleepstudy_lmm.json sigma, asserted in fit_sleepstudy_slope_varcorr_matches_lme4) — and 1.0 for binomial/Poisson (where dispersion is fixed by the family, not estimated). NaN on any fit reporting converged == false — no route reports a dispersion off a fit that never reached an endpoint, including the fixed families’ 1.0 and a φ held fixed through FitOptions::dispersion.

diagnostics: Diagnostics

Everything the fit reports about itself: convergence, singularity, the aliased-column mask, the θ boundary state, which variance components were pinned there, and any Notes the solver raised. Single storage location — Fit::converged, Fit::singular and Fit::aliased are one-hop forwarders onto it, not copies.

varcorr: Vec<Vec<f64>>

RE (co)variance per grouping: one vech-packed lower-triangular covariance block D̂ = σ̂²·Λ̂Λ̂' per grouping, in declaration order (primary, then each extra). σ̂² is the residual scale for an LMM and the free GLMM scale pwrss/n for dispersion families (Gamma) — exactly the factor lme4’s VarCorr stddevs carry; it is ≡ 1 for binomial/Poisson/NB, and the same scale tau2 reports, so the two accessors agree. Vech order is column-major lower-triangular (matching the θ vech convention): for a q×q block, (0,0),(1,0),…,(q-1,0),(1,1),…,(q-1,q-1). Validated against lme4 VarCorr (validation/goldens/sleepstudy_lmm.json; Gamma scale: validation/goldens/sim_gamma_glmm.json). Empty for OLS/GLM (no random effects). This is the q≥2-valid replacement for tau2’s per-component variances; tau2 is retained for back-compat.

stddev_se: Vec<f64>

SE of each RE standard deviation, laid out like tau2 (per θ coordinate, length n_theta; primary block then each extra, in declaration order). Populated ONLY on a converged GLMM WaldSe::Hessian fit, from the θ block of the joint (θ,β) Hessian covariance that joint_hessian_cov inverts and otherwise discards (exact Hessian on every shape derivative::supports_shape accepts, finite-difference elsewhere). NaN under WaldSe::Rx, on the Hessian RX fallback, and for OLS/LMM (no Hessian machinery). Correct for SCALAR groupings only (the reachable GLMM case), where the RE stddev equals its θ so the θ-scale SE is the stddev SE directly; a q≥2 block would need a delta-method Jacobian.

Scale caveat for dispersion families (Gamma): this SE stays on the θ scale — lme4’s θ-Hessian convention — while varcorr/tau2 carry the σ̂² factor. Deliberately NOT multiplied by σ̂: the joint Hessian does not carry cov(σ̂, θ̂), so a σ̂-rescaled SE would be a delta-method value that matches no oracle. For the φ≡1 families the two scales coincide and there is no split. (The validation suite skips sd_se gating for dispersion families accordingly.)

n_eval: usize

Objective evaluations consumed by the θ (LMM) / joint [θ|β] (GLMM) BOBYQA search — GLMM counts both stages. 0 where no derivative-free optimizer runs (OLS/GLM closed-form or IRLS paths). Deterministic and clock-independent; the optimizer-grid campaign’s primary metric.

deviance: f64

Minimized optimizer criterion at the accepted point. LMM: the profiled REML deviance as computed by reml_deviance — equals lme4’s REMLcrit minus the data-independent constant df·(1 + ln 2π), df = n − p (validated against the frozen lme4 sleepstudy reference). GLMM: the marginal Laplace deviance d(y,ũ) + ‖ũ‖² + log|A| (see GlmmFit::deviance), which differs from −2·logLik by a data-only saturated constant. NaN for OLS/GLM and on optimizer/numerical failure (GLMM non-convergence surfaces as +∞ internally; mapped to NaN here). An LMM fit that hits MaxFunReached still reports the finite endpoint deviance here with converged == false — the plateau policy: a MaxFunReached cap-out reports its finite endpoint with converged == false rather than NaN-filling.

loglik: f64

The actual log-likelihood at the fitted parameters — deviance with its dropped data-only constants restored, on the logLik() scale (R/lme4):

• OLS/GLM — the standard closed forms (R logLik.lm/logLik.glm, MASS::glm.nb), including under prior weights.
• LMM — the REML criterion −REMLcrit/2 (this path is REML-only): −½(deviance + (n−p)·(1 + ln 2π)). REML criteria are comparable only between models with IDENTICAL fixed effects — an AIC/LRT across different fixed parts is meaningless; check Fit::reml before comparing. Matches lme4::logLik on a REML=TRUE fit.
• GLMM — the marginal Laplace/AGQ log-likelihood: −½·deviance + saturated_loglik (binomial/Poisson/NB, see family::saturated_loglik), −½·deviance for Gamma (lme4’s logLik(glmer) is −devfun/2 verbatim, gamma_aic’s +2 included — see fit::common::glmm_loglik). Matches lme4::logLik on the same fit, including the aggregated-binomial cbind(s, m−s) form under weights=.

NaN wherever deviance’s failure modes apply (non-converged/degenerate fits); finite on an LMM MaxFunReached endpoint, like deviance. AIC = 2·df − 2·loglik, BIC = df·ln(n) − 2·loglik with Fit::df.

df: usize

Parameters counted for AIC/BIC: retained fixed effects (p minus aliased columns, matching lme4’s NA-coefficient handling) + n_theta RE parameters + 1 if the family estimates a dispersion/scale (Gaussian σ², Gamma φ unless held fixed via FitOptions::dispersion, NB θ). 0 on degenerate NaN-fill paths.

reml: bool

true iff loglik is a REML criterion (the Gaussian LMM paths — REML is this engine’s locked LMM objective) rather than an ML log-likelihood. Model comparisons (AIC/LRT) across fits with different fixed effects are invalid when this is set — mirror of lme4’s REML-fit anova warning.

fitted: Vec<f64>

Fitted means μ̂ per row (length n): the conditional means through the inverse link — g⁻¹(Xβ̂ + Zb̂) for mixed fits (lme4 fitted()), g⁻¹(Xβ̂) for fixed-only. Empty on non-converged fits. Always on the original y scale, including on the Gaussian LMM paths — see fit::common::lmm_fitted for how the offset is restored there.

ranef: Vec<f64>

Random-effect conditional modes b̂ = Λ̂û on the natural (link) scale — lme4’s ranef() values. One block per grouping in declaration order (primary, then each extra — same order as varcorr), each block level-major: level l’s q values (intercept, then slopes in declaration order) at [l·q .. (l+1)·q]. η̂ = Xβ̂ + Zb̂ reproduces the linear predictor behind fitted. Level counts per grouping are in Fit::ranef_levels; q per grouping is recoverable from varcorr (vech length). For a nested grouping the block spans n_parents·n_per_parent slots (child id = parent·n_per_parent + within), padded with zero modes for parents with fewer observed children. Empty on non-converged fits.

This is the numbers only. To get them LABELLED — which level each row belongs to, with the padded slots dropped — use crate::formula::label_ranef, the one place the block layout is interpreted. Do not re-derive the slicing in a consumer: which layout a grouping lands in is a data-dependent speed decision, so a consumer that inferred it would be wrong on a dataset nobody tested.

The Gaussian LMM paths never form these during the fit — the profiled REML criterion does not need them — and recover them once at θ̂ by back-substitution afterwards, which is why they are here without the per-evaluation cost of forming them.

ranef_levels: Vec<usize>

Level count per grouping for slicing Fit::ranef, declaration order (primary, then each extra). ranef.len() = Σ_g levels[g]·q_g. Empty exactly when ranef is.

Implementations

impl Fit

pub fn converged(&self) -> bool

Whether the optimizer reached its convergence criterion — the most-read field in the API, forwarded so that moving it behind diagnostics costs its callers a () and nothing else. See Diagnostics::converged.

pub fn singular(&self) -> bool

Whether the fit is singular (lme4’s isSingular). Forwarder — see Diagnostics::singular.

pub fn aliased(&self) -> &[bool]

The length-p rank-deficiency mask. Forwarder — see Diagnostics::aliased.

pub fn stddev_corr(&self, group_idx: usize) -> (Vec<f64>, Vec<Vec<f64>>)

Reduce varcorr[group_idx] — a vech-packed (column-major lower-triangular, see Fit::varcorr) q×q covariance block — into per-dimension standard deviations and a full symmetric q×q correlation matrix, mirroring lme4’s VarCorr stddev/corr split. q is recovered from vech.len() = q(q+1)/2 via the quadratic formula. Column c’s entries start at vech offset c*q - c*(c-1)/2 (the same running cursor varcorr_block builds — for r >= c, idx(r,c) = c*q - c*(c-1)/2 + (r-c)); change together.