modewise_stats/
permutation.rs

1//! Permutation inference with family-wise error by the maximum statistic.
2
3use modewise_method::MethodError;
4use ndarray::{Array1, Array2, ArrayView2, Axis};
5use rand::SeedableRng;
6use rand::seq::SliceRandom;
7use rand_chacha::ChaCha8Rng;
8use rayon::prelude::*;
9
10use crate::dist::t_upper;
11use crate::glm::{Glm, welch_df, welch_t};
12use crate::nbs::{Component, Edge, components, largest};
13
14/// Which differences count as evidence.
15#[derive(Debug, Clone, Copy, PartialEq, Eq)]
16pub enum Tail {
17    /// The tested effect is positive.
18    Greater,
19    /// The tested effect is negative.
20    Less,
21    /// Either sign.
22    Both,
23}
24
25/// Rearranged statistics within this relative distance of the observed one count as ties (as
26/// large). Rearrangements that are equal in exact arithmetic, such as flipping every sign in a
27/// two-sided test, must not be split by rounding.
28const TIE: f64 = 1e-12;
29
30fn at_least(value: f64, observed: f64) -> bool {
31    value >= observed - TIE * observed.abs().max(1.0)
32}
33
34impl Tail {
35    fn stat(self, t: f64) -> f64 {
36        match self {
37            Self::Greater => t,
38            Self::Less => -t,
39            Self::Both => t.abs(),
40        }
41    }
42}
43
44/// What the observations may be rearranged by under the null hypothesis.
45#[derive(Debug, Clone, Copy, PartialEq, Eq)]
46pub enum Exchange {
47    /// Observations are exchangeable: shuffle them.
48    Permute,
49    /// Errors are symmetric around zero: flip signs (paired differences).
50    SignFlip,
51}
52
53/// The statistic tested.
54#[derive(Debug, Clone)]
55pub enum Statistic {
56    /// t of one regressor of a GLM; nuisance regressors are handled by Freedman–Lane.
57    Glm(Glm),
58    /// Welch's t of group a against the rest, for two groups without covariates.
59    Welch { in_a: Vec<bool> },
60}
61
62#[derive(Debug, Clone)]
63pub struct PermutationTest {
64    pub statistic: Statistic,
65    pub exchange: Exchange,
66    pub tail: Tail,
67    pub n_permutations: usize,
68    pub seed: u64,
69    /// Network-based statistic over edge features, at this primary threshold of the tailed t.
70    pub nbs: Option<(Vec<Option<Edge>>, f64)>,
71}
72
73/// Per-feature results.
74#[derive(Debug, Clone, PartialEq)]
75pub struct TestResult {
76    pub t: Array1<f64>,
77    pub df: Array1<f64>,
78    /// From Student's t (Welch–Satterthwaite df for Welch).
79    pub p_parametric: Array1<f64>,
80    /// (b + 1) / (m + 1) over m rearrangements (Phipson & Smyth 2010); exact when enumerated.
81    pub p_permutation: Array1<f64>,
82    /// Family-wise error over all features, by the null distribution of the maximum statistic.
83    pub p_fwer: Array1<f64>,
84    pub components: Vec<Component>,
85    pub n_rearrangements: usize,
86    pub exhaustive: bool,
87}
88
89impl PermutationTest {
90    fn t(&self, y: ArrayView2<f64>) -> Array1<f64> {
91        match &self.statistic {
92            Statistic::Glm(glm) => glm.t(y),
93            Statistic::Welch { in_a } => welch_t(y, in_a),
94        }
95    }
96
97    /// Runs the test on `y` (observations × features). Results do not depend on the number of
98    /// threads: each rearrangement draws from its own random stream.
99    pub fn run(&self, y: ArrayView2<f64>) -> Result<TestResult, MethodError> {
100        let n = y.nrows();
101        let expected = match &self.statistic {
102            Statistic::Glm(glm) => glm.n_observations(),
103            Statistic::Welch { in_a } => in_a.len(),
104        };
105        if n != expected {
106            return Err(MethodError::ShapeMismatch(format!(
107                "{n} observations but the design has {expected}"
108            )));
109        }
110        if self.n_permutations == 0 {
111            return Err(MethodError::invalid("n_permutations", "must be at least 1"));
112        }
113        if let (Statistic::Welch { .. }, Exchange::SignFlip) = (&self.statistic, self.exchange) {
114            return Err(MethodError::invalid(
115                "exchange",
116                "Welch's t compares groups; use permutation",
117            ));
118        }
119        let base = match &self.statistic {
120            Statistic::Glm(glm) => glm.nuisance_residuals(y),
121            Statistic::Welch { .. } => y.to_owned(),
122        };
123        let t = self.t(y);
124        let observed: Vec<f64> = t.iter().map(|&v| self.tail.stat(v)).collect();
125        let exhaustive = self.exchange == Exchange::SignFlip
126            && n < 63
127            && (1u64 << n) - 1 <= self.n_permutations as u64;
128        let m = if exhaustive {
129            (1usize << n) - 1
130        } else {
131            self.n_permutations
132        };
133
134        let rearranged = |r: usize| -> Array2<f64> {
135            match (self.exchange, exhaustive) {
136                (Exchange::SignFlip, true) => {
137                    let flips = r + 1; // 1..2^n − 1; 0 would be the observed data
138                    let mut out = base.clone();
139                    for (i, mut row) in out.axis_iter_mut(Axis(0)).enumerate() {
140                        if flips >> i & 1 == 1 {
141                            row.mapv_inplace(|v| -v);
142                        }
143                    }
144                    out
145                }
146                (Exchange::SignFlip, false) => {
147                    let mut rng = ChaCha8Rng::seed_from_u64(self.seed);
148                    rng.set_stream(r as u64);
149                    let mut out = base.clone();
150                    for mut row in out.axis_iter_mut(Axis(0)) {
151                        if rand::Rng::random::<bool>(&mut rng) {
152                            row.mapv_inplace(|v| -v);
153                        }
154                    }
155                    out
156                }
157                (Exchange::Permute, _) => {
158                    let mut rng = ChaCha8Rng::seed_from_u64(self.seed);
159                    rng.set_stream(r as u64);
160                    let mut order: Vec<usize> = (0..n).collect();
161                    order.shuffle(&mut rng);
162                    base.select(Axis(0), &order)
163                }
164            }
165        };
166
167        let features = observed.len();
168        let null: Vec<(Vec<bool>, f64, usize)> = (0..m)
169            .into_par_iter()
170            .map(|r| {
171                let stats: Vec<f64> = self
172                    .t(rearranged(r).view())
173                    .iter()
174                    .map(|&v| self.tail.stat(v))
175                    .collect();
176                let exceeds = stats
177                    .iter()
178                    .zip(&observed)
179                    .map(|(s, o)| at_least(*s, *o))
180                    .collect();
181                let max = stats
182                    .iter()
183                    .copied()
184                    .filter(|v| !v.is_nan())
185                    .fold(f64::NEG_INFINITY, f64::max);
186                let biggest = self
187                    .nbs
188                    .as_ref()
189                    .map_or(0, |(edges, threshold)| largest(&stats, edges, *threshold));
190                (exceeds, max, biggest)
191            })
192            .collect();
193
194        let ratio = |count: usize| (count + 1) as f64 / (m + 1) as f64;
195        let nan_or = |o: f64, value: f64| if o.is_nan() { f64::NAN } else { value };
196        let p_permutation = Array1::from_shape_fn(features, |f| {
197            nan_or(
198                observed[f],
199                ratio(null.iter().filter(|(e, _, _)| e[f]).count()),
200            )
201        });
202        let p_fwer = Array1::from_shape_fn(features, |f| {
203            nan_or(
204                observed[f],
205                ratio(
206                    null.iter()
207                        .filter(|(_, max, _)| at_least(*max, observed[f]))
208                        .count(),
209                ),
210            )
211        });
212        let found = match &self.nbs {
213            Some((edges, threshold)) => components(&observed, edges, *threshold)
214                .into_iter()
215                .map(|(network, edges)| {
216                    let size = edges.len();
217                    Component {
218                        network,
219                        edges,
220                        p: ratio(null.iter().filter(|(_, _, b)| *b >= size).count()),
221                    }
222                })
223                .collect(),
224            None => Vec::new(),
225        };
226
227        let df = match &self.statistic {
228            Statistic::Glm(glm) => Array1::from_elem(features, glm.degrees_of_freedom()),
229            Statistic::Welch { in_a } => welch_df(y, in_a),
230        };
231        let p_parametric = Array1::from_shape_fn(features, |f| {
232            let (tf, d) = (t[f], df[f]);
233            match self.tail {
234                Tail::Greater => t_upper(tf, d),
235                Tail::Less => t_upper(-tf, d),
236                Tail::Both => (2.0 * t_upper(tf.abs(), d)).min(1.0),
237            }
238        });
239        Ok(TestResult {
240            t,
241            df,
242            p_parametric,
243            p_permutation,
244            p_fwer,
245            components: found,
246            n_rearrangements: m,
247            exhaustive,
248        })
249    }
250}
modewise_stats/permutation.rs

modewise_stats/
permutation.rs
