regression_diagnostics/logistic/
influence.rs

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

regression_diagnostics/logistic/
influence.rs
