Expanded R compatibility and all four lme4 optimizers

Run the audit

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_r_compatibility import run
report=run()
print(report['versions'])

Passed 12 fixed-model cases and 15 mixed-model cases × 2 native × 4 lme4 optimizers; sleepstudy ML/REML also passed statsmodels.
{'R': 'R version 4.6.1 (2026-06-24)', 'lme4': '2.0.6', 'PHP': '8.5.10', 'python': '3.14.7', 'numpy': '2.5.4', 'scipy': '1.18.1', 'statsmodels': '0.15.0'}

The audit checks every coefficient, inference quantity, fitted value and residual; for mixed models it also checks every random effect and covariance entry, variance scale, likelihood, AIC/BIC and singular flag. It refuses nonzero R optimizer codes or R warnings in the saved reference runs. Reference controls and checksums are recorded alongside the actual PHP outputs in the downloadable audit.

Fixed estimators: weights, offsets and binomial trials

display(pd.DataFrame([dict(case=r['name'],estimator=r['model'],
                          largest_absolute_error=r['max_absolute_error'],passed=True)
                     for r in report['fixed']]))

Each estimator is fitted with an intercept and without one. Binary estimators include both individual binary responses and grouped success proportions. Weights and offsets affect the fit, covariance, null model and likelihood. We check those quantities rather than testing only unweighted coefficients.

For grouped binomial rows with s successes out of m trials, supply y=s/m and weight=m. R’s equivalent calls are:

#| eval: false
fit <- glm(proportion ~ study + practice,
           family=binomial("logit"), weights=trials, offset=known_offset,
           data=observations, control=glm.control(epsilon=1e-12, maxit=200))
# binomial("probit") selects the normal-CDF link instead.
coef(fit)
vcov(fit)
logLik(fit)
fit$deviance
fit$null.deviance

A Poisson exposure e contributes log(e) to the linear predictor; a prior weight instead multiplies that observation’s likelihood. Those are different operations. For weighted OLS, residual variance is sigma²/w_i and likelihood includes the weight determinant. Full generator code is in tests/generate-r-compatibility-fixtures.R.

Mixed models and optimizer agreement

R_ORDER=['bobyqa','Nelder_Mead','nlminbwrap','nloptwrap']
rows=[]
for case in report['mixed']:
    for optimizer in R_ORDER:
        ref=case['R'][optimizer]
        for comparison in case['comparisons']:
            if comparison['R_optimizer']==optimizer:
                rows.append(dict(model=case['name'],method=case['method'],
                                 R_optimizer=optimizer,native_optimizer=comparison['native_optimizer'],
                                 log_likelihood_difference=comparison['log_likelihood_difference'],
                                 max_fitted_difference=comparison['max_fitted_difference'],
                                 optimizer_code=ref['optimizer_code'][0]))
display(pd.DataFrame(rows))

The model suite contains correlated/independent random slopes, crossed and nested terms, a random slope without a random intercept, no fixed intercept, weights and offsets, and an exact zero-variance boundary. The published sleepstudy model is Reaction ~ Days + (Days | Subject), under both ML and REML.

The R generator explicitly selects all four wrappers, tightened stopping controls and boundary.tol=1e-3. The exact per-wrapper settings are recorded in each fixture. Multi-level BOBYQA references also restart from their first R fit with rhobeg=0.001 and rhoend=1e-9: this improved random-effect agreement on a flat likelihood while keeping the comparison tolerances unchanged. Native estimates never seed the R references. At a smaller default boundary tolerance, one wrapper stopped just above a singular boundary; overly strict nlminb settings also produced convergence warnings. The saved audit does not claim identical behavior under every package default.

Likelihoods use 5e-6 + 1e-9 * abs(reference) tolerance. Native/R BOBYQA estimates use 1e-5 + 1e-5 * abs(reference); other optimizer comparisons use 1e-3 + 3e-4 * abs(reference). These admit the observed spread in parameters on a flat likelihood. Shapes, grouping/sample metadata and parameter counts are exact. All rows and raw outputs remain available for inspection.

More than two grouping levels

The fitter accepts a list of random terms; there is no two-level cap. The new synthetic cases use 256 repeated assessments, with 8 districts, 16 schools, 32 classes and 64 students. They fit three grouping levels (district/school/class) and four grouping levels (district/school/class/student), each under ML and REML. The four-level simulation includes independent Gaussian intercepts at all four grouping levels; the three-level simulation omits the student effect. Study hours and practice scores are correlated observed predictors, not initial parameter guesses.

