regression-diagnostics 0.2.0

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

/// Overall GLM goodness-of-fit statistics.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct GoodnessOfFit {
    /// Deviance of the intercept-only (mean-only) model.
    pub null_deviance: f64,
    /// Residual deviance of the fitted model (the sum of unit deviances, and of
    /// squared deviance residuals).
    pub residual_deviance: f64,
    /// Degrees of freedom of the null deviance (`n − 1`).
    pub df_null: f64,
    /// Degrees of freedom of the residual deviance (`n − p`).
    pub df_residual: f64,
    /// Estimated dispersion `φ` (`1` for Poisson / negative binomial).
    pub dispersion: f64,
    /// McFadden's pseudo-R², `1 − ℓ/ℓ₀`.
    pub mcfadden_r2: f64,
    /// Akaike information criterion, `−2ℓ + 2k`.
    pub aic: f64,
    /// Bayesian information criterion, `−2ℓ + ln(n)·k`.
    pub bic: f64,
}

impl<F: Family> GlmFit<F> {
    /// Overall goodness-of-fit summary: null/residual deviance, dispersion,
    /// McFadden's pseudo-R², and AIC/BIC.
    ///
    /// The **residual deviance** is the sum of per-observation unit deviances
    /// (equivalently the sum of squared deviance residuals). The **null model**
    /// is the mean-only fit, whose MLE mean is the sample mean `ȳ` for all
    /// log-link families here, so the null deviance is computed directly without
    /// a second IRLS solve.
    ///
    /// The information criteria count `k = p` parameters when the dispersion is
    /// fixed (Poisson, negative binomial) and `k = p + 1` when it is estimated
    /// (Gamma), charging the extra degree of freedom for `φ̂`. McFadden's
    /// pseudo-R² compares the fitted log-likelihood to the mean-only model at the
    /// same dispersion.
    pub fn goodness_of_fit(&self) -> GoodnessOfFit {
        let n = self.n_observations();
        let p = self.n_parameters();
        let y = self.response();
        let mu = self.fitted_means();
        let family = self.family();
        let dispersion = self.dispersion();

        let residual_deviance: f64 = (0..n).map(|i| family.unit_deviance(y[i], mu[i])).sum();

        // Null (mean-only) model: μ ≡ ȳ for every log-link family here.
        let ybar = y.sum() / n as f64;
        let null_deviance: f64 = (0..n).map(|i| family.unit_deviance(y[i], ybar)).sum();

        let ll = self.log_likelihood();
        let ll_null: f64 = (0..n).map(|i| family.loglik(y[i], ybar, dispersion)).sum();
        let mcfadden_r2 = if ll_null != 0.0 {
            1.0 - ll / ll_null
        } else {
            f64::NAN
        };

        // Estimated dispersion counts as one extra parameter in the criteria.
        let k = p as f64 + if family.dispersion_known() { 0.0 } else { 1.0 };
        let aic = -2.0 * ll + 2.0 * k;
        let bic = -2.0 * ll + (n as f64).ln() * k;

        GoodnessOfFit {
            null_deviance,
            residual_deviance,
            df_null: (n - 1) as f64,
            df_residual: (n - p) as f64,
            dispersion,
            mcfadden_r2,
            aic,
            bic,
        }
    }
}