Mixed Models¶

Linear mixed models (LMM) and generalized linear mixed models (GLMM). Random intercepts and slopes, nested and crossed random effects, REML/ML estimation, Satterthwaite degrees of freedom, BLUPs, ICC, likelihood ratio tests.

Mixed models: Linear Mixed Models (LMM) and Generalized Linear Mixed Models (GLMM).

Public API:: lmm() — fit a linear mixed model (REML or ML) glmm() — fit a generalized linear mixed model (Laplace approximation) grm_lmm() — fit a low-rank / GRM mixed model (CPU/GPU; genomics regime) LMMSolution — result wrapper for LMM GLMMSolution — result wrapper for GLMM GRMSolution — result wrapper for the low-rank / GRM mixed model

pystatistics.mixed.lmm(y, X, groups, *, random_effects=None, random_data=None, reml=True, tol=1e-08, max_iter=200, compute_satterthwaite=True, conf_level=0.95)[source]¶

Fit a linear mixed model.

Estimates fixed effects β, random effects variance components, and conditional modes (BLUPs) of random effects using the profiled REML/ML deviance approach from Bates et al. (2015).

Parameters:

y (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]) – Response vector (n,).
X (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]) – Fixed effects design matrix (n, p). Should include an intercept column if desired.
groups (dict[str, _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]]) – Dict mapping grouping factor names to group label arrays. Example: {‘subject’: subject_ids}.
random_effects (dict[str, list[str]] | None) – Optional dict mapping group names to lists of random effect terms. Default: random intercept per group. Example: {‘subject’: [‘1’, ‘time’]} for (1 + time | subject).
random_data (dict[str, _Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]] | None) – Optional dict mapping variable names to data arrays for random slope variables. Example: {‘time’: time_array}.
reml (bool) – If True (default), use REML estimation. If False, use ML. Use ML (reml=False) for likelihood ratio tests between models with different fixed effects.
tol (float) – Convergence tolerance for the optimizer. Default 1e-8.
max_iter (int) – Maximum optimizer iterations. Default 200.
compute_satterthwaite (bool) – If True (default), compute Satterthwaite denominator df for fixed effects. Set to False for speed if p-values are not needed.
conf_level (float)

Returns:

LMMSolution with fixed effects, random effects, variance components, model fit statistics, and R-style summary().

Return type:

LMMSolution

Examples

# Random intercept model >>> result = lmm(y, X, groups={‘subject’: subject_ids})

# Random intercept + slope >>> result = lmm(y, X, groups={‘subject’: subject_ids}, … random_effects={‘subject’: [‘1’, ‘time’]}, … random_data={‘time’: time_array})

# Crossed random effects >>> result = lmm(y, X, groups={‘subject’: subj, ‘item’: item})

pystatistics.mixed.glmm(y, X, groups, *, family='binomial', random_effects=None, random_data=None, tol=1e-08, max_iter=200, conf_level=0.95)[source]¶

Fit a generalized linear mixed model.

Uses the Laplace approximation to the marginal likelihood (equivalent to lme4::glmer(..., nAGQ=1), R’s default): the random-effects modes are profiled by Penalized IRLS (PIRLS) in an inner loop, while BOTH the variance parameters θ and the fixed effects β are optimized in the outer loop (L-BFGS-B) over the Laplace deviance. Optimizing β in the outer loop — rather than solving it jointly inside PIRLS (the cruder nAGQ=0 scheme) — is what makes the fixed effects match glmer’s Laplace fit rather than a biased approximation.

Parameters:

y (ArrayLike) – Response vector (n,).
X (ArrayLike) – Fixed effects design matrix (n, p).
groups (dict[str, ArrayLike]) – Dict mapping grouping factor names to group label arrays.
family (str | Family) – GLM family specification. String (‘binomial’, ‘poisson’) or a Family instance from pystatistics.regression.families.
random_effects (dict[str, list[str]] | None) – Optional random effects specification.
random_data (dict[str, ArrayLike] | None) – Optional data for random slope variables.
tol (float) – Convergence tolerance.
max_iter (int) – Maximum optimizer iterations.
conf_level (float)

Returns:

GLMMSolution with fixed effects, random effects, and model fit.

Return type:

GLMMSolution

pystatistics.mixed.grm_lmm(y, X, W, *, backend=None, reml=True, names=None, tol=1e-08, max_iter=200, conf_level=0.95, force=False)[source]¶

Fit a low-rank / GRM mixed model: y = Xβ + g + ε, g ~ N(0, σ²_g·WW’/M).

Parameters:

y (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]) – Response vector (n,).
X (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]) – Fixed-effects design (n, p) — include an intercept column if wanted.
W (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]) – Low-rank factor (n, M) defining K = WW’/M (e.g. a standardized genotype matrix). M is the rank.
backend (str | None) – Compute backend — 'cpu' (float64 reference), 'gpu' (float32 speed path), 'gpu_fp64' (CUDA-only exact), or 'auto'. None resolves to 'cpu' for a numpy input.
reml (bool) – REML (default) or ML.
names (tuple[str, ...] | None) – Optional fixed-effect names (length p).
tol (float) – Optimizer tolerance for the θ search.
max_iter (int) – Maximum optimizer evaluations.
conf_level (float) – Confidence level for conf_int.
force (bool) – Bypass the float32 Gram-conditioning gate on the GPU path. Use only when you know W is well-conditioned; results on an ill-conditioned W will be unreliable.

Returns:

GRMSolution with fixed effects, variance components, heritability, genetic-value BLUPs, and fit statistics.

Raises:

ValidationError – bad shapes / non-finite input / unknown backend.
RuntimeError – a GPU backend requested but unavailable (or gpu_fp64 on MPS).
NumericalError – the float32 GPU path refused an ill-conditioned W (raise the precision with backend='gpu_fp64' or use backend='cpu').

