regression_diagnostics/logistic/
goodness.rs

1use statrs::distribution::{ChiSquared, ContinuousCDF};
2
3use super::LogisticFit;
4
5/// Result of the Hosmer–Lemeshow goodness-of-fit test.
6#[derive(Debug, Clone, Copy, PartialEq)]
7pub struct HosmerLemeshow {
8    /// The Ĥ statistic.
9    pub statistic: f64,
10    /// Degrees of freedom `g − 2` (`g` = number of groups).
11    pub df: usize,
12    /// Upper-tail p-value under χ²(df). A **small** p-value indicates **poor**
13    /// fit (observed and expected event counts disagree across risk groups).
14    pub p_value: f64,
15    /// Number of groups actually used.
16    pub groups: usize,
17}
18
19/// Bundle of overall logistic goodness-of-fit statistics.
20#[derive(Debug, Clone, Copy, PartialEq)]
21pub struct GoodnessOfFit {
22    /// Deviance of the intercept-only (null) model, `−2ℓ₀`.
23    pub null_deviance: f64,
24    /// Residual deviance of the fitted model, `−2ℓ`.
25    pub residual_deviance: f64,
26    /// Degrees of freedom of the null deviance (`n − 1`).
27    pub df_null: f64,
28    /// Degrees of freedom of the residual deviance (`n − p`).
29    pub df_residual: f64,
30    /// McFadden's pseudo-R², `1 − ℓ/ℓ₀`.
31    pub mcfadden_r2: f64,
32    /// Akaike information criterion, `deviance + 2p`.
33    pub aic: f64,
34    /// Bayesian information criterion, `deviance + ln(n)·p`.
35    pub bic: f64,
36}
37
38impl LogisticFit {
39    /// Overall goodness-of-fit summary: null/residual deviance, McFadden's
40    /// pseudo-R², and AIC/BIC.
41    ///
42    /// The **residual deviance** is `−2ℓ` (the saturated-model log-likelihood is
43    /// zero for ungrouped binary data), and equals the sum of squared deviance
44    /// residuals. **McFadden's pseudo-R²** compares the fitted log-likelihood to
45    /// the intercept-only model; it is bounded in `[0, 1)` but runs lower than an
46    /// OLS R² for comparable fits, so judge it on that scale.
47    pub fn goodness_of_fit(&self) -> GoodnessOfFit {
48        let n = self.n_observations() as f64;
49        let p = self.n_parameters() as f64;
50        let y = self.response();
51
52        let residual_deviance = -2.0 * self.log_likelihood();
53
54        // Null model: constant probability = mean(y).
55        let ybar = y.sum() / n;
56        let ll_null: f64 = y
57            .iter()
58            .map(|&yi| yi * ybar.ln() + (1.0 - yi) * (1.0 - ybar).ln())
59            .sum();
60        let null_deviance = -2.0 * ll_null;
61
62        let mcfadden_r2 = if ll_null != 0.0 {
63            1.0 - self.log_likelihood() / ll_null
64        } else {
65            f64::NAN
66        };
67
68        GoodnessOfFit {
69            null_deviance,
70            residual_deviance,
71            df_null: n - 1.0,
72            df_residual: n - p,
73            mcfadden_r2,
74            aic: residual_deviance + 2.0 * p,
75            bic: residual_deviance + n.ln() * p,
76        }
77    }
78
79    /// Hosmer–Lemeshow goodness-of-fit test with `groups` risk deciles
80    /// (`groups = 10` is the customary choice).
81    ///
82    /// Observations are ordered by fitted probability and split into `groups`
83    /// near-equal bins; the statistic compares observed and expected event counts
84    /// per bin,
85    ///
86    /// `Ĥ = Σ_g (O_g − E_g)² / (n_g · π̄_g · (1 − π̄_g))`,
87    ///
88    /// which is asymptotically χ²(groups − 2). Unlike most tests here, a **small
89    /// p-value means the model fits poorly**. Bins whose mean probability is
90    /// exactly 0 or 1 contribute nothing (their variance term is undefined).
91    ///
92    /// The test is only meaningful when `groups ≥ 3` and there are enough
93    /// observations to populate the bins; with fewer than `groups + 1`
94    /// observations the statistic is returned as `NaN`.
95    pub fn hosmer_lemeshow(&self, groups: usize) -> HosmerLemeshow {
96        let n = self.n_observations();
97        let y = self.response();
98        let p = self.fitted_probabilities();
99
100        if groups < 3 || n < groups + 1 {
101            return HosmerLemeshow {
102                statistic: f64::NAN,
103                df: groups.saturating_sub(2),
104                p_value: f64::NAN,
105                groups,
106            };
107        }
108
109        // Order indices by fitted probability.
110        let mut order: Vec<usize> = (0..n).collect();
111        order.sort_by(|&a, &b| p[a].partial_cmp(&p[b]).unwrap_or(std::cmp::Ordering::Equal));
112
113        let mut statistic = 0.0;
114        let mut used_groups = 0usize;
115        for g in 0..groups {
116            let start = g * n / groups;
117            let end = (g + 1) * n / groups;
118            if end <= start {
119                continue;
120            }
121            let ng = (end - start) as f64;
122            let mut observed = 0.0;
123            let mut expected = 0.0;
124            for &idx in &order[start..end] {
125                observed += y[idx];
126                expected += p[idx];
127            }
128            let pbar = expected / ng;
129            if pbar <= 0.0 || pbar >= 1.0 {
130                continue;
131            }
132            let diff = observed - expected;
133            statistic += diff * diff / (ng * pbar * (1.0 - pbar));
134            used_groups += 1;
135        }
136
137        let df = used_groups.saturating_sub(2);
138        let p_value = match ChiSquared::new(df as f64) {
139            Ok(dist) if df >= 1 && statistic.is_finite() => 1.0 - dist.cdf(statistic),
140            _ => f64::NAN,
141        };
142
143        HosmerLemeshow {
144            statistic,
145            df,
146            p_value,
147            groups: used_groups,
148        }
149    }
150}
regression_diagnostics/logistic/goodness.rs

regression_diagnostics/logistic/
goodness.rs