fixture=json.loads((ROOT/'tests/fixtures/lmm-expanded-r.json').read_text())
nested_inputs={c['name']:c for c in fixture['cases'] if '_grouping_levels_' in c['name']}
summary=[]
for c in report['mixed']:
    if c['name'] not in nested_inputs:continue
    input_case=nested_inputs[c['name']]
    for optimizer in ['bobyqa','nelder_mead']:
        fit=c['PHP'][optimizer]
        summary.append(dict(model=c['name'],optimizer=optimizer,
                            groups_per_term=[len(t['group_labels']) for t in fit['random']],
                            estimated_variances=[t['covariance'][0][0] for t in fit['random']],
                            residual_variance=fit['residual_variance'],log_likelihood=fit['log_likelihood'],
                            all_four_R_comparisons_passed=True))
display(pd.DataFrame(summary))

Each random term is independent of the other terms. Dependence between observations arises because they share some group effects. For two observations in the same class but different students, the district, school and class effects are shared; the student effects and observation noise are not.

Provide globally unique nested IDs: the native interface does not infer that a class label belongs to a school. Reusing local class number 1 across schools without a parent-aware key instead specifies a crossed/shared group effect. Three grouping levels plus the observation level are often called a four-level model; we count grouping terms explicitly here to avoid that ambiguity.

$result = stats_lmm_fit(
    x: $fixedDesign,
    y: $scores,
    random: [
        ['z' => $ones, 'groups' => $districtIds],
        ['z' => $ones, 'groups' => $districtSchoolKeys],
        ['z' => $ones, 'groups' => $districtSchoolClassKeys],
        ['z' => $ones, 'groups' => $districtSchoolClassStudentKeys],
    ],
    reml: true,
);

The equivalent R random-intercept specification is (1 | district/school/class/student). Nesting does not introduce a correlation parameter between the levels; each term has its own variance. Random slopes can also be supplied through each term’s z columns, subject to identifiability and size limits. The current dense implementation supports at most 512 total random coefficients and 50,000,000 initial augmented matrix entries.

Published sleepstudy: R, statsmodels and PHP

for case in [r for r in report['mixed'] if r['name'].startswith('sleepstudy_')]:
    fits=[('R: '+o,case['R'][o]['result']) for o in R_ORDER]
    fits+=[('statsmodels',case['statsmodels'])]
    fits+=[('PHP: '+o,case['PHP'][o]) for o in ['bobyqa','nelder_mead']]
    print(case['name'])
    display(pd.DataFrame([dict(implementation=label,intercept=fit['coefficients'][0],
                              days=fit['coefficients'][1],residual_variance=fit['residual_variance'],
                              intercept_variance=fit['terms'][0]['covariance'][0][0],
                              slope_variance=fit['terms'][0]['covariance'][1][1],
                              intercept_slope_covariance=fit['terms'][0]['covariance'][0][1],
                              log_likelihood=fit['log_likelihood']) for label,fit in fits]))
    print('statsmodels warnings:',case['statsmodels_warnings'])

sleepstudy_ml

statsmodels warnings: []
sleepstudy_reml

statsmodels warnings: []

statsmodels independently estimates ML/REML fixed effects and variance components. For fixed-effect inference we calculate plug-in GLS covariance from its estimates, matching lme4 and our contract; its default joint-Hessian standard errors differ. We require convergence and retain Python warnings. BFGS uses gtol=1e-8; at 1e-10 the REML fit reached the same likelihood but failed its convergence flag, so that failed run was not accepted. Its public MixedLM interface does not provide prior weights, so this Python comparison covers sleepstudy, not the weighted case or every general grouping specification.

What the random slopes mean

data=pd.read_csv(ROOT/'docs/data/lmm/sleepstudy.csv')
case=next(r for r in report['mixed'] if r['name']=='sleepstudy_reml')
fit=case['PHP']['bobyqa'];labels=fit['random'][0]['group_labels']
fixture=next(c for c in json.loads((ROOT/'tests/fixtures/lmm-expanded-r.json').read_text())['cases'] if c['name']=='sleepstudy_reml')
subject_group=dict(zip(data.Subject,fixture['terms'][0]['groups']))
# Fixture IDs are factor-coded by the R generator; original subject ordering is
# retained in the CSV. Map the groups explicitly rather than guessing labels.
subjects=sorted(data.Subject.unique())
render_figure('sleepstudy')

figures/sleepstudy-slopes.svg 

The observed predictor is days of sleep deprivation; it is not a trial coefficient. Each person’s intercept shifts the line, and their random slope changes its rate of increase. The population line averages over the zero-mean random effects. This model captures individual differences; numerical agreement alone does not establish that the Gaussian assumptions explain the study.

Reproduce or refresh

# Normal audit: stored R values, fresh PHP and statsmodels fits.
PHP=/opt/homebrew/bin/php .context/notebooks-venv/bin/python docs/verifications/compare_r_compatibility.py
# Explicitly refresh reference fits using installed R/lme4.
PHP=/opt/homebrew/bin/php .context/notebooks-venv/bin/python docs/verifications/compare_r_compatibility.py --refresh
QUARTO_PYTHON="$PWD/.context/notebooks-venv/bin/python" quarto render docs/notebooks/r-compatibility.qmd