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. Further model families each get their own module with the diagnostics that are actually well-defined for them:

• Regularized regression ([regularized]) — [RidgeFit], [LassoFit], [ElasticNetFit] and [PenalizedLogisticFit], 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.
• Other GLM families ([glm]) — [GlmFit] over a [Family] trait (Poisson, negative binomial, Gamma) behind one IRLS solver, with deviance / Pearson residuals, the dispersion estimate, McFadden pseudo-R², AIC/BIC, and GLM leverage / Cook's distance.
• Categorical responses ([categorical]) — [MultinomialFit] (baseline-category logit) and [OrdinalFit] (proportional odds), each reducing to binary logistic at K = 2.
• Survival / time-to-event ([survival]) — [CoxFit] proportional hazards (Efron/Breslow ties, plus stratified and time-varying forms), parametric [AftFit] models, and the [KaplanMeier] estimator, with martingale, deviance and Schoenfeld residuals and the concordance index.
• Mixed / hierarchical models ([mixed]) — random-intercept [LinearMixedModel], the general [MixedModel] (random slopes, crossed / nested factors), and generalized [GlmmFit] (Laplace) models, with variance components, intraclass correlation, and BLUPs.

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); elastic net at α = 1 is the lasso and at λ = 0 is OLS, and penalized logistic at λ = 0 is ordinary logistic (see tests/penalized.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).

Other GLM families (Poisson / negative binomial / Gamma)

One IRLS solver serves every family through the [Family] trait; each family is just its link, variance function, deviance and log-likelihood. Anchored on exact identities: the Poisson MLE satisfies Xᵀ(y − μ) = 0, a group-indicator design reproduces the group means on the log scale, deviance residuals square to the residual deviance, GLM leverage sums to p, and the negative binomial collapses onto Poisson as θ → ∞ (see tests/glm.rs).

The dispersion is fixed at 1 for Poisson and the negative binomial and estimated (Pearson χ²/(n − p)) for Gamma; that choice drives the Wald reference distribution (normal vs Student's t) and whether φ̂ is charged as a parameter in AIC/BIC, following R's glm(). The negative-binomial θ can be supplied fixed to NegativeBinomial, or estimated jointly with the coefficients by glm::fit_negative_binomial (a profile-likelihood search over θ).

Categorical responses (multinomial / ordinal)

Anchored on the exact reduction to binary logistic at K = 2 (both models reproduce the LogisticFit coefficients and log-likelihood), plus the softmax / cumulative score equations and the deviance identity (see tests/categorical.rs).

Survival / time-to-event

Validated against the published Freireich 6-MP coxph result — Efron β ≈ 1.572, Breslow β ≈ 1.509 — with the score-zero, residual-sum and hand-computed Kaplan–Meier anchors; the exponential AFT reproduces the closed-form MLE β₀ = ln(Σt/d), and the stratified and counting-process forms reduce exactly to the ordinary fit in the degenerate case (see tests/survival.rs).

Mixed / hierarchical models

Anchored on the exact result that, for a balanced one-way design, the REML variance components equal the ANOVA estimators (σ̂²_e = MSE, σ̂²_b = (MSB − MSE)/m); the GLS intercept is the grand mean, the fit collapses to OLS with no group signal, the general MixedModel reduces to the closed-form LinearMixedModel for a single intercept term, and the Poisson GLMM approaches the plain GLM as σ_b → 0 (see tests/mixed.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.
• Elastic net minimizes (1/2n)‖y−Xβ‖² + λ[α‖β‖₁ + ½(1−α)‖β‖²] on internally-standardized predictors (so α = 1 is exactly the lasso); its effective df is the ridge-shrunk active-set trace, between |active set| and the full ridge df.
• Penalized logistic shrinks a GLM: its effective df is tr[(XᵀWX+λP)⁻¹XᵀWX], falling from p toward the intercept as λ grows, and its standard errors come from the sandwich covariance, not the naïve inverse information.

Limitations

Not blockers, but the boundaries of what the current implementations target:

• Scale. The models that need a general optimizer (parametric survival, the general/generalized mixed models) use dense O(n³) linear algebra — sized for the grouped, diagnostic-scale datasets these tools are for, not for very large n. The OLS/ridge/lasso/logistic/GLM paths remain lightweight.
• GLMM accuracy. [GlmmFit] uses a first-order Laplace approximation (not higher-order adaptive Gaussian quadrature), and covers random intercepts for Poisson/Bernoulli responses.
• Survival extras. Interval-censoring, competing risks, and frailty models are beyond the Cox (standard/stratified/time-varying), parametric-AFT and Kaplan–Meier estimators provided.

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

cargo run --example glm_diagnostics             # Poisson + negative-binomial, over-dispersion, AIC

cargo run --example categorical_diagnostics     # multinomial vs ordinal (proportional-odds) logit

cargo run --example elastic_net_diagnostics     # elastic net + ridge-penalized logistic, effective df

cargo run --example survival_diagnostics        # Kaplan-Meier + Cox PH, Schoenfeld PH check

cargo run --example mixed_model                 # random-intercept LMM, variance components, BLUPs

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.