The same poisson models in R, statsmodels and PHP · phit

We fit the same saved observations and design columns in all three runtimes. See the math walkthrough for how likelihood, Newton updates and normal Wald inference produce these quantities.

from pathlib import Path
import sys
import pandas as pd
ROOT = next(p for p in [Path.cwd(), *Path.cwd().parents] if (p / "Cargo.toml").exists())
sys.path.insert(0, str(ROOT / "docs/verifications"))
from compare_poisson import cases, compare_case
results = [compare_case(case) for case in cases()]
first = results[0]
print(first["versions"])
rows = []
for implementation in ["R", "statsmodels", "php-stats"]:
    fit = first[implementation]
    for j, label in enumerate(["intercept", "visitors", "staff"]):
        rows.append({"implementation": implementation, "coefficient": label,
                     "estimate": fit["coefficients"][j], "SE": fit["std_errors"][j],
                     "z": fit["z_statistics"][j], "p": fit["p_values"][j],
                     "95% lower": fit["conf_int"][j][0], "95% upper": fit["conf_int"][j][1]})
display(pd.DataFrame(rows))
display(pd.DataFrame([{ "case": r["name"], "retained n": r["php-stats"]["n_observations"],
                        "iterations (PHP)": r["php-stats"]["iterations"],
                        "max absolute difference vs R": max(e for values in r["max_absolute_errors_vs_R"].values() for e in values.values()),
                        "warnings": len(r["warnings"])} for r in results]))

{'R': 'R version 4.6.1 (2026-06-24)', 'statsmodels': '0.15.0', 'numpy': '2.5.4', 'Python': '3.14.7', 'PHP': '8.5.10'}

	implementation	coefficient	estimate	SE	z	p	95% lower	95% upper
0	R	intercept	0.351754	0.243133	1.446754	0.147966	-0.124778	0.828286
1	R	visitors	0.005835	0.001839	3.172650	0.001511	0.002230	0.009440
2	R	staff	-0.234004	0.087265	-2.681524	0.007329	-0.405042	-0.062967
3	statsmodels	intercept	0.351754	0.243133	1.446754	0.147966	-0.124778	0.828286
4	statsmodels	visitors	0.005835	0.001839	3.172650	0.001511	0.002230	0.009440
5	statsmodels	staff	-0.234004	0.087265	-2.681524	0.007329	-0.405042	-0.062967
6	php-stats	intercept	0.351754	0.243133	1.446754	0.147966	-0.124778	0.828286
7	php-stats	visitors	0.005835	0.001839	3.172650	0.001511	0.002230	0.009440
8	php-stats	staff	-0.234004	0.087265	-2.681524	0.007329	-0.405042	-0.062967

	case	retained n	iterations (PHP)	max absolute difference vs R	warnings
0	simulated_clinic	120	6	6.528972e-12	0
1	simple	120	6	3.046319e-11	0
2	no_intercept	120	6	6.183248e-13	0
3	intercept_only	120	6	2.965983e-11	0
4	correlated	120	6	1.098599e-11	0
5	sparse	120	7	9.897534e-13	0
6	listwise_missing	118	6	5.690559e-12	1

All complete statistical outputs, including every fitted mean count and residual, are checked at tolerance 1e-9 + 1e-7 × abs(reference). Differences in iteration counts are expected because initial values and convergence rules differ.

Equivalent calls

R adds an intercept in this formula; the other two designs supply it explicitly.

students <- read.csv("docs/data/clinic.csv")
fit <- glm(calls ~ visitors + staff, data=students,
           family=poisson(), control=glm.control(epsilon=1e-16, maxit=100))
coef(fit)
vcov(fit)
confint.default(fit) # normal Wald, not profile likelihood

import statsmodels.api as sm
students = pd.read_csv("docs/data/clinic.csv")
x = sm.add_constant(students[["visitors", "staff"]])
fit = sm.GLM(students["calls"], x, family=sm.families.Poisson()).fit(tol=1e-11, tol_criterion="params", maxiter=100)
fit.params, fit.bse, fit.conf_int()

$x = []; $y = [];
foreach ($students as $student) {
    $x[] = [1.0, $student['visitors'], $student['staff']];
    $y[] = $student['calls'];
}
$result = stats_poisson(x: $x, y: $y);

Here $students contains already numeric records; see examples/poisson.php for the complete CSV conversion and call. The native interface uses no formula parser, implicit column construction or model classes.

Missing rows and model statistics

missing = results[-1]
print("Retained rows:", missing["php-stats"]["retained_rows"])
print("Dropped rows:", missing["php-stats"]["dropped_rows"])
print("PHP warnings:", missing["warnings"])
display(pd.DataFrame([{ "implementation": name,
                       **{key: first[name][key] for key in ["log_likelihood", "deviance", "null_deviance", "aic", "df_residual"]}}
                      for name in ["R", "statsmodels", "php-stats"]]))

Retained rows: [0, 2, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119]
Dropped rows: [1, 4]
PHP warnings: [{'severity': 2, 'message': 'stats_poisson(): stats_poisson dropped 2 of 120 observations with missing values (listwise deletion)'}]

	implementation	log_likelihood	deviance	null_deviance	aic	df_residual
0	R	-198.381638	144.530662	154.468908	402.763276	117
1	statsmodels	-198.381638	144.530662	154.468908	402.763276	117
2	php-stats	-198.381638	144.530662	154.468908	402.763276	117

The runner applies the identical retained sample before fitting reference models. For no-intercept cases it explicitly calculates the fixed mean-one null baseline, rather than using statsmodels’ intercept-only null fit. Deviance uses individual Poisson counts, and no Gaussian residual variance is estimated. See the verification record for rejection cases and limits.