regression_diagnostics/residuals/white_test.rs
1use super::{lm_heteroskedasticity_test, non_intercept_columns};
2use crate::OlsFit;
3
4/// Result of White's heteroskedasticity test (LM form).
5#[derive(Debug, Clone, Copy, PartialEq)]
6pub struct WhiteTest {
7 /// Lagrange-multiplier statistic `n·R²_aux`.
8 pub statistic: f64,
9 /// Degrees of freedom (number of auxiliary regressors).
10 pub df: usize,
11 /// Upper-tail p-value under χ²(df).
12 pub p_value: f64,
13}
14
15/// White's test for heteroskedasticity.
16///
17/// Like [`breusch_pagan`](super::breusch_pagan) it regresses the squared
18/// residuals on auxiliary regressors, but the auxiliary set is richer: the
19/// original regressors, their **squares**, and all **pairwise cross-products**.
20/// This makes White's test sensitive to more general (including nonlinear) forms
21/// of heteroskedasticity and, because it includes cross-products, it doubles as a
22/// specification check.
23///
24/// # Tradeoff versus Breusch-Pagan
25///
26/// The extra terms cost degrees of freedom. In **small samples** that spreads the
27/// test's power thin, so White's can fail to flag heteroskedasticity that the
28/// tighter Breusch-Pagan catches. White's is the more general test; Breusch-Pagan
29/// has more power when you already suspect a linear variance relationship. The
30/// crate offers both rather than picking one for you.
31///
32/// With `k` non-intercept regressors the auxiliary regression has `k` linear
33/// terms, `k` squares, and `k(k−1)/2` cross-products; duplicate columns that can
34/// arise (e.g. the square of a 0/1 dummy equals the dummy) are dropped to keep
35/// the auxiliary design full rank.
36///
37/// Because that regressor count grows quadratically in `k`, small samples can
38/// leave the auxiliary regression without enough residual degrees of freedom to
39/// estimate. When there are too few observations the statistic and p-value are
40/// returned as `NaN` rather than a fabricated number.
41pub fn white_test(fit: &OlsFit) -> WhiteTest {
42 let n = fit.n_observations();
43 let cols = non_intercept_columns(fit);
44 let x = fit.design_matrix();
45
46 let base: Vec<Vec<f64>> = cols.iter().map(|&j| x.column(j).to_vec()).collect();
47 let k = base.len();
48
49 let mut regressors: Vec<Vec<f64>> = Vec::new();
50 // Linear terms.
51 for col in &base {
52 regressors.push(col.clone());
53 }
54 // Squares.
55 for col in &base {
56 regressors.push(col.iter().map(|v| v * v).collect());
57 }
58 // Pairwise cross-products.
59 for a in 0..k {
60 for b in (a + 1)..k {
61 let cross: Vec<f64> = base[a]
62 .iter()
63 .zip(base[b].iter())
64 .map(|(u, v)| u * v)
65 .collect();
66 regressors.push(cross);
67 }
68 }
69
70 // Drop near-duplicate / near-constant auxiliary columns to protect the
71 // auxiliary design's rank (e.g. the square of a 0/1 dummy).
72 let regressors = dedup_regressors(n, regressors);
73
74 let resid_sq: Vec<f64> = fit.residuals().iter().map(|e| e * e).collect();
75
76 let (statistic, df, p_value) = lm_heteroskedasticity_test(n, ®ressors, &resid_sq);
77 WhiteTest {
78 statistic,
79 df,
80 p_value,
81 }
82}
83
84/// Remove constant columns and columns that duplicate an earlier one (up to a
85/// tight relative tolerance), which would otherwise make the auxiliary design
86/// rank-deficient.
87fn dedup_regressors(n: usize, cols: Vec<Vec<f64>>) -> Vec<Vec<f64>> {
88 let mut kept: Vec<Vec<f64>> = Vec::new();
89 for col in cols {
90 let first = col[0];
91 let is_constant = col
92 .iter()
93 .all(|&v| (v - first).abs() <= 1e-12 * first.abs().max(1.0));
94 if is_constant {
95 continue;
96 }
97 let dup = kept
98 .iter()
99 .any(|k| (0..n).all(|i| (k[i] - col[i]).abs() <= 1e-12 * col[i].abs().max(1.0)));
100 if !dup {
101 kept.push(col);
102 }
103 }
104 kept
105}