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) LMMSolution — result wrapper for LMM GLMMSolution — result wrapper for GLMM

pystatistics.mixed.lmm(y, X, groups, *, random_effects=None, random_data=None, reml=True, tol=1e-08, max_iter=200, compute_satterthwaite=True)[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 (ArrayLike) – Response vector (n,).

  • X (ArrayLike) – Fixed effects design matrix (n, p). Should include an intercept column if desired.

  • groups (dict[str, ArrayLike]) – 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, ArrayLike] | 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.

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)[source]

Fit a generalized linear mixed model.

Uses Laplace approximation to the marginal likelihood with Penalized IRLS (PIRLS) for the inner loop and L-BFGS-B for the outer optimization over variance components.

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.

Returns:

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

Return type:

GLMMSolution

class pystatistics.mixed.LMMSolution(_result)[source]

Bases: object

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])

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 se: 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 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
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)[source]

Bases: object

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])

property params: GLMMParams
property coefficients: ndarray[tuple[Any, ...], dtype[_ScalarT]]
property fixef: dict[str, float]
property se: 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 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
summary()[source]

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

Return type:

str