regression_diagnostics/glm/
influence.rs

1use ndarray::Array1;
2
3use super::family::Family;
4use super::{pearson_residuals, GlmFit};
5
6/// GLM leverage — the diagonal of the weighted hat matrix
7/// `H = W^{1/2}X(XᵀWX)⁻¹XᵀW^{1/2}`, i.e. `hᵢ = wᵢ · xᵢᵀ(XᵀWX)⁻¹xᵢ` with the IRLS
8/// working weight `wᵢ = (dμ/dη)² / V(μᵢ)`.
9///
10/// The GLM analogue of OLS leverage: how much observation `i`'s own fitted value
11/// is determined by its predictors, weighted by the family variance. The
12/// dispersion cancels (it scales `W` and `(XᵀWX)⁻¹` inversely), so leverage is a
13/// pure geometric quantity; the values sum to `p`. Computed from the stored
14/// inverse information without forming the `n × n` hat matrix.
15pub fn leverage<F: Family>(fit: &GlmFit<F>) -> Array1<f64> {
16    let x = fit.design_matrix();
17    let cov = fit.cov_unscaled();
18    let w = fit.weights();
19    let p = fit.n_parameters();
20
21    Array1::from_shape_fn(fit.n_observations(), |i| {
22        // quad = xᵢᵀ (XᵀWX)⁻¹ xᵢ
23        let mut quad = 0.0;
24        for a in 0..p {
25            let mut inner = 0.0;
26            for b in 0..p {
27                inner += cov[(a, b)] * x[(i, b)];
28            }
29            quad += x[(i, a)] * inner;
30        }
31        w[i] * quad
32    })
33}
34
35/// Cook's-distance analogue for a GLM (Pregibon):
36///
37/// `Cᵢ = r_pᵢ² · hᵢ / (p · (1 − hᵢ)²)`,
38///
39/// where `r_pᵢ` is the Pearson residual and `hᵢ` the GLM leverage. As in OLS it
40/// combines residual size and leverage into one per-observation influence
41/// measure — large when a point is both poorly fit and has unusual,
42/// well-weighted predictor values — flagging observations whose removal would
43/// most move the coefficients.
44pub fn cooks_distance<F: Family>(fit: &GlmFit<F>) -> Array1<f64> {
45    let h = leverage(fit);
46    let rp = pearson_residuals(fit);
47    let p = fit.n_parameters() as f64;
48
49    Array1::from_shape_fn(fit.n_observations(), |i| {
50        let one_minus_h = 1.0 - h[i];
51        if one_minus_h <= 0.0 || p <= 0.0 {
52            return f64::NAN;
53        }
54        rp[i] * rp[i] * h[i] / (p * one_minus_h * one_minus_h)
55    })
56}
regression_diagnostics/glm/influence.rs

regression_diagnostics/glm/
influence.rs
