regression-diagnostics 0.2.0

//! [`Summary`] — a one-call, R/`statsmodels`-style diagnostic report.

use std::fmt;

use statrs::distribution::{ContinuousCDF, StudentsT};

use crate::coefficients::standardized_coefficients;
use crate::fit_statistics::{
    adjusted_r_squared, aic, bic, f_statistic, log_likelihood, r_squared, FStatistic,
};
use crate::influence::cooks_distance;
use crate::multicollinearity::{condition_number, vif};
use crate::residuals::white_test;
use crate::residuals::{
    breusch_pagan, durbin_watson, jarque_bera, BreuschPagan, JarqueBera, WhiteTest,
};
use crate::OlsFit;

/// One row of the coefficient table.
#[derive(Debug, Clone, PartialEq)]
pub struct CoefficientRow {
    /// Display name (`const` for the intercept, otherwise `x1`, `x2`, …).
    pub name: String,
    /// Point estimate `βⱼ`.
    pub estimate: f64,
    /// Standard error `sⱼ`.
    pub std_error: f64,
    /// t statistic `βⱼ / sⱼ`.
    pub t_value: f64,
    /// Two-sided p-value under the t distribution with `n − p` degrees of freedom.
    pub p_value: f64,
    /// Standardized (beta) coefficient (`NaN` for the intercept).
    pub std_coefficient: f64,
    /// Variance Inflation Factor (`NaN` for the intercept).
    pub vif: f64,
}

/// A structured snapshot of every diagnostic in the crate for one fit.
///
/// `Summary` is **structured data first**: every field is a real number you can
/// read programmatically. The [`Display`](fmt::Display) impl is a convenience
/// layer that renders it as a `statsmodels`-style table — it is not the primary
/// interface, so you never have to parse text to get a value back out.
///
/// Build it with [`OlsFit::summary`].
#[derive(Debug, Clone)]
pub struct Summary {
    /// Number of observations.
    pub n_observations: usize,
    /// Number of parameters.
    pub n_parameters: usize,
    /// Residual degrees of freedom.
    pub df_residual: f64,
    /// Whether an intercept is present.
    pub has_intercept: bool,
    /// One row per coefficient.
    pub coefficients: Vec<CoefficientRow>,
    /// `R²`.
    pub r_squared: f64,
    /// Adjusted `R²`.
    pub adj_r_squared: f64,
    /// Overall F-test.
    pub f_statistic: FStatistic,
    /// Residual standard error `s`.
    pub residual_std_error: f64,
    /// Gaussian log-likelihood.
    pub log_likelihood: f64,
    /// Akaike information criterion.
    pub aic: f64,
    /// Bayesian information criterion.
    pub bic: f64,
    /// Design-matrix condition number.
    pub condition_number: f64,
    /// Durbin-Watson statistic.
    pub durbin_watson: f64,
    /// Jarque-Bera normality test.
    pub jarque_bera: JarqueBera,
    /// Breusch-Pagan heteroskedasticity test.
    pub breusch_pagan: BreuschPagan,
    /// White's heteroskedasticity test.
    pub white: WhiteTest,
}

impl OlsFit {
    /// Compute the full [`Summary`] for this fit — every statistic in Milestones
    /// 2–6 in one call.
    ///
    /// ```
    /// use ndarray::array;
    /// use regression_diagnostics::OlsFit;
    ///
    /// let x = array![[1.0, 0.0], [1.0, 1.0], [1.0, 2.0], [1.0, 3.0], [1.0, 4.0]];
    /// let y = array![1.0, 3.1, 4.9, 7.0, 9.1];
    /// let fit = OlsFit::new(x, y).unwrap();
    /// let s = fit.summary();
    /// println!("{s}");
    /// assert!(s.r_squared > 0.99);
    /// ```
    pub fn summary(&self) -> Summary {
        let se = self.coefficient_standard_errors();
        let coef = self.coefficients();
        let vifs = vif(self);
        let std_coefs = standardized_coefficients(self);
        let df = self.df_residual();
        let t_dist = StudentsT::new(0.0, 1.0, df).ok();

        let mut predictor_counter = 0usize;
        let coefficients = (0..self.n_parameters())
            .map(|j| {
                let name = if self.intercept_column() == Some(j) {
                    "const".to_string()
                } else {
                    predictor_counter += 1;
                    format!("x{predictor_counter}")
                };
                let est = coef[j];
                let s = se[j];
                let t = if s > 0.0 { est / s } else { f64::NAN };
                let p = match &t_dist {
                    Some(d) if t.is_finite() => 2.0 * (1.0 - d.cdf(t.abs())),
                    _ => f64::NAN,
                };
                CoefficientRow {
                    name,
                    estimate: est,
                    std_error: s,
                    t_value: t,
                    p_value: p,
                    std_coefficient: std_coefs[j],
                    vif: vifs[j],
                }
            })
            .collect();

        Summary {
            n_observations: self.n_observations(),
            n_parameters: self.n_parameters(),
            df_residual: df,
            has_intercept: self.has_intercept(),
            coefficients,
            r_squared: r_squared(self),
            adj_r_squared: adjusted_r_squared(self),
            f_statistic: f_statistic(self),
            residual_std_error: self.residual_standard_error(),
            log_likelihood: log_likelihood(self),
            aic: aic(self),
            bic: bic(self),
            condition_number: condition_number(self),
            durbin_watson: durbin_watson(self),
            jarque_bera: jarque_bera(self),
            breusch_pagan: breusch_pagan(self),
            white: white_test(self),
        }
    }
}

