Generalized linear mixed models¶

A GLMM combines a non-Gaussian response family, link function, fixed effects, and Gaussian random effects:

\[ g\{\operatorname{E}(Y_{ij}\mid b_i)\} =x_{ij}^{T}\beta+z_{ij}^{T}b_i,\qquad b_i\sim N(0,G). \]

PyMixEF’s current reference engine integrates supported random effects with a first-order Laplace approximation around the joint conditional mode.

Supported fit subset¶

Family	Typical endpoint	Link/parameterization
`Bernoulli`	one binary outcome per row	supported canonical logit path
`Binomial`	successes out of trials	validated supported binomial path
`Poisson`	count	supported canonical log path
`NegativeBinomial2`	overdispersed count	supported NB2 mean/dispersion path

The probability family catalog is broader. Catalog objects such as beta, ordinal, zero-inflated, hurdle, survival, NB1, and COM-Poisson can support standalone calculations or model representation without being accepted by this engine.

Fit a Bernoulli random-intercept model¶

import pymixef

model = pymixef.Model.from_formula(
    "response ~ treatment + baseline + (1 | clinic)",
    family=pymixef.families.Bernoulli(),
)

validation = model.validate(data, engine="glmm", method="laplace")
validation.raise_for_errors()

plan = model.compile(
    data,
    engine="glmm",
    method="laplace",
    maxiter=400,
)
print(plan.explain())
result = plan.fit()

pymixef.model.Model.validate() returns a structured pymixef.model.ValidationReport; its raise_for_errors() method prevents an invalid declaration from reaching pymixef.model.Model.compile().

Laplace calculation¶

For population parameters \(\theta\), the engine finds the random-effect conditional mode and uses local curvature to approximate the random-effect integral. The result retains:

approximation name: first-order Laplace at the joint conditional mode;
quadrature order 1 as an explicit calculation label;
normalization/objective convention;
inner conditional-mode failure count;
outer optimizer state, gradient, boundaries, and warnings.

Adaptive Gauss–Hermite quadrature (AGHQ) is unavailable. Requesting aghq raises a stable compatibility error; Laplace is not silently substituted.

Writing the joint log density as \(h_i(b_i;\theta)=\log p(y_i\mid b_i,\theta)+\log p(b_i\mid\theta)\), the first-order Laplace contribution is

\[ \log p(y_i\mid\theta) \approx h_i(\widehat b_i;\theta) +\frac{q_i}{2}\log(2\pi) -\frac12\log\left|-\nabla^2_{b_i}h_i(\widehat b_i;\theta)\right|. \]

The reported approximation metadata on FitResult.extra makes this first-order calculation auditable.

Interpret coefficients on the correct scale¶

For a logit model, a coefficient is a conditional log odds ratio:

import numpy as np

conditional_odds_ratio = np.exp(result.parameters["treatment"])

It is not a risk ratio, probability ratio, or percentage-point change. A probability curve with the random effect set to zero describes a typical cluster conditional on \(b=0\), not a population-marginal curve integrated over the random-effect distribution.

For a log link, exponentiated coefficients are conditional multiplicative effects on the mean.

Random-effect modes¶

Conditional modes are shrunk estimates. Extreme observed cluster rates can produce moderated modes because the model pools information. Modes should not be treated as permanent clinic, lot, or facility performance scores; their uncertainty, exposure, case mix, and estimation on the same outcomes matter.

Diagnostics¶

Review both optimization levels:

outer parameter optimization;
cluster-level conditional-mode convergence;
scaled gradient and warning codes;
covariance boundary/singularity;
fitted probabilities/means versus observed summaries;
residual diagnostics appropriate to the endpoint;
cluster-level influence and held-out-cluster performance;
approximation sensitivity for consequential conclusions.

Binned calibration computed on the same data and fitted conditional modes is an in-sample diagnostic. It is not external predictive validation.

Simulation and predictions¶

FitResult.simulate can generate seeded replicate data under archived fixed and random-effect covariance. Preserve whether new random effects and residual sampling were included. Simulation checks consistency with the fitted model; it does not prove that the data-generating assumptions are true.

Reference-engine scope¶

The engine is a dense experimental first-order Laplace implementation. It does not claim AGHQ order sensitivity, sparse production scaling, zero-inflation or hurdle fitting, noncanonical-link generality, or a full glmmTMB-class feature surface.

Worked examples¶

Both show conditional probability curves, calibration, random modes, exact executed results, and explicit interpretation boundaries.