OLS verification in R, statsmodels and PHP

Reproduction

Follow notebook setup to install development tools and build the PHP extension, then run:

.context/notebooks-venv/bin/python docs/verifications/compare_ols.py

The script fits every case afresh using Rscript, statsmodels and PHP, checks all statistical fields against R with absolute tolerance 1e-9 and relative tolerance 1e-7, and writes ols-results.json. It consumes the ordinary deterministic cases from tests/fixtures/ols-r.json plus an explicit missing-data case. R and Python are optional development tools; cargo test and tests/ols.php still run without them. See verification policy for the tolerance’s limits.

Equivalent specifications

For the shared simulated housing dataset in examples/ols.php, run from the repository root:

observations <- read.csv("docs/data/housing.csv")
fit <- lm(monthly_rent ~ area_m2, data = observations, na.action = na.omit)
summary(fit)
confint(fit, level = 0.95)
logLik(fit)
# A second model uses correlated predictors:
summary(lm(monthly_rent ~ area_m2 + bedrooms, data = observations))

import pandas as pd
import statsmodels.api as sm
observations = pd.read_csv("docs/data/housing.csv")
X = sm.add_constant(observations[["area_m2"]])
y = observations["monthly_rent"]
fit = sm.OLS(y, X, missing='drop', hasconst=True).fit(use_t=True)
print(fit.params, fit.bse, fit.tvalues, fit.pvalues, fit.conf_int(), fit.llf)

// The verification runner prepares a complete design from the shared housing CSV.
// Each row of $x is [1.0, area_m2]; $y contains monthly_rent.
$result = Phit\fit_ols(x: $x, y: $y, intercept: false);

These are synthetic apartments, not empirical housing records. The math notebook explains the generating process and reports sample correlations. The fixtures retain small hand-checkable numerical cases alongside these teaching examples.

These use an intercept, conventional homoskedastic covariance, residual variance RSS/(n-p), Student t tests, and 95% intervals. Likelihood uses ML variance RSS/n. The general R runner receives the explicit design columns with formula y ~ . - 1, then calculates R-squared using the same ones-column convention as the extension; this avoids a duplicate intercept. R’s coefficients, covariance, standard errors, tests, intervals, residuals and log-likelihood come from lm(). The notebook displays both the user-facing calls and this runner convention.

Sources: R lm, statsmodels OLS and statsmodels missing data.

Verified cases

All fields pass the stated tolerance for all seven cases. The executable notebook displays the housing coefficients and inference from all three implementations. The saved report records the actual environment: R 4.6.1, statsmodels 0.15.0, NumPy 2.5.4, Python 3.14, and PHP 8.5.10 on the development Mac. This is not evidence of testing Linux or other PHP builds.

For the missing-data case, original positions 1 and 3 are omitted; positions 0, 2, 4 and 5 are retained. PHP emits one E_WARNING with the dropped/original counts and returns sample metadata. The runner checks those positions and counts, and exposes the captured warning in both the saved report and notebook.

Perfect/zero fits and invalid/singular inputs are covered by the Rust/PHP tests, not finite three-way comparisons. Decomposition-dependent tiny residuals make perfect-fit t/p and likelihood comparisons unstable. References may also accept singular designs by aliasing coefficients while our implementation rejects them.