MixedEffectRegressor#

class mmer.MixedEffectRegressor(fixed_effects_model, max_iter=30, tol=1e-06, patience=3, correction_method='bste', slq_steps=30, n_probes=60, preconditioner=True, cg_maxiter=1000, n_jobs=-1, backend='threading')[source]

Bases: object

Multivariate Mixed Effects Regression (MMER) using Expectation-Maximization.

Fits mixed model with multiple responses, supporting arbitrary grouping factors and linear random slopes. Solves for random effects and residual covariances using EM algorithm with stochastic log-determinant estimation.

Parameters:

  • fixed_effects_model (RegressorMixin) – Base regressor for fixed effects (must support multi-output).

  • max_iter (int, default=20) – Maximum number of EM iterations.

  • tol (float, default=1e-6) – Convergence tolerance on log-likelihood relative change.

  • patience (int, default=3) – Number of iterations to wait for likelihood improvement before early stopping. Setting to a high value effectively disables early stopping and relies solely on tol.

  • correction_method (str, default='bste') –

    Method for variance correction in M-step:

    • ’bste’: block stochastic trace estimation

    • ’de’: deterministic estimation

  • slq_steps (int, default=30) – Number of Lanczos steps for Stochastic Lanczos Quadrature (log-det estimation). A range of 30-50 is typically sufficient. Higher values yield slightly more accurate estimates but increase computation time and risk numerical instability.

  • n_probes (int, default=60) – Number of random probes used for both SLQ log-determinant estimation and stochastic variance correction. A fixed target around 50-60 is usually optimal independent of matrix dimension for O(1/sqrt(p)) error convergence.

  • preconditioner (bool, default=True) – Whether to use residual-based preconditioner for CG solver.

  • n_jobs (int, default=-1) – Number of parallel jobs for SLQ and trace estimation (-1 uses all cores). Setting to number of outputs (m) is recommended for optimal performance.

  • backend (str, default='threading') – Joblib parallel backend (‘threading’, ‘loky’). Setting to ‘threading’ is highly recommended as almsot always the solver is Woodbury-based. Setting to ‘loky’ can be used rarely for cg-based solvers under very specific conditions of large m and small n where woodbury is not beneficial.

  • cg_maxiter (int)

Variables:

  • fe_model (RegressorMixin) – Fitted fixed effects model.

  • random_effect_terms (tuple of RandomEffectTerm) – Fitted random effect terms containing covariance matrices.

  • residual_term (ResidualTerm) – Fitted residual term containing residual covariance matrix.

  • log_likelihood (list of float) – Log-likelihood values across EM iterations.

  • n (int) – Number of observations.

  • m (int) – Number of output dimensions.

  • k (int) – Number of grouping factors.

Examples

>>> from sklearn.linear_model import Ridge
>>> model = MixedEffectRegressor(fixed_effects_model=Ridge())
>>> results = model.fit(X, y, groups, random_slopes=([0, 1], None))
>>> predictions = model.predict(X_new)
prepare_data(X, y, groups, validation_split=0.0, validation_group=0)[source]

Prepare data for EM algorithm by creating realized objects.

Generates transient realized random effects and residual for the current dataset. Optionally splits data into training and validation sets based on group membership.

Parameters:

  • X (np.ndarray) – Covariates, shape (n, p).

  • y (np.ndarray) – Multi-output targets, shape (n, m).

  • groups (np.ndarray) – Grouping factors, shape (n, k).

  • validation_split (float, default=0.0) – Fraction of groups to use for fixed-effects validation (0.0 means no validation). Setting to a non-zero value means fixed effects can accept validation data.

  • validation_group (int, default=0) – Column index in groups to use for group-wise validation splitting.

Returns:

  • marginal_residual (np.ndarray) – Initial marginal residual, raveled shape (m*n,).

  • realized_effects (tuple of RealizedRandomEffect) – Realized random effect objects.

  • realized_residual (RealizedResidual) – Realized residual term.

fit(X, y, groups, random_slopes=None, validation_split=0.0, validation_group=0)[source]

Fit the MMER model using the EM algorithm.

Parameters:

  • X (np.ndarray) – Covariates, shape (n, p) where n is number of observations and p is number of features.

  • y (np.ndarray) – Multi-output targets, shape (n, m) where m is number of outputs.

  • groups (np.ndarray) – Grouping factors, shape (n, k) where k is number of grouping factors. Each column represents a different grouping structure.

  • random_slopes (tuple of list of int, optional) – Tuple of lists specifying random slopes for each grouping factor. Each list contains column indices in X for random slopes corresponding to that group. None or empty list implies random intercept only for that group. If None, all groups get random intercepts only.

  • validation_split (float, default=0.0) – Fraction of groups to use for fixed-effects validation (early stopping). Must be between 0.0 and 1.0. Set to 0.0 to disable validation. Setting to a non-zero value means fixed effects can accept validation data.

  • validation_group (int, default=0) – Column index in groups to use for group-wise validation splitting.

Returns:

Fitted result object containing covariance estimates and diagnostics.

Return type:

MixedEffectResults

Examples

>>> # Fit model with random intercepts only
>>> results = model.fit(X, y, groups)

>>> # Fit with random slopes on features 0 and 1 for first group
>>> results = model.fit(X, y, groups, random_slopes=([0, 1], None))

>>> # Fit with validation split
>>> results = model.fit(X, y, groups, validation_split=0.2)