The same random-intercept ML model in R, Python and PHP

Equivalent model calls

library(lme4)
students <- read.csv('docs/data/schools.csv')
fit <- lmer(score ~ study_hours + prior_score + (1|school),
            data=students, REML=FALSE,
            control=lmerControl(optimizer='bobyqa',optCtrl=list(rhoend=1e-10)))
fixef(fit)
vcov(fit)
ranef(fit)

import statsmodels.api as sm
students = pd.read_csv('docs/data/schools.csv')
x = sm.add_constant(students[['study_hours','prior_score']])
fit = sm.MixedLM(students.score, x, groups=students.school).fit(
    reml=False, method='bfgs', gtol=1e-10, maxiter=1000)
fit.fe_params, fit.scale, fit.cov_re

$x = $y = $groups = [];
foreach ($students as $student) {
    $x[] = [1.0, $student['study_hours'], $student['prior_score']];
    $y[] = $student['score'];
    $groups[] = $student['school'];
}
$result = stats_lmm(x: $x, y: $y, groups: $groups);

These PHP records are already numeric; examples/lmm.php supplies the complete CSV conversion. No formulas or implicit intercepts are interpreted by the native call. The comparison runner gives R the exact design columns with formula y ~ 0 + x0 + ... + (1|group).

Conventions that must match

ML is explicit: R and statsmodels otherwise support/default to REML, which fits a different variance criterion. Likelihood here integrates over group effects. Fixed-effect standard errors use plug-in GLS covariance (X’ V^-1 X)^-1, matching lme4. statsmodels’ reported joint observed-Hessian covariance can differ, so the Python runner explicitly recalculates the declared GLS convention using its fitted variance components; it does not claim its unmodified bse_fe matches every case. Normal z-based Wald intervals are constructed explicitly for every reference.

Our df_residual=n-p is informational. R mixed-model df methods may include the two variance parameters; no finite-sample t degrees of freedom are claimed here. AIC counts p+2 parameters in every case, including a zero variance boundary.

The exact tau²=0 fixture is checked against lme4 and PHP, with a separately labeled NumPy OLS-limit reference. It is not represented as a converged statsmodels mixed fit. At this boundary, variance is RSS/n and standard errors differ from OLS’s unbiased-scale inference. This explicit check avoids treating an approximate positive-variance optimizer solution as proof of exact boundary behavior.

keys=['residual_variance','random_intercept_variance','icc','log_likelihood','aic']
display(pd.DataFrame([dict(implementation=name,**{k:first[name][k] for k in keys}) for name in ['R','Python','php-stats']]))
print('Python fit warnings:', first['python_diagnostics']['warnings'])
missing=results[-1]
print('Dropped original rows:',missing['php-stats']['dropped_rows'])
print('Missing-data warning:',missing['warnings'])
print('Schools and effects:',list(zip(first['php-stats']['group_labels'],first['php-stats']['random_effects'])))

Python fit warnings: []
Dropped original rows: [1, 4, 8]
Missing-data warning: [{'severity': 2, 'message': 'stats_lmm(): stats_lmm dropped 3 of 250 observations with missing values (listwise deletion)'}]
Schools and effects: [(107, 0.13413303204091775), (114, 1.6086975147685796), (121, -1.7846684515271147), (128, -5.564364724188796), (135, 0.24907018967277064), (142, 2.2564389505860314), (149, 10.331238035927292), (156, -5.104699176036904), (163, 2.102747547050482), (170, 0.5690811113270031), (177, 0.07903267727850147), (184, -4.87670670689903)]

See the verification record for rejection tests, finite-search limitations and reproduction instructions. A converged numerical fit does not establish that the model assumptions or small-group Wald inference are adequate for an application.