# 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.
```rust
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`):
```text
============================== 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
| Coefficients, std errors, t, p | `OlsFit` + `summary()` | `.params`, `.bse`, `.tvalues`, `.pvalues` |
| R² / adjusted R² | `fit_statistics::r_squared` / `adjusted_r_squared` | `.rsquared` / `.rsquared_adj` |
| F-statistic (overall) | `fit_statistics::f_statistic` | `.fvalue`, `.f_pvalue` |
| Log-likelihood, AIC, BIC | `fit_statistics::{log_likelihood, aic, bic}` | `.llf`, `.aic`, `.bic` |
| VIF (per predictor) | `multicollinearity::vif` | `statsmodels.stats.outliers_influence.variance_inflation_factor` |
| Condition number | `multicollinearity::condition_number` | `.condition_number` (see scaling note) |
| Durbin-Watson | `residuals::durbin_watson` | `statsmodels.stats.stattools.durbin_watson` |
| Breusch-Pagan | `residuals::breusch_pagan` | `het_breuschpagan` |
| White's test | `residuals::white_test` | `het_white` |
| Jarque-Bera | `residuals::jarque_bera` | `jarque_bera` |
| Standardized / studentized residuals | `residuals::{standardized, internally_studentized, externally_studentized}_residuals` | `OLSInfluence.resid_studentized_{internal,external}` |
| Leverage | `influence::leverage` | `OLSInfluence.hat_matrix_diag` |
| Cook's distance | `influence::cooks_distance` | `OLSInfluence.cooks_distance` |
| DFFITS | `influence::dffits` | `OLSInfluence.dffits` |
| QQ-plot data | `residuals::qq_plot_data` | `qqplot` (data only; rendering is left to a plotting crate) |
| Standardized (beta) coefficients | `coefficients::standardized_coefficients` | (not built in) |
| One-call report | `OlsFit::summary()` | `.summary()` |
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
| Ridge fit (closed form via SVD) | `regularized::RidgeFit` |
| Effective degrees of freedom `Σ dⱼ²/(dⱼ²+λ)` | `RidgeFit::effective_df` |
| Ridge leverage (`diag H_λ`) | `RidgeFit::leverage` |
| GCV score / GCV λ-selection | `RidgeFit::gcv`, `regularized::select_lambda_gcv` |
| Effective AIC/BIC | `RidgeFit::{aic, bic}` |
| Ridge VIF (→ OLS VIF at λ=0) | `RidgeFit::ridge_vif` |
| Lasso fit (coordinate descent) | `regularized::LassoFit` |
| Active set / non-zero count (df estimate) | `LassoFit::{active_set, n_nonzero}` |
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
| Logistic fit (IRLS) | `logistic::LogisticFit` | `Logit` / `GLM(Binomial)` |
| Coefficients, std errors, z, p | `LogisticFit::{coefficients, coefficient_standard_errors, z_values, p_values}` | `.params`, `.bse`, `.tvalues`, `.pvalues` |
| Deviance / Pearson residuals | `logistic::{deviance_residuals, pearson_residuals}` | `.resid_deviance`, `.resid_pearson` |
| Null/residual deviance, McFadden R², AIC/BIC | `LogisticFit::goodness_of_fit` | `.null_deviance`, `.deviance`, `.prsquared`, `.aic` |
| Hosmer–Lemeshow test | `LogisticFit::hosmer_lemeshow` | (statsmodels contrib) |
| Logistic leverage / Cook's distance | `logistic::{leverage, cooks_distance}` | `GLMInfluence.hat_matrix_diag`, `.cooks_distance` |
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`](https://crates.io/crates/nalgebra) (pure-Rust
QR/SVD) while the **public API speaks [`ndarray`](https://crates.io/crates/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
```bash
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.