regression-diagnostics 0.2.0

use ndarray::Array1;

use super::family::Family;
use super::GlmFit;

/// Pearson residuals `(yᵢ − μᵢ) / √V(μᵢ)`.
///
/// The raw residual scaled by the family's standard deviation (the `√φ` factor
/// is left out, matching R). Their sum of squares is the Pearson χ² statistic,
/// which — divided by `n − p` — is exactly the dispersion estimate the Gamma fit
/// reports. Pearson residuals are simple but skewed for the count families; for
/// a more symmetric scale prefer [`deviance_residuals`].
pub fn pearson_residuals<F: Family>(fit: &GlmFit<F>) -> Array1<f64> {
    let y = fit.response();
    let mu = fit.fitted_means();
    let family = fit.family();
    Array1::from_shape_fn(fit.n_observations(), |i| {
        let denom = family.variance(mu[i]).sqrt();
        if denom > 0.0 {
            (y[i] - mu[i]) / denom
        } else {
            f64::NAN
        }
    })
}

/// Deviance residuals `sign(yᵢ − μᵢ) · √d(yᵢ, μᵢ)`.
///
/// The signed square root of each observation's unit 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 GLM,
/// which is why they are the default for GLM residual plots.
pub fn deviance_residuals<F: Family>(fit: &GlmFit<F>) -> Array1<f64> {
    let y = fit.response();
    let mu = fit.fitted_means();
    let family = fit.family();
    Array1::from_shape_fn(fit.n_observations(), |i| {
        let dev = family.unit_deviance(y[i], mu[i]);
        let sign = if y[i] - mu[i] >= 0.0 { 1.0 } else { -1.0 };
        sign * dev.max(0.0).sqrt()
    })
}