regression-diagnostics 0.2.0

use ndarray::Array1;

use crate::OlsFit;

/// Standardized (beta) coefficients — the linear model's version of feature
/// importance.
///
/// Each raw coefficient is rescaled to the units it would take on a z-scored
/// design and z-scored response:
///
/// `βⱼ* = βⱼ · (sⱼ / s_y)`
///
/// where `sⱼ` is the standard deviation of predictor `j` and `s_y` that of the
/// response. Because the predictors are put on a common (unit-variance) scale,
/// the **relative magnitudes** of the standardized coefficients are comparable
/// across predictors in a way the raw coefficients — each in its own units — are
/// not.
///
/// This rescales the existing fit's coefficients rather than refitting; the two
/// are algebraically identical for OLS. Sample standard deviations use the `n − 1`
/// (unbiased) denominator.
///
/// # Return layout
///
/// Aligned to the design columns: entry `j` is the standardized coefficient of
/// column `j`. The **intercept column's entry is [`f64::NAN`]** — a standardized
/// intercept is not meaningful (it is zero by construction on centered data).
///
/// # Example
///
/// ```
/// use ndarray::array;
/// use regression_diagnostics::{OlsFit, coefficients::standardized_coefficients};
///
/// // x1 drives y an order of magnitude harder than x2 does.
/// let x = array![
///     [1.0, 1.0, 5.0],
///     [1.0, 2.0, 4.0],
///     [1.0, 3.0, 6.0],
///     [1.0, 4.0, 5.0],
///     [1.0, 5.0, 7.0],
/// ];
/// let y = array![10.0, 20.5, 29.5, 40.5, 50.0];
/// let fit = OlsFit::new(x, y).unwrap();
/// let b = standardized_coefficients(&fit);
/// assert!(b[0].is_nan());              // intercept
/// assert!(b[1].abs() > b[2].abs());    // x1 dominates
/// ```
pub fn standardized_coefficients(fit: &OlsFit) -> Array1<f64> {
    let n = fit.n_observations() as f64;
    let x = fit.design_matrix();
    let y = fit.response();

    let sd = |col: ndarray::ArrayView1<f64>| -> f64 {
        let mean = col.sum() / n;
        let var = col.iter().map(|v| (v - mean).powi(2)).sum::<f64>() / (n - 1.0);
        var.sqrt()
    };

    let s_y = sd(y);
    let coef = fit.coefficients();

    Array1::from_shape_fn(fit.n_parameters(), |j| {
        if fit.intercept_column() == Some(j) || s_y <= 0.0 {
            return f64::NAN;
        }
        let s_j = sd(x.column(j));
        coef[j] * s_j / s_y
    })
}