regression-diagnostics 0.2.0

//! Logistic-regression diagnostics. The correctness anchors are exact: the MLE
//! satisfies the score equations `Xᵀ(y − p) = 0`, and for a single binary
//! predictor the coefficients equal the log-odds and log-odds-ratio in closed
//! form.

use ndarray::{Array1, Array2};
use regression_diagnostics::logistic::{
    cooks_distance, deviance_residuals, leverage, pearson_residuals, LogisticFit,
};
use regression_diagnostics::RegressionError;

/// A moderately-separated but non-degenerate dataset: x drives the log-odds.
fn make_dataset() -> (Array2<f64>, Array1<f64>) {
    let n = 40usize;
    let mut x = Array2::<f64>::ones((n, 2));
    let mut y = Array1::<f64>::zeros(n);
    for i in 0..n {
        let xi = (i as f64) * 0.25 - 5.0;
        x[(i, 1)] = xi;
        // Deterministic labels that overlap in the middle (not separable).
        y[i] = if (i % 3 == 0) ^ (xi > 0.0) { 1.0 } else { 0.0 };
    }
    (x, y)
}

#[test]
fn score_equations_hold_at_mle() {
    let (x, y) = make_dataset();
    let fit = LogisticFit::new(x.clone(), y.clone()).unwrap();
    let p = fit.fitted_probabilities();
    // Xᵀ(y − p) must be ~0 componentwise at the maximum likelihood estimate.
    let resid: Array1<f64> = &y - &p;
    for j in 0..fit.n_parameters() {
        let g: f64 = (0..fit.n_observations())
            .map(|i| x[(i, j)] * resid[i])
            .sum();
        assert!(g.abs() < 1e-6, "score[{j}] = {g}");
    }
}

#[test]
fn single_binary_predictor_matches_log_odds_ratio() {
    // x = 0 group: 4 ones, 6 zeros -> odds 4/6. x = 1 group: 7 ones, 3 zeros.
    // Closed-form MLE: intercept = ln(4/6), slope = ln((7/3)/(4/6)) = ln(3.5).
    let mut rows: Vec<[f64; 2]> = Vec::new();
    let mut labels: Vec<f64> = Vec::new();
    for k in 0..10 {
        rows.push([1.0, 0.0]);
        labels.push(if k < 4 { 1.0 } else { 0.0 });
    }
    for k in 0..10 {
        rows.push([1.0, 1.0]);
        labels.push(if k < 7 { 1.0 } else { 0.0 });
    }
    let n = rows.len();
    let mut x = Array2::<f64>::zeros((n, 2));
    for i in 0..n {
        x[(i, 0)] = rows[i][0];
        x[(i, 1)] = rows[i][1];
    }
    let y = Array1::from(labels);

    let fit = LogisticFit::new(x, y).unwrap();
    let want_intercept = (4.0f64 / 6.0).ln();
    let want_slope = 3.5f64.ln();
    assert!(
        (fit.coefficients()[0] - want_intercept).abs() < 1e-6,
        "intercept {} vs {want_intercept}",
        fit.coefficients()[0]
    );
    assert!(
        (fit.coefficients()[1] - want_slope).abs() < 1e-6,
        "slope {} vs {want_slope}",
        fit.coefficients()[1]
    );
}

#[test]
fn deviance_residuals_square_to_residual_deviance() {
    let (x, y) = make_dataset();
    let fit = LogisticFit::new(x, y).unwrap();
    let dr = deviance_residuals(&fit);
    let sum_sq: f64 = dr.iter().map(|d| d * d).sum();
    let gof = fit.goodness_of_fit();
    assert!(
        (sum_sq - gof.residual_deviance).abs() < 1e-8,
        "Σdev² = {sum_sq} vs residual_deviance = {}",
        gof.residual_deviance
    );
    // Pearson residuals are all finite.
    assert!(pearson_residuals(&fit).iter().all(|v| v.is_finite()));
}

#[test]
fn goodness_of_fit_is_sane() {
    let (x, y) = make_dataset();
    let fit = LogisticFit::new(x, y).unwrap();
    let gof = fit.goodness_of_fit();
    assert!(gof.residual_deviance <= gof.null_deviance + 1e-9);
    assert!((0.0..1.0).contains(&gof.mcfadden_r2));
    assert!(gof.aic.is_finite() && gof.bic.is_finite());
}

#[test]
fn logistic_leverage_sums_to_parameters() {
    let (x, y) = make_dataset();
    let fit = LogisticFit::new(x, y).unwrap();
    let total: f64 = leverage(&fit).sum();
    assert!(
        (total - fit.n_parameters() as f64).abs() < 1e-6,
        "Σh = {total}, expected {}",
        fit.n_parameters()
    );
    // Cook's distances are finite and non-negative.
    assert!(cooks_distance(&fit)
        .iter()
        .all(|v| v.is_finite() && *v >= 0.0));
}

#[test]
fn hosmer_lemeshow_runs() {
    let (x, y) = make_dataset();
    let fit = LogisticFit::new(x, y).unwrap();
    let hl = fit.hosmer_lemeshow(5);
    assert!(hl.statistic.is_finite());
    assert!(hl.p_value.is_finite());
    assert!(hl.df >= 1);
}

#[test]
fn perfect_separation_is_reported_not_hidden() {
    // y is exactly determined by the sign of x -> perfectly separable ->
    // coefficients diverge -> NotConverged, rather than a huge bogus fit.
    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 as f64) - 9.5;
        x[(i, 1)] = xi;
        y[i] = if xi > 0.0 { 1.0 } else { 0.0 };
    }
    let err = LogisticFit::new(x, y).unwrap_err();
    assert!(matches!(err, RegressionError::NotConverged { .. }));
}

#[test]
fn invalid_response_is_rejected() {
    let x = Array2::<f64>::ones((4, 2));
    // Non-binary value.
    let bad = Array1::from(vec![0.0, 1.0, 2.0, 1.0]);
    assert!(matches!(
        LogisticFit::new(x.clone(), bad).unwrap_err(),
        RegressionError::InvalidResponse { .. }
    ));
    // Single class.
    let one_class = Array1::from(vec![1.0, 1.0, 1.0, 1.0]);
    assert!(matches!(
        LogisticFit::new(x, one_class).unwrap_err(),
        RegressionError::InvalidResponse { .. }
    ));
}