prlmpy (Pseudo-observation Regularized LM)

An exact Python port of the R package prlm.

Fit OLS while adding pseudo-observations to regularize every coefficient toward zero, scaling so the implied prior is roughly β ~ N(0, 1) on each coefficient. It does this by adding a scaled identity matrix to the X variables and a vector of 0s to y.

See Bayesian Data Analysis Third Edition (Gelman et al.), Section 14.8 — "Including Numerical Prior Information" — for the method.

The insight: you can add a single pseudo-observation to regularize a coefficient. This is weak in the sense that a single observation is not much — and if a coefficient doesn't survive adding a single observation, you should be skeptical of it.

This package adds a bit of regularization to OLS in a no-nonsense way; for more complex models reach for numpyro, STAN, brms, etc.

Method

Fit plain OLS to get the regression standard error σ = sqrt(SSR / df_residual), then solve the augmented least-squares system

X_aug = [ X ;  σ·I_p ]      y_aug = [ y ;  0_p ]

This is exactly ridge regression with λ = σ² applied to every coefficient (including the intercept — the identity is the full p × p), which corresponds to a β ~ N(0, 1) prior. σ is computed from the unweighted OLS fit even when weights are supplied (faithful to the R source); weights apply only to the augmented refit, with the pseudo-rows getting weight 1.

Engine

Uses statsmodels (the Python equivalent of R's lm) for the augmented fit, so the returned object exposes the full inference: params, bse, tvalues, pvalues, rsquared, conf_int(), summary(), predict(), … Inference is computed on the augmented system, exactly as R's prlm returns an lm fit on the augmented data.

Install

pip install -e .            # from this folder

Example

import pandas as pd
from prlmpy import prlm, prlm_fit

# Formula interface — mirrors R's prlm(y ~ x, df)
fit = prlm("y ~ x1 + x2", df)            # coefs regularized to ~ N(0,1)
print(fit.summary())
print(fit.params, fit.pvalues, fit.sigma, fit.lam)

# Matrix interface (handy for residual-on-features regressions)
fit = prlm_fit(X, y, weights=w)          # adds an intercept by default
beta = fit.params
resid = fit.resid_original               # residuals on the original n rows only

prlm(formula, data, weights=None) mirrors R; prlm_fit(X, y, weights=None, names=None, add_intercept=True) is the array interface. Both return a PrlmResults that delegates to the underlying statsmodels result and adds sigma, lam (= σ²), fitted_original, and resid_original.