regression-diagnostics 0.2.0

use ndarray::Array1;

use super::{pearson_residuals, LogisticFit};

/// Logistic 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
/// `wᵢ = pᵢ(1 − pᵢ)`.
///
/// The GLM analogue of OLS leverage: it measures how much observation `i`'s own
/// fitted value is determined by its predictors, but weighted by the binomial
/// variance, so points where the model is already near-certain (`pᵢ ≈ 0` or `1`)
/// carry little leverage. Computed from the stored covariance without forming the
/// `n × n` hat matrix, and the values sum to `p`.
pub fn leverage(fit: &LogisticFit) -> Array1<f64> {
    let x = fit.design_matrix();
    let cov = fit.covariance();
    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 logistic regression (Pregibon):
///
/// `Cᵢ = r_pᵢ² · hᵢ / (p · (1 − hᵢ)²)`,
///
/// where `r_pᵢ` is the Pearson residual and `hᵢ` the logistic leverage. Like its
/// OLS counterpart it combines residual size and leverage into a single
/// per-observation influence measure — large when a point is both poorly fit and
/// has unusual, well-weighted predictor values — and flags observations whose
/// removal would most move the coefficients.
pub fn cooks_distance(fit: &LogisticFit) -> 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)
    })
}