The same probit 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, Fisher-scoring 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_probit 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	-1.251649	0.704093	-1.777676	0.075457	-2.631645	0.128347
1	R	study hours	0.150439	0.061547	2.444279	0.014514	0.029808	0.271070
2	R	practice score	-0.000054	0.014593	-0.003682	0.997062	-0.028655	0.028547
3	statsmodels	intercept	-1.251649	0.704093	-1.777676	0.075457	-2.631645	0.128347
4	statsmodels	study hours	0.150439	0.061547	2.444279	0.014514	0.029808	0.271070
5	statsmodels	practice score	-0.000054	0.014593	-0.003682	0.997062	-0.028655	0.028547
6	php-stats	intercept	-1.251649	0.704093	-1.777676	0.075457	-2.631645	0.128347
7	php-stats	study hours	0.150439	0.061547	2.444279	0.014514	0.029808	0.271070
8	php-stats	practice score	-0.000054	0.014593	-0.003682	0.997062	-0.028655	0.028547

	case	retained n	iterations (PHP)	max absolute difference vs R	warnings
0	simulated_learning	120	9	6.375807e-10	0
1	simple	120	8	2.718537e-10	0
2	no_intercept	120	8	5.628209e-11	0
3	intercept_only	120	4	1.999630e-12	0
4	correlated	120	10	2.713989e-10	0
5	imbalanced	120	16	1.369056e-10	0
6	listwise_missing	118	10	3.723812e-10	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/probit-learning.csv")
fit <- glm(passed ~ study_hours + practice_score, data=students,
           family=binomial(link="probit"), 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/probit-learning.csv")
x = sm.add_constant(students[["study_hours", "practice_score"]])
fit = sm.GLM(students["passed"], x, family=sm.families.Binomial(link=sm.families.links.Probit())).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['study_hours'], $student['practice_score']];
    $y[] = $student['passed'];
}
$result = stats_probit(x: $x, y: $y);

Here $students contains already numeric records; see examples/probit.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_probit(): stats_probit dropped 2 of 120 observations with missing values (listwise deletion)'}]

	implementation	log_likelihood	deviance	null_deviance	aic	df_residual
0	R	-72.927605	145.85521	162.299087	151.85521	117
1	statsmodels	-72.927605	145.85521	162.299087	151.85521	117
2	php-stats	-72.927605	145.85521	162.299087	151.85521	117

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.

Probit comparisons use expected-information covariance from GLM, rather than the observed-Hessian covariance of statsmodels discrete Probit.