single_statistics/testing/effect/
mod.rs

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

single_statistics/testing/effect/
mod.rs
