r2rs-stats 0.1.1

Statistics programming for Rust based on R's stats package
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

// "Whatever you do, work at it with all your heart, as working for the Lord,
// not for human masters, since you know that you will receive an inheritance
// from the Lord as a reward. It is the Lord Christ you are serving."
// (Col 3:23-24)

use r2rs_base::func::{order, pmin};

use crate::traits::StatArray;

#[derive(Copy, Clone, Debug, PartialEq)]
pub enum PAdjustMethod {
    FDR,
    Bonferroni,
    Hommel,
    Hochberg,
    BH,
    BY,
}

pub fn p_adjust(p: &[f64], mut method: Option<PAdjustMethod>) -> Vec<f64> {
    if p.len() <= 1 {
        p.to_vec()
    } else {
        if p.len() == 2 && method == Some(PAdjustMethod::Hommel) {
            method = Some(PAdjustMethod::Hochberg);
        }

        match method {
            None => p.to_vec(),
            Some(PAdjustMethod::FDR) | Some(PAdjustMethod::BH) => pmin(
                &[1.0],
                &p.iter()
                    .map(|&p_i| p.len() as f64 * p_i)
                    .collect::<Vec<_>>(),
            ),
            Some(PAdjustMethod::Bonferroni) => {
                let o = order(p);
                let ro = order(&o);
                let intermediate = o
                    .into_iter()
                    .enumerate()
                    .map(|(i, o)| (p.len() + 1 - i) as f64 * p[o])
                    .collect::<Vec<_>>()
                    .cummax();
                let mins = pmin(&[1.0], &intermediate);
                ro.into_iter().map(|ro| mins[ro]).collect()
            }
            Some(PAdjustMethod::Hommel) => {
                // let i = (0..p.len()).collect::<Vec<_>>();
                let o = order(p);
                let p = o.iter().map(|&o| p[o]).collect::<Vec<_>>();
                let ro = order(&o);
                let mut q = vec![
                    p.iter()
                        .enumerate()
                        .map(|(i, &p_i)| p.len() as f64 * p_i / i as f64)
                        .min_by(|f1, f2| f1.partial_cmp(f2).unwrap())
                        .unwrap();
                    p.len()
                ];
                let mut pa = q.clone();

                for j in (2..(p.len() - 1)).rev() {
                    let ij = (0..(p.len() - j + 1)).collect::<Vec<_>>();
                    let i2 = ((p.len() - j + 2)..p.len()).collect::<Vec<_>>();
                    let q1 = i2
                        .iter()
                        .zip(2..j)
                        .map(|(&i2, _2j)| j as f64 * p[i2] / _2j as f64)
                        .min_by(|f1, f2| f1.partial_cmp(f2).unwrap())
                        .unwrap();
                    for ij in ij {
                        q[ij] = q1.min(j as f64 * p[ij]);
                    }
                    for i2 in i2 {
                        q[i2] = q[p.len() - j - 1];
                    }
                    for i in 0..pa.len() {
                        pa[i] = pa[i].max(q[i]);
                    }
                }

                let mut ret = Vec::new();
                for i in 0..p.len() {
                    let j = ro[i];
                    ret.push(pa[j].max(p[j]));
                }
                ret
            }
            Some(PAdjustMethod::Hochberg) => {
                let o = order(p).into_iter().rev().collect::<Vec<_>>();
                let ro = order(&o);
                let intermediate = o
                    .into_iter()
                    .enumerate()
                    .map(|(i, o)| (p.len() + 1 - i) as f64 * p[o])
                    .collect::<Vec<_>>()
                    .cummin();
                let mins = pmin(&[1.0], &intermediate);
                ro.into_iter().map(|ro| mins[ro]).collect()
            }
            Some(PAdjustMethod::BY) => {
                let o = order(p).into_iter().rev().collect::<Vec<_>>();
                let ro = order(&o);
                let q = (0..p.len()).map(|n| 1.0 / n as f64).sum::<f64>();
                let intermediate = o
                    .into_iter()
                    .enumerate()
                    .map(|(i, o)| q * (p.len() as f64 / i as f64) * p[o])
                    .collect::<Vec<_>>()
                    .cummin();
                let mins = pmin(&[1.0], &intermediate);
                ro.into_iter().map(|ro| mins[ro]).collect()
            }
        }
    }
}