glmm

Standalone f64 GLM(M) fit kernels — OLS → GLM → LMM → GLMM — in pure Rust on faer.

Fits fixed-effect and mixed (random-intercept/random-slope) models for Gaussian, Binomial (logit/probit/cloglog), Poisson, Gamma, Negative-Binomial, and (fixed-effect GLM only) Inverse-Gaussian outcomes, validated against R/lme4 and Julia/MixedModels.jl goldens.

Beta. Handles the usual mixed models, supports AGQ more widely than lme4, and is typically several times faster on them. Some lme4 features are still missing (see docs). Version 1.0 will add a built-in test for whether a random effect matters — currently you need the separate RLRsim package, and it only works for normal (Gaussian) outcomes.

New to the crate? Start with tutorial-rust.md — a single-page, three-layer walkthrough (cold fit → warm fit → advanced loop) plus a short section on parsing an R-style formula string instead of building inputs by hand. Python and R packages are also available — see Python and R below.

One crate, one fit

A single entry point covers the whole linear-regression family. fit reads ModelSpec and routes on family × random effects × weights:

flowchart LR
    F["fit(x, y, ModelSpec)"] --> RE{"re?"}
    RE -- "None" --> FAM1{"family?"}
    FAM1 -- "Gaussian" --> OLS["OLS - WLS with weights"]
    FAM1 -- "other" --> GLM["GLM - IRLS"]
    RE -- "Some" --> FAM2{"family?"}
    FAM2 -- "Gaussian" --> LMM["LMM - profiled REML"]
    FAM2 -- "other" --> GLMM["GLMM - PIRLS, Laplace/AGQ"]

This is the pitch view. Every branch point, solver path, and tuning knob is traced to code in the full algorithm map: documentation/algorithms.md, with the LMM and GLMM legs detailed in documentation/algorithms-lmm.md and documentation/algorithms-glmm.md.

Documentation

The fuller map — entry point, tutorials, examples, migration guides, reference, and internals, all in one tiered table — lives in documentation/index.md. The table below is this README's own copy of the essentials, not the whole set.

Scope and stability (0.3.x)

The semver-covered surface is fit_cold/fit_warm + ModelSpec + GroupIds.

Every wired family, except Inverse-Gaussian (fixed-effects only), fits through both routes — there is no reachable panic for falling outside the dense solver's envelope (too many extra groupings, or an extra grouping too wide); classification just routes to the sparse-Z solver instead. See tutorial-rust.md and documentation/algorithms-glmm.md for the dense/sparse routing envelope.

The loop_advanced cargo feature (off by default) exposes an unstable scratch-explicit hot-path surface for warm-start callers like MCPower's simulation loop — no semver guarantees; do not depend on it outside a pinned revision.

The orchestrate cargo feature (off by default) exposes the string-typed fit orchestration the Python and R packages are built on — one definition of the family/link string vocabulary, the formula-and-data lowering, and the flattened result both ports publish. No semver guarantees, for the same reason as loop_advanced: its shape follows the ports' needs.

The parallel cargo feature (off by default, experimental) enables in-fit parallelism — the AGQ cluster loop and the FD-Hessian SE grid — via rayon, and is additionally gated at runtime by FitOptions::parallel_inner (also off by default): both the feature and the flag must be on for any thread to spawn. Parallel results are bit-identical to serial ones, but the kernels are new and their performance envelope isn't characterized yet — treat as opt-in only. A no-op on wasm32 (compile-time excluded, rayon never pulled in).

Quick example

use glmm::{fit_cold, Family, FitOptions, GroupIds, ModelSpec, ReStructure, Sizing};

// y ~ x + (1 | group), Gaussian — a random-intercept LMM.
let model = ModelSpec {
    family: Family::Gaussian,
    re: Some(ReStructure {
        sizing: Sizing::FixedClusters { n_clusters: 6 },
        slopes: vec![],
        extra_groupings: vec![],
    }),
};
let ids = GroupIds { primary: vec![0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5], extra: vec![] };
// target_indices = which coefficients get a standard error; here, both.
let opts = FitOptions { target_indices: vec![0, 1], ..Default::default() };

let fit = fit_cold(&x, &y, n, p, &model, &ids, &opts);
assert!(fit.converged());

See tutorial-rust.md for the full walkthrough: warm starts, the advanced hot-loop surface, and building x/ModelSpec/GroupIds from a formula string with glmm::formula (the formula feature, on by default) instead of by hand.

Python and R

The same kernel ships as a Python package and an R package. Both take a data table and an R-style formula — no design matrices by hand.

Python (glmm on PyPI; Python 3.10+, NumPy is the only dependency):

pip install glmm

import glmm

fit = glmm.fit(data, "y ~ x1 + (1 | group)")   # data: dict / pandas / polars
fit.summary()

The public surface is eight names: glmm.fit, glmm.Fit, and the six warning categories the diagnostics channel raises. See the Python README and tutorial-python.md.

R (fastglmm, via r-universe):

install.packages("fastglmm", repos = c("https://pawlenartowicz.r-universe.dev", getOption("repos")))

library(fastglmm)

fit <- fastglmm(y ~ x1 + (1 | group), data)
summary(fit)   # plus fixef, ranef, fitted, vcov, VarCorr, confint, logLik, isSingular

Deliberately scoped to fast fitting — anything the engine cannot compute honestly (predict, residuals) errors with the reason instead of guessing. See the R README and tutorial-r.md.

Design

Origin

glmm was carved out of MCPower's simulation engine — the numerics here are the same validation-pinned kernels that power MCPower's Monte Carlo fits, split out so they're usable standalone. The loop_advanced feature above exists specifically to serve MCPower's warm-start hot loop as a consumer of this crate.

License

GPL-3.0-or-later (coupled to the GPL-3 MCPower flagship).

Paweł Lenartowicz — Freestyler Scientist · GitHub · ORCID