1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
//! Diagnostics for **generalized linear models** beyond binary logistic — the
//! exponential-dispersion families that share one iteratively-reweighted
//! least-squares (IRLS) engine but differ in link, variance, and residual scale.
//!
//! Where [`logistic`](crate::logistic) is a single hand-written model, this
//! module factors the GLM machinery into a [`Family`] trait and one generic
//! solver [`GlmFit`], so each family is just its link / variance / deviance
//! formulas. Three families are provided:
//!
//! * [`Poisson`] — counts, `V(μ) = μ`, fixed dispersion.
//! * [`NegativeBinomial`] — over-dispersed counts, `V(μ) = μ + μ²/θ` for a fixed
//! `θ`.
//! * [`Gamma`] — positive skewed continuous responses, `V(μ) = μ²`, with an
//! **estimated** dispersion.
//!
//! All three use the **log link** (see [`family`] for why), and each carries the
//! same diagnostic surface:
//!
//! * [`GlmFit`] — the maximum-likelihood fit via IRLS, with coefficient standard
//! errors, Wald statistics, and p-values (normal or Student's *t* depending on
//! whether the dispersion is estimated).
//! * [`deviance_residuals`] / [`pearson_residuals`] — the two GLM residual
//! scales; deviance residuals square to the residual deviance, Pearson
//! residuals to the χ² that defines the dispersion estimate.
//! * [`GoodnessOfFit`] via [`GlmFit::goodness_of_fit`] — null/residual deviance,
//! dispersion, McFadden's pseudo-R², and AIC/BIC.
//! * [`leverage`] and [`cooks_distance`] — the weighted-hat-matrix influence
//! measures.
//!
//! # Response conventions
//!
//! As with OLS and logistic, the **caller owns the design matrix**, intercept
//! column included. The response must lie in the family's support: non-negative
//! counts for [`Poisson`] and [`NegativeBinomial`], strictly positive reals for
//! [`Gamma`]. Out-of-support values are rejected with
//! [`RegressionError::InvalidResponse`](crate::RegressionError::InvalidResponse).
//!
//! # Relationship to the `logistic` module
//!
//! Binary logistic regression is itself a GLM (binomial family, logit link) and
//! could be expressed here, but it keeps its own module: its separation
//! diagnostics and the Hosmer–Lemeshow test are specific to the binary case, and
//! rewriting a shipped, tested API in terms of this trait would buy nothing. This
//! module covers the families logistic does *not*.
//!
//! # What is still out of scope
//!
//! **Multinomial and ordinal** logistic have a vector-valued linear predictor
//! and block covariance — a different solver shape, not another [`Family`] — and
//! live in [`categorical`](crate::categorical), not here. The negative-binomial
//! `θ` can be supplied fixed to [`NegativeBinomial`] or estimated jointly with
//! [`fit_negative_binomial`] (a profile-likelihood search over `θ`).
//!
//! # Quick start
//!
//! ```
//! use ndarray::array;
//! use regression_diagnostics::glm::{GlmFit, Poisson, deviance_residuals};
//!
//! // Caller supplies the intercept column (first column of ones).
//! let x = array![
//! [1.0, 0.0],
//! [1.0, 1.0],
//! [1.0, 2.0],
//! [1.0, 3.0],
//! [1.0, 4.0],
//! [1.0, 5.0],
//! ];
//! let y = array![1.0, 2.0, 3.0, 5.0, 8.0, 13.0]; // roughly exponential growth
//! let fit = GlmFit::new(Poisson, x, y).unwrap();
//!
//! let gof = fit.goodness_of_fit();
//! assert!(gof.residual_deviance <= gof.null_deviance + 1e-9);
//! let dr = deviance_residuals(&fit);
//! let sum_sq: f64 = dr.iter().map(|d| d * d).sum();
//! assert!((sum_sq - gof.residual_deviance).abs() < 1e-8);
//! ```
pub use ;
pub use GlmFit;
pub use fit_negative_binomial;
pub use GoodnessOfFit;
pub use ;
pub use ;