regression-diagnostics 0.2.0

use ndarray::Array1;

use super::LogisticFit;

/// Pearson residuals `(yᵢ − pᵢ) / √(pᵢ(1 − pᵢ))`.
///
/// The raw residual scaled by the binomial standard deviation of that
/// observation. Their sum of squares is the Pearson χ² goodness-of-fit statistic.
/// Pearson residuals are simple but skewed for probabilities near 0 or 1; for a
/// more symmetric scale prefer [`deviance_residuals`].
pub fn pearson_residuals(fit: &LogisticFit) -> Array1<f64> {
    let y = fit.response();
    let p = fit.fitted_probabilities();
    Array1::from_shape_fn(fit.n_observations(), |i| {
        let denom = (p[i] * (1.0 - p[i])).sqrt();
        if denom > 0.0 {
            (y[i] - p[i]) / denom
        } else {
            f64::NAN
        }
    })
}

/// Deviance residuals
/// `sign(yᵢ − pᵢ)·√(−2[yᵢ ln pᵢ + (1 − yᵢ) ln(1 − pᵢ)])`.
///
/// The signed square-root of each observation's contribution to the model
/// deviance; their sum of squares **is** the residual deviance. Deviance
/// residuals are the more symmetric, better-behaved scale for spotting poorly-fit
/// points in a logistic model, which is why they are the default for GLM residual
/// plots.
pub fn deviance_residuals(fit: &LogisticFit) -> Array1<f64> {
    let y = fit.response();
    let p = fit.fitted_probabilities();
    Array1::from_shape_fn(fit.n_observations(), |i| {
        let yi = y[i];
        let pi = p[i];
        // Per-observation deviance contribution (always >= 0).
        let dev = -2.0 * (yi * pi.ln() + (1.0 - yi) * (1.0 - pi).ln());
        let sign = if yi - pi >= 0.0 { 1.0 } else { -1.0 };
        sign * dev.max(0.0).sqrt()
    })
}