regression-diagnostics 0.2.0

//! Multinomial and ordinal logistic regression on a three-category outcome:
//! fit both, print the coefficient tables and goodness-of-fit, and contrast what
//! each says — the multinomial gives a separate slope per class, the ordinal a
//! single proportional-odds slope with ordered thresholds.
//!
//! Run with: `cargo run --example categorical_diagnostics`

use ndarray::{Array1, Array2};
use regression_diagnostics::categorical::{MultinomialFit, OrdinalFit};

fn main() {
    // Ordered outcome (0 = low, 1 = medium, 2 = high) driven by one covariate.
    let n = 90usize;
    let mut y = Array1::<f64>::zeros(n);
    // Multinomial design includes an intercept column; ordinal must not.
    let mut x_multi = Array2::<f64>::ones((n, 2));
    let mut x_ord = Array2::<f64>::zeros((n, 1));
    for i in 0..n {
        let s = (i as f64) * 0.1 - 4.5; // covariate ranges roughly -4.5..4.5
        x_multi[(i, 1)] = s;
        x_ord[(i, 0)] = s;
        let base = if s < -1.5 {
            0
        } else if s < 1.5 {
            1
        } else {
            2
        };
        // Deterministic "noise" nudging some observations up a level.
        y[i] = ((base + (i % 6 == 0) as usize).min(2)) as f64;
    }

    // --- Multinomial (baseline-category) ----------------------------------
    let multi = MultinomialFit::new(x_multi, y.clone()).unwrap();
    println!("Multinomial logit (baseline = class 0), {} iters:", multi.iterations());
    let se = multi.coefficient_standard_errors();
    let z = multi.z_values();
    println!("{:<14}{:>10}{:>10}{:>9}", "", "coef", "std err", "z");
    for kk in 1..multi.n_classes() {
        for (a, name) in ["const", "x"].iter().enumerate() {
            println!(
                "class {kk} {:<6}{:>10.3}{:>10.3}{:>9.2}",
                name,
                multi.coefficients()[(kk - 1, a)],
                se[(kk - 1, a)],
                z[(kk - 1, a)]
            );
        }
    }
    println!(
        "  residual deviance = {:.2}, McFadden R² = {:.3}, AIC = {:.2}",
        multi.residual_deviance(),
        multi.mcfadden_r2(),
        multi.aic()
    );

    // --- Ordinal (proportional odds) --------------------------------------
    let ord = OrdinalFit::new(x_ord, y).unwrap();
    println!("\nProportional-odds ordinal logit, {} iters:", ord.iterations());
    println!(
        "  slope β(x) = {:.3}  (se {:.3}, z {:.2}, p {:.2e})",
        ord.coefficients()[0],
        ord.coefficient_standard_errors()[0],
        ord.z_values()[0],
        ord.p_values()[0]
    );
    print!("  thresholds α = [");
    for (k, a) in ord.thresholds().iter().enumerate() {
        print!("{}{:.3}", if k > 0 { ", " } else { "" }, a);
    }
    println!("]  (must be increasing)");
    println!(
        "  residual deviance = {:.2}, McFadden R² = {:.3}, AIC = {:.2}",
        ord.residual_deviance(),
        ord.mcfadden_r2(),
        ord.aic()
    );
    println!(
        "\n  The ordinal model spends {} parameters vs the multinomial's {} — one\n  \
         shared slope plus ordered cutpoints, rather than a slope per class.",
        ord.n_parameters(),
        multi.n_parameters()
    );
}