fit¶

The primary entry point for fitting linear mixed models with crossed random effects. Accepts both random intercepts (via groups) and random slopes (via random). Works with any narwhals-compatible DataFrame (pandas, polars, etc.).

fit(formula, data, groups=None, method='REML', random=None, optimizer='lbfgsb', theta0=None, weights=None, offset=None, correlation=None, df_method='satterthwaite')[source]¶

Fit a linear mixed model with crossed random effects via profiled REML.

Parameters:

formula – Patsy-syntax fixed-effects formula, e.g. "y ~ x1 + x2".
data – DataFrame containing all variables. Any narwhals-compatible frame (pandas, polars, …) is accepted.
groups – Column name (str) or list of column names for crossed random intercepts. Shorthand for random=["(1|g1)", "(1|g2)", ...]. The first entry is the primary grouping factor.
random – List of lme4-style random effect specifications, e.g. ["(1 + x | g1)", "(1 | g2)"]. Supports correlated (|) and independent (||) parameterisations. Takes precedence over groups when both are provided.
method – Estimation method. "REML" (default) maximises the restricted likelihood and is recommended for variance-component estimation. "ML" maximises the full likelihood and is required for likelihood-ratio tests (LRT) comparing models with different fixed effects via interlace.anova().
optimizer – Optimizer for the REML criterion. "lbfgsb" (default) uses scipy.optimize.minimize with L-BFGS-B. "bobyqa" uses pybobyqa (requires the bobyqa optional extra), a gradient-free trust-region method that is more robust near variance-parameter boundaries and matches lme4’s default.
theta0 – Initial theta for the optimizer. Defaults to np.ones(n_theta). Pass the theta attribute of a previously fitted model to warm-start the optimizer (e.g. for case-deletion refits).
weights – Observation-level prior weights, shape (n,). Each observation’s contribution to the log-likelihood is scaled by its weight. Equivalent to pre-multiplying all data by sqrt(diag(weights)). Defaults to ones (unweighted).
offset – Offset vector, shape (n,). A known term added to the linear predictor that is not estimated. For LMMs (identity link) this is equivalent to fitting y - offset ~ …. Defaults to zero.

Return type:

CrossedLMEResult

Returns:

CrossedLMEResult – Drop-in replacement for statsmodels MixedLMResults.

Parameters:

formula (str)
data (Any)
groups (str | list[str] | None)
method (str)
random (list[str] | None)
optimizer (str)
theta0 (ndarray | None)
weights (ndarray | None)
offset (ndarray | None)
correlation (Any | None)
df_method (str)

Examples

>>> import interlace
>>> import pandas as pd
>>> df = pd.DataFrame({
...     "y": [1.0, 2.0, 3.0, 1.5, 2.5, 3.5],
...     "x": [0.1, 0.2, 0.3, 0.1, 0.2, 0.3],
...     "g": ["a", "a", "a", "b", "b", "b"],
... })
>>> result = interlace.fit("y ~ x", df, groups="g")
>>> result.fe_params
Intercept    ...
x            ...
dtype: float64

Key parameters¶

Parameter	Type	Description
`formula`	`str`	Fixed-effects formula in Wilkinson notation (e.g. `"y ~ x + z"`)
`data`	`DataFrame`	Input data
`groups`	`str \| list[str]`	Column name(s) for crossed random intercepts (shorthand)
`random`	`list[str]`	lme4-style random-effect specs, e.g. `["(1 + x \| g)"]`
`method`	`"REML"` or `"ML"`	Estimator. Use `"ML"` for model comparison via LRT
`optimizer`	`"lbfgsb"` or `"bobyqa"`	Optimizer. `"bobyqa"` gives better R/lme4 parity
`weights`	`ndarray`	Observation-level prior weights, shape (n,)
`offset`	`ndarray`	Known term added to the linear predictor (not estimated)
`correlation`	`AR1 \| CompoundSymmetry`	Residual correlation structure for longitudinal data
`df_method`	`str`	`"satterthwaite"` (default) or `"kenward-roger"` for denominator DFs
`theta0`	`ndarray`	Initial variance parameters; defaults to ones

Examples¶

Random intercepts (shorthand)¶

import interlace

result = interlace.fit(
    "rt ~ condition",
    data=df,
    groups=["subject", "item"],  # crossed random intercepts
)

print(result.fe_params)           # fixed-effect coefficients
print(result.variance_components) # σ² per grouping factor
print(result.aic, result.bic)

Random slopes¶

result = interlace.fit(
    "rt ~ condition",
    data=df,
    random=[
        "(1 + condition | subject)",  # correlated intercept + slope
        "(1 | item)",                 # intercept only
    ],
)

# random_effects["subject"] is a DataFrame: one column per term
print(result.random_effects["subject"])
print(result.varcov)  # full random-effect covariance matrix

Model comparison with ML¶

# Fit with ML for likelihood ratio test
m1 = interlace.fit("y ~ x",      data=df, groups="g", method="ML")
m2 = interlace.fit("y ~ x + z",  data=df, groups="g", method="ML")

import scipy.stats
lrt_stat = 2 * (m2.llf - m1.llf)
p_value  = scipy.stats.chi2.sf(lrt_stat, df=1)