Return type:

GRMSolution

class pystatistics.mixed.LMMSolution(_result, _conf_level=0.95)[source]¶

Bases: SolutionReprMixin

Solution wrapper for a fitted linear mixed model.

Provides R-style summary output matching lmerTest::summary(), property accessors for fixed effects, random effects, ICC, and model comparison via likelihood ratio test.

Parameters:

_result (Result[LMMParams])
_conf_level (float)

property params: LMMParams¶

property coefficients: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶: Fixed effect estimates β̂.

property fixef: dict[str, float]¶: Fixed effects as name → value dict.

property standard_errors: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶: Standard errors of fixed effects.

property t_values: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶: t-statistics for fixed effects.

property p_values: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶: p-values for fixed effects (Satterthwaite df).

property df_satterthwaite: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶: Satterthwaite denominator df for each fixed effect.

property conf_level: float¶: Confidence level for conf_int (default 0.95).

property conf_int: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

Wald confidence intervals for the fixed effects, shape (p, 2).

β ± t * SE using the Student-t quantile at each coefficient’s Satterthwaite denominator df (the finite-sample reference LMM uses for its p-values, matching lmerTest).

property ranef: dict[str, ndarray[tuple[Any, ...], dtype[_ScalarT]]]¶: Random effects (BLUPs / conditional modes) per grouping factor.

property var_components: tuple[VarCompSummary, ...]¶: Variance component summaries.

property icc: dict[str, float]¶

Intraclass correlation coefficient per grouping factor.

ICC = σ²_group / (σ²_group + σ²_residual)

For models with random slopes, uses the intercept variance only.

property log_likelihood: float¶

property aic: float¶

property bic: float¶

property fitted_values: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property residuals: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property converged: bool¶

property n_iter: int¶

property is_singular: bool¶

Whether this is a boundary (singular) fit.

True when a random-effects variance has collapsed to (near) zero or an implied correlation has reached ±1 — the fit sits on the boundary of the feasible region. Mirrors lme4’s isSingular(). The estimates are still the correct (boundary) MLE; a singular fit is a signal the random-effects structure may be too complex for the data, not an error.

compare(other)[source]¶

Likelihood ratio test between two nested models.

Both models should be fit with ML (reml=False) for valid LRT.

Parameters:: other (LMMSolution) – The other model to compare against.
Returns:: Formatted LRT summary string.
Return type:: str

summary()[source]¶

R-style summary matching lmerTest::summary(lmer(…)).

Return type:: str

class pystatistics.mixed.GLMMSolution(_result, _conf_level=0.95)[source]¶

Bases: SolutionReprMixin

Solution wrapper for a fitted generalized linear mixed model.

Same interface as LMMSolution plus family-specific properties. Uses Wald z-statistics (not Satterthwaite t) for inference.

Parameters:

_result (Result[GLMMParams])
_conf_level (float)

property params: GLMMParams¶

property coefficients: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property fixef: dict[str, float]¶

property standard_errors: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property z_values: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶: Wald z-statistics for fixed effects.

property p_values: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property conf_level: float¶: Confidence level for conf_int (default 0.95).

property conf_int: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

Wald confidence intervals for the fixed effects, shape (p, 2).

β ± z * SE with the normal quantile for conf_level (GLMM inference is asymptotic-normal).

property ranef: dict[str, ndarray[tuple[Any, ...], dtype[_ScalarT]]]¶

property var_components: tuple[VarCompSummary, ...]¶

property icc: dict[str, float]¶

ICC on the latent (link) scale.

For GLMM, ICC is computed on the link scale: ICC = σ²_group / (σ²_group + π²/3) for logistic ICC = σ²_group / (σ²_group + 1) for probit

property log_likelihood: float¶

property deviance: float¶

property aic: float¶

property bic: float¶

property fitted_values: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶: Fitted values on the response scale (μ̂).

property linear_predictor: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶: Linear predictor (η̂ = Xβ̂ + Zb̂).

property residuals: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property converged: bool¶

property n_iter: int¶

summary()[source]¶

R-style summary matching lme4::summary(glmer(…)).

Return type:: str

class pystatistics.mixed.GRMSolution(_result, _conf_level=0.95)[source]¶

Bases: SolutionReprMixin

Solution wrapper for a fitted low-rank / GRM mixed model.

Exposes the uniform accessors (coefficients, standard_errors, z_values, p_values, conf_int, fitted_values, residuals, converged, n_iter, backend_name) plus the quantitative-genetics quantities (heritability, var_genetic, var_residual, variance_ratio, genetic_values).

Parameters:

_result (Result[GRMParams])
_conf_level (float)

property params: GRMParams¶

property coefficients: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property coef: dict[str, float]¶

property standard_errors: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property z_values: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶: Wald z-statistics for fixed effects (asymptotic-normal reference).

property p_values: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property conf_level: float¶

property conf_int: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

Wald confidence intervals for the fixed effects, shape (p, 2).

β ± z * SE with the normal quantile for conf_level (GRM REML fixed-effect inference is asymptotic-normal).

property heritability: float¶: Narrow-sense heritability h² = σ²_g / (σ²_g + σ²_e).

property var_genetic: float¶: Genetic variance component σ²_g.

property var_residual: float¶: Residual variance component σ²_e.

property variance_ratio: float¶: Variance ratio γ = σ²_g / σ²_e.

property genetic_values: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶: BLUP of the genetic values g (n,).

property log_likelihood: float¶

property aic: float¶

property bic: float¶

property fitted_values: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property residuals: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property converged: bool¶

property n_iter: int¶

property backend_name: str¶

summary()[source]¶

R-style summary for the low-rank / GRM mixed model.

Return type:: str