The same logistic 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_logistic 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", "study hours", "practice score"]):
        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	-3.170981	1.478959	-2.144062	0.032028	-6.069687	-0.272274
1	R	study hours	0.063675	0.113182	0.562592	0.573713	-0.158157	0.285508
2	R	practice score	0.034248	0.028941	1.183392	0.236654	-0.022474	0.090970
3	statsmodels	intercept	-3.170981	1.478959	-2.144062	0.032028	-6.069687	-0.272274
4	statsmodels	study hours	0.063675	0.113182	0.562592	0.573713	-0.158157	0.285508
5	statsmodels	practice score	0.034248	0.028941	1.183392	0.236654	-0.022474	0.090970
6	php-stats	intercept	-3.170981	1.478959	-2.144062	0.032028	-6.069687	-0.272274
7	php-stats	study hours	0.063675	0.113182	0.562592	0.573713	-0.158157	0.285508
8	php-stats	practice score	0.034248	0.028941	1.183392	0.236654	-0.022474	0.090970

	case	retained n	iterations (PHP)	max absolute difference vs R	warnings
0	simulated_learning	80	5	1.134662e-11	0
1	simple	80	5	1.767218e-12	0
2	no_intercept	80	5	3.753844e-12	0
3	intercept_only	80	4	6.502732e-13	0
4	correlated	80	5	1.284750e-12	0
5	imbalanced	80	8	3.527378e-08	0
6	listwise_missing	78	5	1.369999e-11	1

All complete statistical outputs, including every fitted probability 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/learning.csv")
fit <- glm(passed ~ study_hours + practice_score, data=students,
           family=binomial(), control=glm.control(epsilon=1e-14, maxit=100))
coef(fit)
vcov(fit)
confint.default(fit) # normal Wald, not profile likelihood

students = pd.read_csv("docs/data/learning.csv")
x = sm.add_constant(students[["study_hours", "practice_score"]])
fit = sm.Logit(students["passed"], x).fit(method="newton", tol=1e-10, maxiter=100)
fit.params, fit.bse, fit.conf_int()

$x = []; $y = [];
foreach ($students as $student) {
    $x[] = [1.0, $student['study_hours'], $student['practice_score']];
    $y[] = $student['passed'];
}
$result = stats_logistic(x: $x, y: $y);

Here $students contains already numeric records; see examples/logistic.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]
Dropped rows: [1, 4]
PHP warnings: [{'severity': 2, 'message': 'stats_logistic(): stats_logistic dropped 2 of 80 observations with missing values (listwise deletion)'}]

	implementation	log_likelihood	deviance	null_deviance	aic	df_residual
0	R	-49.484325	98.96865	104.774889	104.96865	77
1	statsmodels	-49.484325	98.96865	104.774889	104.96865	77
2	php-stats	-49.484325	98.96865	104.774889	104.96865	77

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