# 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.
```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}` |
| Elastic-net fit (L1+L2, coordinate descent) | `regularized::ElasticNetFit` |
| Shrinkage-aware effective df / AIC | `ElasticNetFit::{effective_df, aic}` |
| Ridge-penalized logistic (penalized IRLS) | `regularized::PenalizedLogisticFit` |
| Penalized effective df `tr[(XᵀWX+λP)⁻¹XᵀWX]`, sandwich SEs | `PenalizedLogisticFit::{effective_df, coefficient_standard_errors}` |
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
| 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`).
### Other GLM families (Poisson / negative binomial / Gamma)
| GLM fit (IRLS / Fisher scoring) | `glm::GlmFit::new(family, x, y)` | `GLM(y, X, family=...)` |
| Families (all log-link) | `glm::{Poisson, NegativeBinomial, Gamma}` | `families.{Poisson, NegativeBinomial, Gamma}` |
| Coefficients, std errors, Wald stat, p | `GlmFit::{coefficients, coefficient_standard_errors, wald_statistics, p_values}` | `.params`, `.bse`, `.tvalues`, `.pvalues` |
| Dispersion `φ` (fixed vs estimated) | `GlmFit::dispersion` | `.scale` |
| Deviance / Pearson residuals | `glm::{deviance_residuals, pearson_residuals}` | `.resid_deviance`, `.resid_pearson` |
| Null/residual deviance, McFadden R², AIC/BIC | `GlmFit::goodness_of_fit` | `.null_deviance`, `.deviance`, `.aic` |
| GLM leverage / Cook's distance | `glm::{leverage, cooks_distance}` | `GLMInfluence.hat_matrix_diag`, `.cooks_distance` |
| Predicted means for new `X` | `GlmFit::predict_mean` | `.predict` |
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)
| Multinomial (baseline-category) logit | `categorical::MultinomialFit` | `nnet::multinom` / `MNLogit` |
| Ordinal (proportional-odds) logit | `categorical::OrdinalFit` | `MASS::polr` / `OrderedModel` |
| Coefficients, std errors, z, p | `{coefficients, coefficient_standard_errors, z_values, p_values}` | `.params`, `.bse`, … |
| Deviance, McFadden R², AIC/BIC, predicted probs | `{residual_deviance, mcfadden_r2, aic, bic, predict_proba}` | `.deviance`, `.aic`, `.predict` |
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
| Cox proportional hazards (Efron/Breslow ties) | `survival::CoxFit`, `survival::Ties` | `coxph` |
| Coefficients (log HR), HR, std errors, z, p | `CoxFit::{coefficients, hazard_ratios, coefficient_standard_errors, z_values, p_values}` | `coef`, `exp(coef)`, … |
| Stratified Cox (baseline per stratum) | `CoxFit::stratified` | `coxph(... + strata())` |
| Time-varying covariates `(start, stop]` | `CoxFit::counting_process` | `coxph(Surv(start,stop,ev)~…)` |
| Breslow baseline cumulative hazard | `CoxFit::baseline_cumulative_hazard` | `basehaz` |
| Concordance (C-index) | `CoxFit::concordance` | `concordance` |
| Martingale / deviance / Schoenfeld residuals | `survival::{martingale_residuals, deviance_residuals, schoenfeld_residuals}` | `residuals(type=…)` |
| Parametric AFT (Weibull/exp/log-normal/log-logistic) | `survival::AftFit`, `survival::AftDistribution` | `survreg` |
| Kaplan–Meier (Greenwood SE, median) | `survival::KaplanMeier` | `survfit` |
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
| Random-intercept LMM (REML / ML, closed form) | `mixed::LinearMixedModel`, `mixed::Method` | `lmer(y ~ x + (1|g))` |
| Random slopes / crossed / nested LMM | `mixed::MixedModel`, `mixed::RandomEffect` | `lmer(y ~ x + (x|g) + (1|h))` |
| Generalized mixed model (Poisson/Bernoulli, Laplace) | `mixed::GlmmFit`, `mixed::GlmmFamily` | `glmer` |
| Fixed effects + GLS std errors | `{coefficients, coefficient_standard_errors}` | `fixef`, `vcov` |
| Variance components / term covariances | `{group_variance, residual_variance, term_covariance}` | `VarCorr` |
| Intraclass correlation | `LinearMixedModel::icc` | (derived) |
| BLUPs (predicted random effects) | `{random_effects}` | `ranef` |
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`](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.
* **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
```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
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.