regression-diagnostics 0.2.0

use ndarray::Array1;

use super::family::Family;
use super::{pearson_residuals, GlmFit};

/// GLM leverage — the diagonal of the weighted hat matrix
/// `H = W^{1/2}X(XᵀWX)⁻¹XᵀW^{1/2}`, i.e. `hᵢ = wᵢ · xᵢᵀ(XᵀWX)⁻¹xᵢ` with the IRLS
/// working weight `wᵢ = (dμ/dη)² / V(μᵢ)`.
///
/// The GLM analogue of OLS leverage: how much observation `i`'s own fitted value
/// is determined by its predictors, weighted by the family variance. The
/// dispersion cancels (it scales `W` and `(XᵀWX)⁻¹` inversely), so leverage is a
/// pure geometric quantity; the values sum to `p`. Computed from the stored
/// inverse information without forming the `n × n` hat matrix.
pub fn leverage<F: Family>(fit: &GlmFit<F>) -> Array1<f64> {
    let x = fit.design_matrix();
    let cov = fit.cov_unscaled();
    let w = fit.weights();
    let p = fit.n_parameters();

    Array1::from_shape_fn(fit.n_observations(), |i| {
        // quad = xᵢᵀ (XᵀWX)⁻¹ xᵢ
        let mut quad = 0.0;
        for a in 0..p {
            let mut inner = 0.0;
            for b in 0..p {
                inner += cov[(a, b)] * x[(i, b)];
            }
            quad += x[(i, a)] * inner;
        }
        w[i] * quad
    })
}

/// Cook's-distance analogue for a GLM (Pregibon):
///
/// `Cᵢ = r_pᵢ² · hᵢ / (p · (1 − hᵢ)²)`,
///
/// where `r_pᵢ` is the Pearson residual and `hᵢ` the GLM leverage. As in OLS it
/// combines residual size and leverage into one per-observation influence
/// measure — large when a point is both poorly fit and has unusual,
/// well-weighted predictor values — flagging observations whose removal would
/// most move the coefficients.
pub fn cooks_distance<F: Family>(fit: &GlmFit<F>) -> Array1<f64> {
    let h = leverage(fit);
    let rp = pearson_residuals(fit);
    let p = fit.n_parameters() as f64;

    Array1::from_shape_fn(fit.n_observations(), |i| {
        let one_minus_h = 1.0 - h[i];
        if one_minus_h <= 0.0 || p <= 0.0 {
            return f64::NAN;
        }
        rp[i] * rp[i] * h[i] / (p * one_minus_h * one_minus_h)
    })
}