modewise_stats/
glm.rs

1//! Ordinary least squares for many features sharing one design.
2
3use modewise_method::MethodError;
4use ndarray::{Array1, Array2, ArrayView1, ArrayView2, Axis};
5
6/// A design matrix (observations × regressors) and the regressor being tested.
7#[derive(Debug, Clone)]
8pub struct Glm {
9    design: Array2<f64>,
10    tested: usize,
11    /// (MᵀM)⁻¹Mᵀ, so coefficients are one product per permutation.
12    pinv: Array2<f64>,
13    /// [(MᵀM)⁻¹]_jj for the tested regressor.
14    scale: f64,
15    /// Projection onto the nuisance regressors (all but the tested one); `None` if there are none.
16    nuisance: Option<Array2<f64>>,
17}
18
19impl Glm {
20    pub fn new(design: Array2<f64>, tested: usize) -> Result<Self, MethodError> {
21        let (n, p) = design.dim();
22        if tested >= p {
23            return Err(MethodError::invalid(
24                "tested",
25                format!("regressor {tested} does not exist; there are {p}"),
26            ));
27        }
28        if n <= p {
29            return Err(MethodError::invalid(
30                "design",
31                format!("{n} observations cannot fit {p} regressors with error left over"),
32            ));
33        }
34        let inverse = invert(&design.t().dot(&design))
35            .ok_or_else(|| MethodError::invalid("design", "regressors are collinear"))?;
36        let nuisance_columns: Vec<usize> = (0..p).filter(|&c| c != tested).collect();
37        let nuisance = if nuisance_columns.is_empty() {
38            None
39        } else {
40            let z = design.select(Axis(1), &nuisance_columns);
41            let zz = invert(&z.t().dot(&z)).ok_or_else(|| {
42                MethodError::invalid("design", "nuisance regressors are collinear")
43            })?;
44            Some(z.dot(&zz).dot(&z.t()))
45        };
46        Ok(Self {
47            pinv: inverse.dot(&design.t()),
48            scale: inverse[[tested, tested]],
49            design,
50            tested,
51            nuisance,
52        })
53    }
54
55    pub fn n_observations(&self) -> usize {
56        self.design.nrows()
57    }
58
59    pub fn degrees_of_freedom(&self) -> f64 {
60        (self.design.nrows() - self.design.ncols()) as f64
61    }
62
63    /// t of the tested coefficient for each feature (column of `y`).
64    pub fn t(&self, y: ArrayView2<f64>) -> Array1<f64> {
65        let beta = self.pinv.dot(&y);
66        let residuals = &y - &self.design.dot(&beta);
67        let df = self.degrees_of_freedom();
68        let rss = residuals.map_axis(Axis(0), |r| r.dot(&r));
69        Array1::from_shape_fn(y.ncols(), |f| {
70            beta[[self.tested, f]] / (rss[f] / df * self.scale).sqrt()
71        })
72    }
73
74    /// Residuals of each feature after the nuisance regressors (Freedman & Lane 1983).
75    pub fn nuisance_residuals(&self, y: ArrayView2<f64>) -> Array2<f64> {
76        match &self.nuisance {
77            Some(h) => &y - &h.dot(&y),
78            None => y.to_owned(),
79        }
80    }
81}
82
83/// Mean, variance of the mean, and size of one group in a column.
84fn group(column: ArrayView1<f64>, in_a: &[bool], member: bool) -> (f64, f64, f64) {
85    let values: Vec<f64> = column
86        .iter()
87        .zip(in_a)
88        .filter(|(_, g)| **g == member)
89        .map(|(v, _)| *v)
90        .collect();
91    let n = values.len() as f64;
92    let mean = values.iter().sum::<f64>() / n;
93    let var = values.iter().map(|v| (v - mean).powi(2)).sum::<f64>() / (n - 1.0);
94    (mean, var / n, n)
95}
96
97/// Welch's t of group `a` minus the rest, for each feature (column of `y`).
98pub fn welch_t(y: ArrayView2<f64>, in_a: &[bool]) -> Array1<f64> {
99    Array1::from_iter(y.columns().into_iter().map(|c| {
100        let ((ma, va, _), (mb, vb, _)) = (group(c, in_a, true), group(c, in_a, false));
101        (ma - mb) / (va + vb).sqrt()
102    }))
103}
104
105/// Welch–Satterthwaite degrees of freedom per feature.
106pub fn welch_df(y: ArrayView2<f64>, in_a: &[bool]) -> Array1<f64> {
107    Array1::from_iter(y.columns().into_iter().map(|c| {
108        let ((_, a, na), (_, b, nb)) = (group(c, in_a, true), group(c, in_a, false));
109        (a + b).powi(2) / (a * a / (na - 1.0) + b * b / (nb - 1.0))
110    }))
111}
112
113/// Inverse of a small symmetric positive definite matrix by Gauss–Jordan elimination with
114/// partial pivoting; `None` when (nearly) singular.
115fn invert(m: &Array2<f64>) -> Option<Array2<f64>> {
116    let n = m.nrows();
117    let mut a = m.clone();
118    let mut inv = Array2::eye(n);
119    let scale = m.iter().fold(0.0f64, |s, v| s.max(v.abs()));
120    for col in 0..n {
121        let pivot = (col..n).max_by(|&i, &j| a[[i, col]].abs().total_cmp(&a[[j, col]].abs()))?;
122        if a[[pivot, col]].abs() <= 1e-12 * scale {
123            return None;
124        }
125        for k in 0..n {
126            a.swap([pivot, k], [col, k]);
127            inv.swap([pivot, k], [col, k]);
128        }
129        let d = a[[col, col]];
130        for k in 0..n {
131            a[[col, k]] /= d;
132            inv[[col, k]] /= d;
133        }
134        for row in 0..n {
135            if row != col {
136                let f = a[[row, col]];
137                for k in 0..n {
138                    a[[row, k]] -= f * a[[col, k]];
139                    inv[[row, k]] -= f * inv[[col, k]];
140                }
141            }
142        }
143    }
144    Some(inv)
145}
146
147#[cfg(test)]
148mod tests {
149    use ndarray::array;
150
151    use super::*;
152
153    #[test]
154    fn inversion_and_rejections() {
155        let inv = invert(&array![[4.0, 2.0], [2.0, 3.0]]).unwrap();
156        assert!(
157            (&inv.dot(&array![[4.0, 2.0], [2.0, 3.0]]) - &Array2::<f64>::eye(2))
158                .iter()
159                .all(|v| v.abs() < 1e-14)
160        );
161        assert!(invert(&array![[1.0, 2.0], [2.0, 4.0]]).is_none());
162        let collinear = array![[1.0, 2.0], [1.0, 2.0], [1.0, 2.0], [1.0, 2.0]];
163        assert!(Glm::new(collinear, 1).is_err());
164        assert!(
165            Glm::new(Array2::ones((2, 2)), 0).is_err(),
166            "no error degrees of freedom"
167        );
168        assert!(Glm::new(Array2::ones((4, 1)), 1).is_err());
169        let with_nuisance = array![
170            [1.0, 0.0, 1.0],
171            [1.0, 1.0, 1.0],
172            [1.0, 0.0, 1.0],
173            [1.0, 1.0, 1.0],
174            [1.0, 0.5, 1.0]
175        ];
176        assert!(
177            Glm::new(with_nuisance, 1).is_err(),
178            "intercept duplicated in the nuisance"
179        );
180    }
181
182    #[test]
183    fn one_sample_t_by_hand() {
184        let y = array![[1.0], [2.0], [3.0], [6.0]];
185        let glm = Glm::new(Array2::ones((4, 1)), 0).unwrap();
186        // mean 3, sd sqrt(14/3), t = 3 / (sd / 2).
187        let expected = 3.0 / ((14.0f64 / 3.0).sqrt() / 2.0);
188        assert!((glm.t(y.view())[0] - expected).abs() < 1e-12);
189        assert_eq!(glm.nuisance_residuals(y.view()), y);
190        assert_eq!((glm.n_observations(), glm.degrees_of_freedom()), (4, 3.0));
191    }
192}
modewise_stats/glm.rs

modewise_stats/
glm.rs
