modewise_stats/
fdr.rs

1//! False discovery rate adjustment.
2
3/// Which false discovery rate procedure.
4#[derive(Debug, Clone, Copy, PartialEq, Eq)]
5pub enum Fdr {
6    /// Benjamini & Hochberg (1995): valid for independent or positively dependent tests.
7    BenjaminiHochberg,
8    /// Benjamini & Yekutieli (2001): valid under any dependence, more conservative.
9    BenjaminiYekutieli,
10}
11
12/// Adjusted p-values (q-values) in the original order. NaN p-values stay NaN and are not counted.
13pub fn adjust(p: &[f64], method: Fdr) -> Vec<f64> {
14    let mut order: Vec<usize> = (0..p.len()).filter(|&i| !p[i].is_nan()).collect();
15    let m = order.len() as f64;
16    let scale = match method {
17        Fdr::BenjaminiHochberg => 1.0,
18        Fdr::BenjaminiYekutieli => (1..=order.len()).map(|i| 1.0 / i as f64).sum(),
19    };
20    order.sort_by(|&a, &b| p[a].total_cmp(&p[b]));
21    let mut q = vec![f64::NAN; p.len()];
22    let mut running = 1.0f64;
23    for (rank, &i) in order.iter().enumerate().rev() {
24        running = running.min(p[i] * m * scale / (rank + 1) as f64);
25        q[i] = running.min(1.0);
26    }
27    q
28}
29
30#[cfg(test)]
31mod tests {
32    use super::*;
33
34    #[test]
35    fn bh_and_by_by_hand() {
36        let p = [0.01, 0.04, 0.03, f64::NAN, 0.5];
37        let bh = adjust(&p, Fdr::BenjaminiHochberg);
38        // Sorted 0.01, 0.03, 0.04, 0.5 with m = 4: 0.04, 0.0533…, 0.0533…, 0.5.
39        assert!((bh[0] - 0.04).abs() < 1e-15);
40        assert!((bh[1] - 0.04 * 4.0 / 3.0).abs() < 1e-15 && (bh[2] - bh[1]).abs() < 1e-15);
41        assert!(bh[3].is_nan() && (bh[4] - 0.5).abs() < 1e-15);
42        let by = adjust(&p, Fdr::BenjaminiYekutieli);
43        let c: f64 = 1.0 + 0.5 + 1.0 / 3.0 + 0.25;
44        assert!((by[0] - (0.04 * c).min(1.0)).abs() < 1e-15);
45        assert_eq!(adjust(&[0.9, 0.95], Fdr::BenjaminiYekutieli), [1.0, 1.0]);
46        assert!(adjust(&[], Fdr::BenjaminiHochberg).is_empty());
47    }
48}
modewise_stats/fdr.rs

modewise_stats/
fdr.rs
