regression-diagnostics 0.2.0

//! Logistic-regression diagnostics: fit by IRLS, print the coefficient table and
//! goodness-of-fit summary, run Hosmer–Lemeshow, and list the most influential
//! observations.
//!
//! Run with: `cargo run --example logistic_diagnostics`

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

fn main() {
    // Synthetic admissions-style data: admit ~ gpa, exam_score. Two predictors
    // with genuine (but noisy) signal, labels overlap so the fit is non-trivial.
    let n = 50usize;
    let mut x = Array2::<f64>::ones((n, 3));
    let mut y = Array1::<f64>::zeros(n);
    for i in 0..n {
        let gpa = 2.0 + (i as f64 % 11.0) / 5.0; // ~2.0 .. 4.0
        let exam = 40.0 + ((i * 7) % 50) as f64; // ~40 .. 90
        x[(i, 1)] = gpa;
        x[(i, 2)] = exam;
        // Latent score with a deterministic "noise" flip pattern.
        let latent = -14.0 + 2.2 * gpa + 0.10 * exam;
        let flip = if i % 5 == 0 { -1.0 } else { 1.0 };
        y[i] = if latent * flip > 0.0 { 1.0 } else { 0.0 };
    }

    let fit = LogisticFit::new(x, y).expect("converges");

    let se = fit.coefficient_standard_errors();
    let z = fit.z_values();
    let pv = fit.p_values();
    let names = ["const", "gpa", "exam"];

    println!(
        "Logistic regression (IRLS, {} iterations)\n",
        fit.iterations()
    );
    println!(
        "{:<8}{:>10}{:>10}{:>9}{:>9}",
        "", "coef", "std err", "z", "P>|z|"
    );
    println!("{:-<46}", "");
    for j in 0..fit.n_parameters() {
        println!(
            "{:<8}{:>10.4}{:>10.4}{:>9.3}{:>9.3}",
            names[j],
            fit.coefficients()[j],
            se[j],
            z[j],
            pv[j]
        );
    }

    let gof = fit.goodness_of_fit();
    println!("\nGoodness of fit:");
    println!(
        "  Null deviance:     {:>8.3}  (df {})",
        gof.null_deviance, gof.df_null as usize
    );
    println!(
        "  Residual deviance: {:>8.3}  (df {})",
        gof.residual_deviance, gof.df_residual as usize
    );
    println!("  McFadden pseudo-R²:{:>8.4}", gof.mcfadden_r2);
    println!("  AIC / BIC:         {:>8.2} / {:.2}", gof.aic, gof.bic);

    let hl = fit.hosmer_lemeshow(10);
    println!(
        "\nHosmer–Lemeshow: Ĥ = {:.3}, df = {}, p = {:.3}  ({})",
        hl.statistic,
        hl.df,
        hl.p_value,
        if hl.p_value < 0.05 {
            "poor fit"
        } else {
            "no evidence of poor fit"
        }
    );

    // Top-3 most influential observations by Cook's distance.
    let cd = cooks_distance(&fit);
    let lev = leverage(&fit);
    let mut idx: Vec<usize> = (0..fit.n_observations()).collect();
    idx.sort_by(|&a, &b| cd[b].partial_cmp(&cd[a]).unwrap());
    println!("\nMost influential observations (Cook's distance):");
    println!("  obs   leverage    cooks_d");
    for &i in idx.iter().take(3) {
        println!("  {:>3}   {:>8.4}   {:>8.4}", i, lev[i], cd[i]);
    }
}