single_statistics/testing/effect/
mod.rs

1use nalgebra_sparse::CsrMatrix;
2use single_utilities::traits::{FloatOps, FloatOpsTS};
3
4/// Calculate log2 fold change between two groups
5pub fn calculate_log2_fold_change<T>(
6    matrix: &CsrMatrix<T>,
7    col: usize,
8    group1_indices: &[usize], // Group of interest
9    group2_indices: &[usize], // Reference group
10    pseudo_count: T,        // Small value like 1e-9 or 1.0
11) -> anyhow::Result<T>
12where
13    T: FloatOps,
14{
15    if group1_indices.is_empty() || group2_indices.is_empty() {
16        return Err(anyhow::anyhow!("Group indices cannot be empty"));
17    }
18
19    let mut sum1 = T::zero();
20    let mut _count1 = 0;
21    for &row in group1_indices {
22        if let Some(entry) = matrix.get_entry(row, col) {
23            let value = entry.into_value();
24            sum1 += value;
25            _count1 += 1;
26        }
27    }
28
29    let mut sum2 = T::zero();
30    let mut _count2 = 0;
31    for &row in group2_indices {
32        if let Some(entry) = matrix.get_entry(row, col) {
33            let value = entry.into_value();
34            sum2 += value;
35            _count2 += 1;
36        }
37    }
38    let go_f64 = T::from(group1_indices.len()).unwrap();
39    let gt_f64 = T::from(group2_indices.len()).unwrap();
40
41    let mean1 = sum1 / go_f64 + pseudo_count;
42    let mean2 = sum2 / gt_f64 + pseudo_count;
43
44    let log2_fc = (mean1 / mean2).log2();
45
46    Ok(log2_fc)
47}
48
49/// Calculate Cohen's d effect size for a row
50pub fn calculate_cohens_d<T>(
51    matrix: &CsrMatrix<T>,
52    row: usize,
53    group1_indices: &[usize],
54    group2_indices: &[usize],
55) -> anyhow::Result<T>
56where
57    T: FloatOpsTS,
58{
59    if group1_indices.len() < 2 || group2_indices.len() < 2 {
60        return Err(anyhow::anyhow!(
61            "Each group must have at least 2 samples for Cohen's d"
62        ));
63    }
64
65    // Extract values for group 1
66    let mut sum_g1 = T::zero();
67    let mut group1_values = Vec::with_capacity(group1_indices.len());
68    for &col in group1_indices {
69        if let Some(entry) = matrix.get_entry(row, col) {
70            let value = entry.into_value();
71            sum_g1 += value;
72            group1_values.push(value);
73        } else {
74            group1_values.push(T::zero());
75        }
76    }
77
78    // Extract values for group 2
79    let mut sum_g2 = T::zero();
80    let mut group2_values = Vec::with_capacity(group2_indices.len());
81    for &col in group2_indices {
82        if let Some(entry) = matrix.get_entry(row, col) {
83            let value = entry.into_value();
84            sum_g2 += value;
85            group2_values.push(value);
86        } else {
87            group2_values.push(T::zero());
88        }
89    }
90
91    let go_t = T::from(group1_indices.len()).unwrap();
92    let gt_t = T::from(group2_indices.len()).unwrap();
93
94    // Calculate means
95    let mean1 = sum_g1 / go_t;
96    let mean2 = sum_g2 / gt_t;
97    // Calculate variances
98    let var1 = group1_values
99        .iter()
100        .map(|&x| num_traits::Float::powi(x - mean1, 2))
101        .sum::<T>()
102        / (go_t - T::one());
103
104    let var2 = group2_values
105        .iter()
106        .map(|&x| num_traits::Float::powi(x - mean2, 2))
107        .sum::<T>()
108        / (gt_t - T::one());
109
110    // Calculate pooled standard deviation
111    let pooled_sd = (((go_t - T::one()) * var1 + (gt_t - T::one()) * var2) / (go_t + gt_t - T::from(2.0).unwrap())).sqrt();
112
113    // Calculate Cohen's d
114    let d = (mean2 - mean1) / pooled_sd;
115
116    Ok(d)
117}
118
119/// Calculate Hedge's g (bias-corrected effect size)
120pub fn calculate_hedges_g<T>(
121    matrix: &CsrMatrix<T>,
122    row: usize,
123    group1_indices: &[usize],
124    group2_indices: &[usize],
125) -> anyhow::Result<T>
126where
127    T: FloatOpsTS,
128{
129    // First calculate Cohen's d
130    let d = calculate_cohens_d(matrix, row, group1_indices, group2_indices)?;
131    let one = T::one();
132    let two = T::from(2.0).unwrap();
133    let three = T::from(3.0).unwrap();
134    let four = T::from(4.0).unwrap();
135    // Apply correction factor
136    let n1 = T::from(group1_indices.len()).unwrap();
137    let n2 = T::from(group2_indices.len()).unwrap();
138    let n = n1 + n2;
139
140    // Correction factor J
141    let j = one - three / (four * (n - two) - one);
142
143    // Calculate Hedge's g
144    let g = j * d;
145
146    Ok(g)
147}
148
149
single_statistics/testing/effect/mod.rs

single_statistics/testing/effect/
mod.rs
