regression-diagnostics 0.2.0

//! Constructed-dataset tests for the diagnostics whose correctness is best shown
//! by triggering the exact condition they detect (influence, normality) rather
//! than by matching a published table.

use ndarray::{Array1, Array2};
use regression_diagnostics::coefficients::standardized_coefficients;
use regression_diagnostics::influence::{cooks_distance, dffits};
use regression_diagnostics::residuals::{
    externally_studentized_residuals, internally_studentized_residuals, jarque_bera, qq_plot_data,
    standardized_residuals,
};
use regression_diagnostics::OlsFit;

/// Fit an intercept-only model, so the residuals equal `y − ȳ` exactly — a clean
/// way to drive the residual-distribution diagnostics with a known sample.
fn mean_model(y: &[f64]) -> OlsFit {
    let n = y.len();
    let x = Array2::<f64>::ones((n, 1));
    OlsFit::new(x, Array1::from(y.to_vec())).unwrap()
}

#[test]
fn injected_outlier_is_flagged_by_all_three_influence_measures() {
    // Line y = 2x on x = 1..=20 with a little well-behaved noise (so the
    // leave-one-out fit is not degenerate), then corrupt one interior point hard.
    let n = 20usize;
    let mut x = Array2::<f64>::ones((n, 2));
    let mut y = Array1::<f64>::zeros(n);
    for i in 0..n {
        let xi = (i + 1) as f64;
        x[(i, 1)] = xi;
        let noise = if i % 2 == 0 { 0.4 } else { -0.4 };
        y[i] = 2.0 * xi + noise;
    }
    let outlier = 10usize;
    y[outlier] += 40.0;

    let fit = OlsFit::new(x, y).unwrap();
    let d = cooks_distance(&fit);
    let f = dffits(&fit);
    let t = externally_studentized_residuals(&fit);

    // The outlier is the single most influential point on every measure.
    let argmax = |v: &Array1<f64>| {
        (0..n)
            .max_by(|&a, &b| v[a].abs().partial_cmp(&v[b].abs()).unwrap())
            .unwrap()
    };
    assert_eq!(
        argmax(&d),
        outlier,
        "Cook's distance should peak at outlier"
    );
    assert_eq!(argmax(&f), outlier, "DFFITS should peak at outlier");
    assert_eq!(
        argmax(&t),
        outlier,
        "studentized residual should peak at outlier"
    );
    // And it clears the common Cook's-distance 4/n flag comfortably.
    assert!(d[outlier] > 4.0 / n as f64);
}

#[test]
fn jarque_bera_detects_skew() {
    // Strongly right-skewed sample (a tight cluster near 0 with a long right
    // tail) -> JB large, small p-value.
    let mut skewed = vec![0.0; 24];
    for (k, v) in skewed.iter_mut().enumerate() {
        *v = (k as f64 % 4.0) * 0.1; // small symmetric-ish jitter in [0, 0.3]
    }
    skewed.extend_from_slice(&[8.0, 10.0, 12.0, 15.0, 20.0, 25.0]);
    let jb = jarque_bera(&mean_model(&skewed));
    assert!(jb.skewness > 1.0, "skewness = {}", jb.skewness);
    assert!(jb.p_value < 0.05, "JB p-value = {}", jb.p_value);
}

#[test]
fn jarque_bera_quiet_on_symmetric() {
    // Symmetric, roughly mesokurtic sample -> JB small, large p-value.
    let sym = [
        -3.0, -2.0, -1.5, -1.0, -0.5, 0.0, 0.0, 0.5, 1.0, 1.5, 2.0, 3.0,
    ];
    let jb = jarque_bera(&mean_model(&sym));
    assert!(jb.skewness.abs() < 1e-6, "skewness = {}", jb.skewness);
    assert!(jb.p_value > 0.2, "JB p-value = {}", jb.p_value);
}

#[test]
fn qq_data_near_diagonal_for_normalish_sample() {
    // A symmetric, bell-ish sample: theoretical vs sample quantiles should be
    // highly correlated (points hug the reference line).
    let sample = [
        -1.8, -1.1, -0.7, -0.4, -0.2, 0.0, 0.0, 0.2, 0.4, 0.7, 1.1, 1.8,
    ];
    let fit = mean_model(&sample);
    let r = standardized_residuals(&fit);
    let pairs = qq_plot_data(r.view());

    let n = pairs.len() as f64;
    let (tx, ty): (Vec<f64>, Vec<f64>) = pairs.iter().cloned().unzip();
    let mx = tx.iter().sum::<f64>() / n;
    let my = ty.iter().sum::<f64>() / n;
    let cov: f64 = tx.iter().zip(&ty).map(|(a, b)| (a - mx) * (b - my)).sum();
    let vx: f64 = tx.iter().map(|a| (a - mx).powi(2)).sum();
    let vy: f64 = ty.iter().map(|b| (b - my).powi(2)).sum();
    let corr = cov / (vx.sqrt() * vy.sqrt());
    assert!(corr > 0.97, "QQ correlation = {corr}");
}

#[test]
fn standardized_coefficients_rank_dominant_predictor() {
    // x1 moves y ten times harder than x2 (per unit of their own spread).
    let n = 12usize;
    let mut x = Array2::<f64>::ones((n, 3));
    let mut y = Array1::<f64>::zeros(n);
    for i in 0..n {
        let a = (i % 6) as f64; // predictor 1
        let b = ((i * 7) % 5) as f64; // predictor 2
        x[(i, 1)] = a;
        x[(i, 2)] = b;
        y[i] = 10.0 * a + 1.0 * b + 3.0;
    }
    let fit = OlsFit::new(x, y).unwrap();
    let beta = standardized_coefficients(&fit);
    assert!(beta[0].is_nan()); // intercept
    assert!(
        beta[1].abs() > beta[2].abs(),
        "x1 should dominate: {} vs {}",
        beta[1],
        beta[2]
    );
}

#[test]
fn studentized_residual_forms_are_consistent() {
    // On a well-behaved fit, internal and external studentized residuals are
    // close but not identical, and both are finite.
    let n = 15usize;
    let mut x = Array2::<f64>::ones((n, 2));
    let mut y = Array1::<f64>::zeros(n);
    for i in 0..n {
        let xi = (i + 1) as f64;
        x[(i, 1)] = xi;
        y[i] = 1.0 + 0.5 * xi + if i % 2 == 0 { 0.3 } else { -0.3 };
    }
    let fit = OlsFit::new(x, y).unwrap();
    let ri = internally_studentized_residuals(&fit);
    let re = externally_studentized_residuals(&fit);
    assert!(ri.iter().all(|v| v.is_finite()));
    assert!(re.iter().all(|v| v.is_finite()));

    // The two forms are tied by the exact identity
    //   t_i = r_i * sqrt((df - 1) / (df - r_i²)),   df = n - p,
    // which is the definitional relationship between external and internal
    // studentization. Verifying it pins both implementations against each other.
    let df = (fit.n_observations() - fit.n_parameters()) as f64;
    for i in 0..n {
        let expected = ri[i] * ((df - 1.0) / (df - ri[i] * ri[i])).sqrt();
        assert!(
            (re[i] - expected).abs() < 1e-9,
            "rstudent[{i}] = {} but identity gives {expected}",
            re[i]
        );
        // Same sign as the raw residual.
        assert_eq!(re[i].signum(), ri[i].signum());
    }
}