Binomial mixed-model verification and approximation accuracy

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

Passed 28 R/PHP comparisons, independent likelihood/integration checks, and statsmodels zero-variance reduction.
Independent cbpp marginal ML fit: {'dataset': 'cbpp', 'coefficients': [-1.399230360516453, -0.9914039350112828, -1.1278198196624616, -1.579470750553407], 'random_intercept_sd': 0.6475184370023838, 'log_likelihood': -91.98336903731132, 'converged': True, 'evaluations': 120, 'message': 'CONVERGENCE: RELATIVE REDUCTION OF F <= FACTR*EPSMCH', 'method': 'SciPy L-BFGS-B + QUADPACK integration', 'start': 'statsmodels binomial GLM, SD=1'}
{'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'}

Seven input datasets are fitted with 1, 9, 25 and 61 nodes: 28 comparisons. They cover binary outcomes, grouped counts with offsets, no fixed intercept, small/unbalanced groups with rare events, large group variation, a zero random variance boundary, and the published cbpp model. R references use BOBYQA with rhoend=1e-9 and tolPwrss=1e-12. Codes must be zero and saved warnings empty. Rust/PHP tests compare every coefficient, covariance/inference value, observation prediction/residual, conditional group mode, scale, normalized likelihood and AIC/BIC.

The independent Python evaluator uses Brent root finding for conditional modes, NumPy Hermite nodes for the same approximation, and QUADPACK adaptive integration for a separate numerical accuracy check. Its matched-approximation likelihood must agree with native output within 2e-8 absolute. QUADPACK error estimates are recorded; they are numerical estimates, not mathematical error guarantees.

Published cbpp: matching estimation methods

rows=[]
for r in report['results']:
    if r['dataset']!='cbpp':continue
    for implementation,fit in [('R: lme4',r['R']),('PHP',r['PHP'])]:
        rows.append(dict(nodes=r['quadrature_points'],implementation=implementation,
                         intercept=fit['coefficients'][0],group_SD=fit['random_intercept_sd'],
                         log_likelihood=fit['log_likelihood'],AIC=fit['aic']))
display(pd.DataFrame(rows))
print('Independent marginal ML fit using SciPy optimization and QUADPACK:',report['exact_fit'])

Independent marginal ML fit using SciPy optimization and QUADPACK: {'dataset': 'cbpp', 'coefficients': [-1.399230360516453, -0.9914039350112828, -1.1278198196624616, -1.579470750553407], 'random_intercept_sd': 0.6475184370023838, 'log_likelihood': -91.98336903731132, 'converged': True, 'evaluations': 120, 'message': 'CONVERGENCE: RELATIVE REDUCTION OF F <= FACTR*EPSMCH', 'method': 'SciPy L-BFGS-B + QUADPACK integration', 'start': 'statsmodels binomial GLM, SD=1'}

The independent SciPy optimizer starts from a statsmodels binomial GLM and SD=1, not from native mixed estimates. It optimizes the numerically integrated cbpp likelihood and checks its estimates against high-node native output.

statsmodels BinomialBayesMixedGLM uses Bayesian posterior estimation with priors on fixed effects and variance parameters. That is a different objective, so we do not claim equivalent mixed ML estimates from it. statsmodels GLM verifies the exact zero-random-variance reduction instead. The full mixed-model Python reference here is explicitly SciPy.

Likelihood normalization matters

display(pd.DataFrame([dict(nodes=r['quadrature_points'],
                         R_reported_log_likelihood=r['R_reference']['reported_log_likelihood'],
                         R_additive_correction=r['R_reference']['log_likelihood_correction'],
                         R_normalized_log_likelihood=r['R']['log_likelihood'],
                         PHP_normalized_log_likelihood=r['PHP']['log_likelihood'])
                    for r in report['results'] if r['dataset']=='cbpp']))

For grouped-binomial responses, lme4 nAGQ>1 computes its reported criterion using conditional deviance relative to the saturated model. Its additive constant differs from nAGQ=1. We retain the raw R value and restore the saturated log likelihood, sum(dbinom(successes, trials, successes/trials, log=TRUE)), before comparing normalized likelihoods and AIC/BIC. For binary observations this constant is zero. This convention is disclosed in the glmer documentation.

The native likelihood always includes binomial coefficients and Gaussian normalization. Changing node counts does not change its additive convention.

Approximation error, separately from implementation agreement

accuracy=pd.DataFrame([dict(dataset=r['dataset'],nodes=r['quadrature_points'],
                           approximation_error=r['approximation_log_likelihood_error'],
                           group_SD=r['PHP']['random_intercept_sd'],
                           max_coefficient_difference_vs_R=r['max_coefficient_difference'])
                      for r in report['results']])
display(accuracy)
render_figure('approximation')

figures/logistic-mixed-integration-error.svg 
figures/logistic-mixed-node-sensitivity.svg 

The likelihood errors are measured at each approximation’s own fitted parameter vector, using independent integration. They are not differences between log likelihoods at different parameters. The rare-event fixture makes the distinction visible: Laplace agrees with R but substantially changes estimated group variance. Higher node counts improve this example; 61 nodes is still an approximation. Check estimate sensitivity to node counts when groups are small or events rare.

Numerical tolerance and failure policy

• Estimates/predictions/group modes: absolute 2e-5 plus relative 2e-4.
• Covariance, standard errors, z tests, p-values and intervals: absolute 2e-5 plus relative 2e-3; R and native use independently computed finite-difference Hessians.
• Likelihood/AIC/BIC: absolute 1e-6 plus relative 1e-8.
• Independent matched-approximation likelihood: absolute 2e-8.
• Sample metadata, group ordering, integer counts and singular flags are exact where checked. Boundary reduction also checks equality with the native GLM.

The coefficient covariance is the fixed block of inverse joint marginal observed information, including scale coupling. A nonpositive Hessian, excessive final score, invalid/aliased design, fixed separation, invalid counts or exhausted budget throws an error. Controls have explicit parameter bounds; convergence is local and does not prove global optimality. Missingness follows the existing listwise policy with one warning and original retained/dropped positions.

Reproduce

cargo test --workspace --locked
cargo build -p php-stats --release --locked
php -n -d "extension=$PWD/target/release/libphp_stats.dylib" tests/logistic-mixed.php
php -n -d "extension=$PWD/target/release/libphp_stats.dylib" examples/logistic-mixed.php
PHP=/opt/homebrew/bin/php .context/notebooks-venv/bin/python docs/verifications/compare_logistic_mixed.py
# Explicit reference refresh, not needed by ordinary Rust/PHP tests:
Rscript tests/generate-logistic-mixed-fixtures.R

Use .so on Linux and select matching PHP/php-config as in the development guide. The saved full audit includes versions, fixture/data/extension hashes, original R likelihoods and corrections, all actual fit outputs, integration error measurements and the independent SciPy fit.