OLS verification: R, statsmodels and PHP

Equivalent calls

R, using the table above:

fit <- lm(monthly_rent ~ area_m2, data = observations, na.action = na.omit)
summary(fit)
confint(fit, level = 0.95)
logLik(fit)

Python statsmodels:

fit = sm.OLS(y, X, missing='drop', hasconst=True).fit(use_t=True)
fit.params, fit.bse, fit.tvalues, fit.pvalues, fit.conf_int(), fit.llf

PHP (after loading the extension):

$result = stats_ols(x: $x, y: $y); // missing defaults to 'drop'

The executable harness uses y ~ . - 1 in R with the explicit design columns to avoid inserting a second intercept. It calculates centered/uncentered R-squared explicitly because that formula is syntactically no-intercept even when the supplied matrix contains ones. All other quantities come from R’s fitted model. This makes the same runner work for multiple and no-intercept designs.

simple_result = compare_case(simple)
print(simple_result['versions'])
rows = []
for name in ['R', 'statsmodels', 'php-stats']:
    result = simple_result[name]
    for i, label in enumerate(['intercept', 'area_m2']):
        rows.append({'implementation': name, 'term': label, 'coefficient': result['coefficients'][i],
                     'SE': result['std_errors'][i], 't': result['t_statistics'][i], 'p': result['p_values'][i],
                     'lower 95%': result['conf_int'][i][0], 'upper 95%': result['conf_int'][i][1]})
display(pd.DataFrame(rows))
display(pd.DataFrame([{ 'implementation': name, **{key: simple_result[name][key] for key in
                     ['rss', 'sigma', 'df_residual', 'r_squared', 'adj_r_squared', 'log_likelihood']}}
                     for name in ['R', 'statsmodels', 'php-stats']]))

{'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'}

Automated agreement, including row deletion

The runner compares all statistical fields using abs(actual-reference) <= 1e-9 + 1e-7*abs(reference) and checks retained sample counts and row positions separately. It covers simple, multiple, correlated full-rank, and no-intercept cases. The missing-data case drops one predictor-missing row and one outcome-missing row: both implementations must use the same four observations. No warning text is silently lost: PHP warnings are captured and shown below.

R’s explicit na.omit and statsmodels’ explicit missing='drop' select complete cases. They do not imply that all libraries warn on deletion by default. PHP deliberately emits a warning and returns metadata.

results = [simple_result] + [compare_case(case) for case in cases() if case['name'] != 'simulated_housing']
summary = []
for result in results:
    largest_error = max(error for errors in result['max_absolute_errors_vs_R'].values() for error in errors.values())
    summary.append({'case': result['name'], 'n input': result['php-stats']['n_input'],
                    'n retained': result['php-stats']['n_observations'],
                    'n dropped': result['php-stats']['n_dropped'], 'largest absolute difference vs R': largest_error})
display(pd.DataFrame(summary))
missing = next(result for result in results if result['name'] == 'listwise_missing')
print('Dropped original positions:', missing['php-stats']['dropped_rows'])
print('Retained original positions:', missing['php-stats']['retained_rows'])
print('PHP warning:', missing['warnings'][0]['message'])
print('All three implementations agree for every case and checked field.')

Dropped original positions: [1, 3]
Retained original positions: [0, 2, 4, 5]
PHP warning: stats_ols(): stats_ols dropped 2 of 6 observations with missing values (listwise deletion)
All three implementations agree for every case and checked field.

Read the limits of this evidence

The harness deliberately excludes perfect and zero fits from three-way finite comparisons: floating-point residuals can differ dramatically relative to their tiny magnitudes, and some inference is undefined. The Rust/PHP test suites cover those conventions separately. We also deliberately reject singular designs instead of asking each reference implementation to make its own aliasing decisions.

The tolerance includes an absolute floor; it does not promise relative accuracy for extremely small p-values. This notebook verifies these deterministic cases in the recorded environment, not every possible dataset or platform. See verification policy for coverage and fixture maintenance.

Sources: R lm, statsmodels OLS, statsmodels missing-data behavior.