regression-diagnostics

Statistical regression diagnostics for Rust — the kind of surface R's car/lmtest and Python's statsmodels expose, packaged as a dependency-light crate. It fills a specific ecosystem gap: no maintained Rust crate offers Variance Inflation Factor, adjusted R², or residual diagnostics (autocorrelation, heteroskedasticity, normality, influence).

The core is ordinary least squares, built on one fitted-model type, [OlsFit], computed once via a single QR factorization and reused by every diagnostic. Two further model families each get their own module with the diagnostics that are actually well-defined for them:

• Regularized regression ([regularized]) — [RidgeFit] and [LassoFit], with effective degrees of freedom, ridge leverage/GCV, and a ridge-VIF that reduces to the OLS VIF at λ = 0 (the "VIF before/after regularization" comparison).
• Binary logistic regression ([logistic]) — [LogisticFit] via IRLS, with deviance/Pearson residuals, McFadden pseudo-R², the Hosmer–Lemeshow test, and logistic leverage / Cook's distance.

use ndarray::array;
use regression_diagnostics::OlsFit;

// Caller supplies the intercept column (leading ones).
let x = array![
    [1.0, 1.0, 2.0],
    [1.0, 2.0, 4.1],
    [1.0, 3.0, 5.9],
    [1.0, 4.0, 8.0],
    [1.0, 5.0, 10.1],
];
let y = array![2.0, 4.1, 6.1, 8.0, 10.2];
let fit = OlsFit::new(x, y).unwrap();

println!("{}", fit.summary());   // full statsmodels-style report

fit.summary() renders like this (from cargo run --example full_diagnostics):

============================== OLS Diagnostics ===============================
No. Observations:     12    Df Residuals:      8    Df Model:      3
R-squared:          0.9797  Adj. R-squared:   0.9720  Resid. SE:  17.1435
F-statistic:      128.4729  Prob(F):          0.0000  Log-Lik:     -48.69
AIC:                105.39  BIC:              107.33  Cond. No.:  2.012e2
------------------------------------------------------------------------------
                coef    std err        t    P>|t|     beta      VIF
------------------------------------------------------------------------------
const       133.3251    28.5980    4.662    0.002     -        -
x1           17.7810     2.7081    6.566    0.000    1.471    19.74
x2           -6.8850     0.5978  -11.517    0.000   -0.793     1.86
x3          -23.1377    19.6447   -1.178    0.273   -0.242    16.63
------------------------------------------------------------------------------
Durbin-Watson:      2.2030   (residual autocorrelation; ~2 is ideal)
Jarque-Bera:        0.4109   Prob:  0.8143   (skew 0.217, kurt 2.204)
Breusch-Pagan:      1.2231   Prob:  0.7475   (heteroskedasticity, LM)
White:              7.1759   Prob:  0.6188   (heteroskedasticity, general)
==============================================================================

Summary is structured data first — every number above is a field you can read programmatically; the table is just a Display convenience on top.

What's covered

OLS

Validated against statsmodels' published Longley-dataset output — coefficients, R²/adjusted R², F, log-likelihood, AIC, BIC, Durbin-Watson, and the six per-predictor VIFs all match (see tests/reference_datasets/longley.rs).

Regularized regression

Anchored on exact identities: RidgeFit/LassoFit at λ = 0 reproduce the OLS coefficients, ridge leverage sums to the effective df, and ridge VIF equals the OLS VIF at λ = 0 (see tests/regularized.rs).

Binary logistic regression

Anchored on exact identities too: the fit satisfies the score equations Xᵀ(y − p) = 0, a single binary predictor recovers the log-odds and log-odds-ratio in closed form, and deviance residuals square to the residual deviance (see tests/logistic.rs).

Linear algebra dependency

The internals use nalgebra (pure-Rust QR/SVD) while the public API speaks ndarray. That split is deliberate:

• ndarray-linalg would be faster but pulls in a system-level BLAS/LAPACK, which cuts against the portability goal (this series is Docker-first and multi-arch).
• nalgebra is pure Rust — no external BLAS — so the crate stays portable, at the cost of a small conversion layer at the API boundary, since the rest of the guide's code is ndarray-based.

We take that trade in favour of portability: nalgebra is used purely internally and never appears in the public API, so callers only ever see ndarray types.

Two numerical choices worth flagging:

• Coefficients come from a QR solve of X, never from inverting XᵀX. Forming XᵀX squares the condition number — exactly the quantity the multicollinearity diagnostics measure — which would make VIF and the condition number untrustworthy on the ill-conditioned designs where they matter most.
• Leverage (diag(H)) is read from the thin Q factor as row norms, so the full n × n hat matrix is never materialized (O(n·p) memory, not O(n²)).

Conventions worth knowing

• Intercept: the caller owns the design matrix. OlsFit::new fits the columns as given (auto-detecting a constant column as the intercept); OlsFit::with_intercept prepends a ones column for you. Getting this wrong silently biases every downstream statistic, so it's explicit at the call site.
• AIC/BIC count the regression parameters (intercept included) as k, not σ² — matching statsmodels' convention so values line up.
• Condition number is computed on the design as fitted (unscaled columns). Tools that normalize columns first report a different number; the convention is documented so the value is reproducible.
• Thresholds (VIF > 5/10, condition number > 30, Cook's D > 4/n, …) are documented as convention, not mathematical fact — they vary by source.
• Undefined-per-column results (VIF and beta of the intercept) are returned as NaN in the column's slot, aligned to the design columns, rather than shifting indices.

Regularization and logistic conventions

Because ridge, lasso and logistic have different structure from OLS, a few choices are worth stating (each is documented on the type):

• Ridge/lasso intercept is an unpenalized, mean-centered constant (the scikit-learn convention). Lasso standardizes predictors internally, so its λ (on the (1/2n)‖y−Xβ‖² + λ‖β‖₁ objective) is not comparable to ridge's. Neither estimator is scale-invariant — standardize predictors if that matters.
• Effective degrees of freedom replace the raw parameter count under shrinkage: Σ dⱼ²/(dⱼ²+λ) for ridge, the active-set size for lasso. The information criteria use these, so a shrunk model is charged for the freedom it actually spends.
• Ridge VIF is the diagonal of (R+λI)⁻¹R(R+λI)⁻¹ on standardized predictors; it can fall below 1, because shrinkage reduces coefficient variance below the orthogonal-OLS baseline — that is the point of the before/after comparison.
• Logistic requires a 0/1 response with both classes present. Perfect or quasi-complete separation makes the MLE diverge; IRLS then returns NotConverged rather than reporting meaningless coefficients.

Out of scope for v0.1

Still noted rather than silently omitted:

• Other GLM families — Poisson/negative-binomial counts, multinomial/ordinal logistic, gamma. Each needs its own link/variance functions and residual scales; only binary logistic is implemented.
• Penalized-GLM and elastic-net diagnostics, survival models, and mixed / hierarchical models — separate frameworks, deferred to possible future crates.

Examples

cargo run --example full_diagnostics
cargo run --example multicollinearity_check
cargo run --example heteroskedasticity_check
cargo run --example qq_plot_integration
cargo run --example ridge_diagnostics       # VIF before/after regularization, GCV
cargo run --example logistic_diagnostics    # IRLS fit, deviance, Hosmer-Lemeshow

qq_plot_integration shows the intended cross-crate pairing: this crate produces QQ-plot data, and a plotting crate (plotters / plotters-statistical) renders it against a y = x normality reference line.

License

MIT.