Gaussian mixed models on four published datasets

Execute all four comparisons

from pathlib import Path
import json,sys
import numpy as np
import pandas as pd

from IPython.display import display
ROOT=next(p for p in [Path.cwd(),*Path.cwd().parents] if (p/'Cargo.toml').exists())
import sys
sys.path.insert(0,str(ROOT/'docs/figures'))
from bridge import render_figure
sys.path.insert(0,str(ROOT/'docs/verifications'))
from compare_lmm import public_cases,compare_case
cases=public_cases()
results=[compare_case(case,include_nlme=True) for case in cases]
ORDER=[('R','R: lme4'),('nlme','R: nlme'),('Python','statsmodels'),('php-stats','PHP')]
print('Versions:',results[0]['versions'])
display(pd.DataFrame([dict(dataset=c['name'],observations=len(c['y']),groups=len(set(c['groups'])),
                           fixed_coefficients=len(c['x'][0]),formula=c['source']['formula'],
                           all_comparisons_passed=True) for c in cases]))

Versions: {'R': 'R version 4.6.1 (2026-06-24)', 'lme4': '2.0.6', 'statsmodels': '0.15.0', 'PHP': '8.5.10', 'Python': '3.14.7', 'nlme': '3.1.169'}

The runner checks every row, not a sample: population and conditional fitted values, both residual definitions, all group effects, coefficients, covariance, standard errors, normal Wald tests/intervals, variance components, ICC, likelihood and AIC. Group labels, ordering, retained rows and observation counts are exact. Both R packages supply independent model estimates. statsmodels supplies its own ML estimates; its fixed-effect covariance is recomputed as plug-in GLS to match lme4 and our contract. All Python mixed fits must converge; warnings are retained. Both R references use vcov() and explicit normal Wald calculations to match our declared inference. Default nlme t-tests and profile-likelihood intervals use different conventions.

Estimates by dataset

Variances and likelihood

rows=[]
for result in results:
    for key,label in ORDER:
        fit=result[key]
        rows.append(dict(dataset=result['name'],implementation=label,
                         residual_variance=fit['residual_variance'],group_variance=fit['random_intercept_variance'],
                         ICC=fit['icc'],log_likelihood=fit['log_likelihood'],AIC=fit['aic']))
display(pd.DataFrame(rows))

Look closely at numerical differences

The graph shows the difference between native and lme4 conditional fitted values, scaled by a million to make rounding differences visible. Zero means exact agreement. Every point is a retained observation; panel units differ by dataset. It is a numerical agreement plot, not a model residual diagnostic.

figures/lmm-public-agreement.svg 

Default public-data tolerances are 1e-7 absolute plus 1e-6 relative. The four observation-length output vectors use 1e-6 absolute in their outcome’s units, plus 1e-6 relative, so residuals close to zero are not tested by relative error alone. Likelihood and AIC use 1e-7 absolute plus 1e-10 relative. Original synthetic fixtures retain their uniform 1e-7 + 1e-6 × abs(reference) tolerance. These policies are recorded with each source fixture and live report. An additional native-versus-statsmodels check keeps the stricter original uniform tolerance for all quantities; reference-optimizer rounding allowances cannot hide divergence between those two estimates.

The High School and Beyond sample illustrates an optimizer precision limit: lme4 and nlme have effectively identical likelihoods but slightly different variance estimates. Conditional residual differences can reach a few 1e-7 when the residual itself is nearly zero. The profile scores below independently check stationarity at each estimate; they are derivatives of negative profile log-likelihood with respect to log(1 + tau²/sigma²). Zero is the interior optimum condition. Raw scores depend on dataset size and curvature and should not be compared across different studies as a universal error metric.

display(pd.DataFrame([dict(dataset=r['name'],implementation=label,
                           profile_score=r['profile_scores'][key],
                           log_likelihood_difference_vs_lme4=r[key]['log_likelihood']-r['R']['log_likelihood'])
                      for r in results for key,label in ORDER]))
for result in results:
    print(result['name'],'Python warnings:',result['python_diagnostics']['warnings'])

sleepstudy Python warnings: []
Orthodont Python warnings: []
Machines Python warnings: []
Hsb82 Python warnings: []

Reproduce without fetching data

The CSVs and lme4/nlme expected outputs are checked into the repository. cargo test --workspace and tests/lmm.php check the published fixtures without R, Python or network access. The live comparison requires lme4, nlme and the pinned Python environment, but not mlmRev: it reads the saved Hsb82 CSV.

cargo test --workspace --locked
PHP=/opt/homebrew/bin/php STATS_EXTENSION="$PWD/target/release/libphp_stats.dylib" \
  .context/notebooks-venv/bin/python docs/verifications/compare_lmm.py --public
QUARTO_PYTHON="$PWD/.context/notebooks-venv/bin/python" PHP=/opt/homebrew/bin/php \
  quarto render docs/notebooks/lmm-public-verification.qmd

Use matching PHP headers to build the extension first; use .so on Linux. The saved live report records full outputs, source checksums, runtime versions, the extension hash, warnings, tolerances and score/error checks. Re-execute this page after statistical-code or source-data changes; merely rendering saved outputs is not new verification.

This strengthens evidence for the current Gaussian random-intercept ML model. It does not establish correctness for random slopes, REML, other response families, assumption checks, or finite-sample group-level inference.