Skip to main content

solow_stats/
contingency.rs

1//! Contingency-table analysis: chi-squared test of nominal association.
2//!
3//! Mirrors the reference `contingency_tables.Table.test_nominal_association`:
4//! the expected counts under independence are the outer product of the row and
5//! column marginals scaled by the grand total; the chi-squared statistic is the
6//! sum of Pearson contributions `(obs − exp)² / exp` with
7//! `(rows − 1)(cols − 1)` degrees of freedom (no continuity correction).
8
9use ndarray::Array2;
10use solow_distributions::chi2_sf;
11
12/// Result of a contingency-table chi-squared test of independence.
13#[derive(Debug, Clone)]
14pub struct ContingencyResult {
15    /// Pearson chi-squared statistic.
16    pub statistic: f64,
17    /// Degrees of freedom `(rows − 1)(cols − 1)`.
18    pub df: usize,
19    /// p-value from the chi-squared distribution.
20    pub pvalue: f64,
21    /// Expected cell counts under independence.
22    pub expected: Array2<f64>,
23}
24
25/// A two-way contingency table of observed counts.
26#[derive(Debug, Clone)]
27pub struct Table {
28    table: Array2<f64>,
29}
30
31impl Table {
32    /// Build a table from a matrix of observed counts.
33    pub fn new(table: Array2<f64>) -> Self {
34        Table { table }
35    }
36
37    /// Estimated marginal probability distributions `(row, col)`.
38    pub fn marginal_probabilities(&self) -> (Vec<f64>, Vec<f64>) {
39        let n = self.table.sum();
40        let row: Vec<f64> = self
41            .table
42            .sum_axis(ndarray::Axis(1))
43            .iter()
44            .map(|&v| v / n)
45            .collect();
46        let col: Vec<f64> = self
47            .table
48            .sum_axis(ndarray::Axis(0))
49            .iter()
50            .map(|&v| v / n)
51            .collect();
52        (row, col)
53    }
54
55    /// Fitted (expected) cell counts under the independence model.
56    pub fn fittedvalues(&self) -> Array2<f64> {
57        let (row, col) = self.marginal_probabilities();
58        let total = self.table.sum();
59        let (r, c) = self.table.dim();
60        let mut fit = Array2::<f64>::zeros((r, c));
61        for i in 0..r {
62            for j in 0..c {
63                fit[[i, j]] = total * row[i] * col[j];
64            }
65        }
66        fit
67    }
68
69    /// Chi-squared test of independence between rows and columns.
70    pub fn test_nominal_association(&self) -> ContingencyResult {
71        let expected = self.fittedvalues();
72        let (r, c) = self.table.dim();
73        let mut statistic = 0.0;
74        for i in 0..r {
75            for j in 0..c {
76                let e = expected[[i, j]];
77                let o = self.table[[i, j]];
78                let resid = (o - e) / e.sqrt();
79                statistic += resid * resid;
80            }
81        }
82        let df = (r - 1) * (c - 1);
83        let pvalue = chi2_sf(statistic, df as f64);
84        ContingencyResult {
85            statistic,
86            df,
87            pvalue,
88            expected,
89        }
90    }
91}
92
93#[cfg(test)]
94mod tests {
95    use super::*;
96    use ndarray::array;
97
98    #[test]
99    fn expected_row_sums_match_observed() {
100        let t = Table::new(array![[10.0, 20.0, 30.0], [6.0, 9.0, 17.0]]);
101        let exp = t.fittedvalues();
102        // Row marginals of expected equal those of observed.
103        for i in 0..2 {
104            let so: f64 = (0..3).map(|j| t.table[[i, j]]).sum();
105            let se: f64 = (0..3).map(|j| exp[[i, j]]).sum();
106            assert!((so - se).abs() < 1e-9);
107        }
108    }
109
110    #[test]
111    fn statistic_nonnegative() {
112        let t = Table::new(array![[10.0, 20.0, 30.0], [6.0, 9.0, 17.0]]);
113        let res = t.test_nominal_association();
114        assert!(res.statistic >= 0.0);
115        assert_eq!(res.df, 2);
116    }
117}