Linear mixed models

A Gaussian linear mixed model (LMM) combines population-level fixed effects with latent group-specific deviations.

\[ y=X\beta+Zb+\epsilon,\qquad b\sim N(0,G),\qquad \epsilon\sim N(0,R). \]

The marginal covariance is \(V=ZGZ^\mathsf{T}+R\). PyMixEF profiles \(\beta\), evaluates the normalized Gaussian likelihood with stable Cholesky solves, and optimizes unconstrained covariance parameters that map to positive-definite natural-scale matrices.

When to use an LMM

An LMM is a useful starting point when:

  • the conditional response is reasonably modeled as continuous Gaussian;

  • observations share groups such as subject, batch, center, lot, device, or laboratory;

  • random intercepts/slopes are a scientifically meaningful dependence model;

  • the target includes population coefficients, group-conditional prediction, variance components, or trajectories.

For repeated visits whose residual covariance is the primary object, compare the MMRM path.

Declare a model

import pymixef

model = pymixef.Model.from_formula(
    "response ~ treatment * time + (1 + time | subject) + (1 | center)",
    family=pymixef.families.Gaussian(),
)

plan = model.compile(data, engine="lmm", method="reml")
print(plan.explain())
result = plan.fit()

Use || for independent random coefficients:

response ~ treatment * time + (1 + time || subject)

Multiple grouping blocks are supported. Their covariance parameters, levels, and design columns are compiled explicitly.

The declaration enters through pymixef.model.Model.from_formula(), compiles to an pymixef.model.ExecutionPlan, and returns a pymixef.results.FitResult.

ML and REML

Under ML, the normalized log likelihood is

\[ \ell(\beta,\theta)=-\frac12 \left[n\log(2\pi)+\log|V_\theta| +(y-X\beta)^\mathsf{T}V_\theta^{-1}(y-X\beta)\right]. \]

REML additionally accounts for estimating fixed effects by integrating their linear space under the package’s archived constant convention:

\[ \ell_\mathrm{REML}(\theta) =-\frac12\left[ (n-p)\log(2\pi)+\log|V_\theta| +\log|X^\mathsf{T}V_\theta^{-1}X| +y^\mathsf{T}P_\theta y \right], \qquad P_\theta=V_\theta^{-1} -V_\theta^{-1}X(X^\mathsf{T}V_\theta^{-1}X)^{-1}X^\mathsf{T}V_\theta^{-1}. \]

  • ML is normally used when comparing nested fixed-effect mean structures.

  • REML is commonly used for covariance estimation when the fixed-effect design is held constant.

  • Do not compare REML objectives from different fixed-effect spaces as if they shared the same likelihood.

What the fit returns

FitResult.parameters contains fixed coefficients and natural-scale covariance summaries such as random-effect SD/correlation and residual SD. FitResult.extra can include:

  • fixed_effect_names and fixed_effect_covariance;

  • objective-convention and likelihood-constant metadata;

  • compiled design/audit summaries;

  • covariance matrices and optimizer provenance.

FitResult.random_effects contains conditional modes where supported. FitResult.fitted_values and FitResult.residuals align to analysis rows.

Prediction modes

conditional = result.prediction(mode="conditional")
population = result.prediction(mode="population")

  • Conditional prediction includes estimated random effects and generally fits observed groups more closely.

  • Population prediction uses fixed effects at the random-effect reference value.

  • Prediction for new groups requires an explicit new-random-effect simulation or integration path; it is not the same as reusing observed group modes.

Both modes are selected explicitly through FitResult.prediction.

Variance and random-effect interpretation

Random-effect SDs describe modeled between-group heterogeneity. A conditional mode is pulled toward zero according to group information and covariance; it is not a directly observed group effect. Near-zero variance or near-perfect correlation can indicate a meaningful boundary, weak information, or an overly rich random structure.

Inspect ConvergenceReport.boundaries, the pymixef.convergence.HessianDiagnostics, scaled gradient, and the covariance pymixef.covariance.singularity_report() before reporting.

Diagnostics

At minimum:

  1. confirm ConvergenceReport.trustworthy;

  2. reconcile source and analysis rows;

  3. inspect conditional and population fits;

  4. plot residuals against fitted values, time, and important predictors;

  5. inspect group sizes and random-mode shrinkage;

  6. review covariance boundaries/conditioning;

  7. run seeded simulation or bootstrap when it answers a prespecified question;

  8. test scientifically plausible mean/covariance/missingness sensitivities.

Reference-engine scope

The initial backend is a dense experimental reference engine intended for small/moderate correctness experiments. It is not the planned compiled sparse million-row backend. Passing unit tests or obtaining a trustworthy fit does not constitute external-software parity for every design.

Worked examples

API