regression-diagnostics 0.2.0

//! Residual diagnostics: the scaled residual forms plus the classic
//! assumption-checking tests.
//!
//! Each test targets a *different* OLS assumption:
//!
//! * [`durbin_watson`] — residual **autocorrelation** (independence).
//! * [`breusch_pagan`] / [`white_test`] — **heteroskedasticity** (constant
//!   variance).
//! * [`jarque_bera`] — residual **normality**.
//!
//! The scaled residual forms ([`standardized_residuals`],
//! [`internally_studentized_residuals`], [`externally_studentized_residuals`])
//! are defined here once and reused by the influence diagnostics
//! ([`crate::influence`]) and QQ-plot data ([`qq_plot_data`]) rather than each
//! module recomputing its own.

mod breusch_pagan;
mod durbin_watson;
mod jarque_bera;
mod qq_plot_data;
mod standardized;
mod studentized;
mod white_test;

pub use breusch_pagan::{breusch_pagan, BreuschPagan};
pub use durbin_watson::durbin_watson;
pub use jarque_bera::{jarque_bera, JarqueBera};
pub use qq_plot_data::qq_plot_data;
pub use standardized::standardized_residuals;
pub use studentized::{externally_studentized_residuals, internally_studentized_residuals};
pub use white_test::{white_test, WhiteTest};

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

use crate::linalg::{aux_r_squared, dmatrix_from_rows, dvector_from_slice};
use crate::OlsFit;

/// Indices of the design columns that are *not* the intercept — the regressors
/// the heteroskedasticity auxiliary regressions are built from.
pub(crate) fn non_intercept_columns(fit: &OlsFit) -> Vec<usize> {
    (0..fit.n_parameters())
        .filter(|&j| fit.intercept_column() != Some(j))
        .collect()
}

/// Lagrange-multiplier heteroskedasticity test shared by Breusch-Pagan and
/// White: regress `target` (the squared residuals) on an intercept plus the
/// supplied `regressors`, then form `LM = n·R²_aux`, which is asymptotically
/// χ²(q) with `q = regressors.len()`.
///
/// Each entry of `regressors` is one auxiliary regressor column of length `n`.
/// Returns `(lm_statistic, df, p_value)`.
pub(crate) fn lm_heteroskedasticity_test(
    n: usize,
    regressors: &[Vec<f64>],
    target: &[f64],
) -> (f64, usize, f64) {
    let q = regressors.len();
    // The auxiliary regression needs positive residual degrees of freedom:
    // `n > q + 1` (regressors plus intercept). With too few observations —
    // e.g. White's test, whose regressor count grows quadratically — the test
    // is not estimable, so report NaN rather than a fabricated statistic.
    if q == 0 || q + 1 >= n {
        return (f64::NAN, q, f64::NAN);
    }
    // Build the [1 | regressors] design in row-major order.
    let mut data = Vec::with_capacity(n * (q + 1));
    for i in 0..n {
        data.push(1.0);
        for reg in regressors {
            data.push(reg[i]);
        }
    }
    let design = dmatrix_from_rows(n, q + 1, &data);
    let target_dv = dvector_from_slice(target);

    let r2 = aux_r_squared(&design, &target_dv).unwrap_or(0.0);
    let lm = n as f64 * r2;
    let p_value = match ChiSquared::new(q as f64) {
        Ok(dist) if lm.is_finite() && lm >= 0.0 && q > 0 => 1.0 - dist.cdf(lm),
        _ => f64::NAN,
    };
    (lm, q, p_value)
}