impl Summary {
    /// Maximum absolute Cook's distance across observations, recomputed from the
    /// fit — a convenience for callers that want a single influence headline
    /// number without walking the full vector. (Not stored on `Summary` because
    /// it is per-observation; use [`crate::influence::cooks_distance`] for the
    /// full vector.)
    pub fn max_cooks_distance(fit: &OlsFit) -> f64 {
        cooks_distance(fit)
            .iter()
            .copied()
            .filter(|v| v.is_finite())
            .fold(0.0_f64, f64::max)
    }
}

fn flag(p: f64, low: f64, high: f64, hi_is_bad: bool) -> &'static str {
    if p.is_nan() {
        return "";
    }
    let bad = if hi_is_bad { p > high } else { p < low };
    if bad {
        " (!)"
    } else {
        ""
    }
}

impl fmt::Display for Summary {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        writeln!(f, "{:=^78}", " OLS Diagnostics ")?;
        writeln!(
            f,
            "No. Observations: {:>6}    Df Residuals: {:>6}    Df Model: {:>6}",
            self.n_observations,
            self.df_residual as usize,
            self.n_parameters - usize::from(self.has_intercept),
        )?;
        writeln!(
            f,
            "R-squared:        {:>8.4}  Adj. R-squared: {:>8.4}  Resid. SE: {:>8.4}",
            self.r_squared, self.adj_r_squared, self.residual_std_error,
        )?;
        writeln!(
            f,
            "F-statistic:      {:>8.4}  Prob(F):        {:>8.4}  Log-Lik:   {:>8.2}",
            self.f_statistic.statistic, self.f_statistic.p_value, self.log_likelihood,
        )?;
        writeln!(
            f,
            "AIC:              {:>8.2}  BIC:            {:>8.2}  Cond. No.: {:>8.3e}",
            self.aic, self.bic, self.condition_number,
        )?;

        writeln!(f, "{:-<78}", "")?;
        writeln!(
            f,
            "{:<8}{:>12}{:>11}{:>9}{:>9}{:>9}{:>9}",
            "", "coef", "std err", "t", "P>|t|", "beta", "VIF",
        )?;
        writeln!(f, "{:-<78}", "")?;
        for row in &self.coefficients {
            let beta = if row.std_coefficient.is_nan() {
                "     -   ".to_string()
            } else {
                format!("{:>9.3}", row.std_coefficient)
            };
            let vif = if row.vif.is_nan() {
                "     -   ".to_string()
            } else if row.vif.is_infinite() {
                "      inf".to_string()
            } else {
                format!("{:>9.2}", row.vif)
            };
            writeln!(
                f,
                "{:<8}{:>12.4}{:>11.4}{:>9.3}{:>9.3}{beta}{vif}",
                row.name, row.estimate, row.std_error, row.t_value, row.p_value,
            )?;
        }
        writeln!(f, "{:-<78}", "")?;

        writeln!(
            f,
            "Durbin-Watson:    {:>8.4}   (residual autocorrelation; ~2 is ideal)",
            self.durbin_watson,
        )?;
        writeln!(
            f,
            "Jarque-Bera:      {:>8.4}   Prob: {:>7.4}{}   (skew {:.3}, kurt {:.3})",
            self.jarque_bera.statistic,
            self.jarque_bera.p_value,
            flag(self.jarque_bera.p_value, 0.05, 0.0, false),
            self.jarque_bera.skewness,
            self.jarque_bera.kurtosis,
        )?;
        writeln!(
            f,
            "Breusch-Pagan:    {:>8.4}   Prob: {:>7.4}{}   (heteroskedasticity, LM)",
            self.breusch_pagan.statistic,
            self.breusch_pagan.p_value,
            flag(self.breusch_pagan.p_value, 0.05, 0.0, false),
        )?;
        writeln!(
            f,
            "White:            {:>8.4}   Prob: {:>7.4}{}   (heteroskedasticity, general)",
            self.white.statistic,
            self.white.p_value,
            flag(self.white.p_value, 0.05, 0.0, false),
        )?;
        write!(f, "{:=<78}", "")?;
        Ok(())
    }
}