scirs2_stats/
multiple_testing.rs

1//! Multiple comparison correction methods
2//!
3//! When performing many simultaneous hypothesis tests, the probability of at least one
4//! false positive (Type I error) increases. Multiple testing corrections adjust p-values
5//! or significance thresholds to control for this.
6//!
7//! ## Methods Provided
8//!
9//! - **Bonferroni** - Controls the family-wise error rate (FWER) by multiplying each p-value by m
10//! - **Holm-Bonferroni** (step-down) - A more powerful step-down variant of Bonferroni
11//! - **Hochberg** (step-up) - A step-up variant that is more powerful than Holm under independence
12//! - **Benjamini-Hochberg** - Controls the false discovery rate (FDR)
13//! - **Benjamini-Yekutieli** - Controls FDR under arbitrary dependence
14//! - **Sidak** - Uses the Sidak correction: 1 - (1 - alpha)^(1/m)
15
16use crate::error::{StatsError, StatsResult};
17
18/// Result of a multiple testing correction
19#[derive(Debug, Clone)]
20pub struct MultipleCorrectionResult {
21    /// Adjusted p-values (in the same order as the input)
22    pub pvalues_corrected: Vec<f64>,
23    /// Boolean array indicating which hypotheses are rejected at the given alpha
24    pub reject: Vec<bool>,
25    /// The method used for correction
26    pub method: String,
27    /// The alpha level used
28    pub alpha: f64,
29}
30
31// ========================================================================
32// Bonferroni correction
33// ========================================================================
34
35/// Applies Bonferroni correction to a set of p-values.
36///
37/// The Bonferroni correction is the simplest and most conservative method for
38/// controlling the family-wise error rate (FWER). Each p-value is multiplied
39/// by the number of tests.
40///
41/// # Arguments
42///
43/// * `pvalues` - Slice of raw p-values
44/// * `alpha` - Significance level (e.g., 0.05)
45///
46/// # Returns
47///
48/// A `MultipleCorrectionResult` with adjusted p-values and rejection decisions.
49///
50/// # Examples
51///
52/// ```
53/// use scirs2_stats::multiple_testing::bonferroni;
54///
55/// let pvals = vec![0.01, 0.04, 0.03, 0.005];
56/// let result = bonferroni(&pvals, 0.05).expect("Bonferroni correction failed");
57/// // p-values are multiplied by 4
58/// assert!((result.pvalues_corrected[0] - 0.04).abs() < 1e-10);
59/// assert!((result.pvalues_corrected[3] - 0.02).abs() < 1e-10);
60/// ```
61pub fn bonferroni(pvalues: &[f64], alpha: f64) -> StatsResult<MultipleCorrectionResult> {
62    validate_inputs(pvalues, alpha)?;
63
64    let m = pvalues.len() as f64;
65    let corrected: Vec<f64> = pvalues.iter().map(|&p| (p * m).min(1.0)).collect();
66    let reject: Vec<bool> = corrected.iter().map(|&p| p <= alpha).collect();
67
68    Ok(MultipleCorrectionResult {
69        pvalues_corrected: corrected,
70        reject,
71        method: "Bonferroni".to_string(),
72        alpha,
73    })
74}
75
76// ========================================================================
77// Holm-Bonferroni (step-down)
78// ========================================================================
79
80/// Applies the Holm-Bonferroni step-down correction.
81///
82/// The Holm method is uniformly more powerful than Bonferroni while still
83/// controlling the FWER. It works by sorting p-values, then multiplying
84/// the i-th smallest p-value by (m - i + 1).
85///
86/// # Arguments
87///
88/// * `pvalues` - Slice of raw p-values
89/// * `alpha` - Significance level
90///
91/// # Returns
92///
93/// A `MultipleCorrectionResult` with adjusted p-values and rejection decisions.
94///
95/// # Examples
96///
97/// ```
98/// use scirs2_stats::multiple_testing::holm_bonferroni;
99///
100/// let pvals = vec![0.01, 0.04, 0.03, 0.005];
101/// let result = holm_bonferroni(&pvals, 0.05).expect("Holm correction failed");
102/// // The smallest p-value (0.005) is multiplied by 4, second (0.01) by 3, etc.
103/// assert!(result.reject[3]); // 0.005 * 4 = 0.02 <= 0.05
104/// ```
105pub fn holm_bonferroni(pvalues: &[f64], alpha: f64) -> StatsResult<MultipleCorrectionResult> {
106    validate_inputs(pvalues, alpha)?;
107
108    let m = pvalues.len();
109
110    // Sort indices by p-value
111    let mut order: Vec<usize> = (0..m).collect();
112    order.sort_by(|&a, &b| {
113        pvalues[a]
114            .partial_cmp(&pvalues[b])
115            .unwrap_or(std::cmp::Ordering::Equal)
116    });
117
118    let mut corrected = vec![0.0_f64; m];
119
120    // Step-down: enforce monotonicity (adjusted values must be non-decreasing)
121    let mut max_so_far = 0.0_f64;
122    for (rank, &idx) in order.iter().enumerate() {
123        let multiplier = (m - rank) as f64;
124        let adjusted = (pvalues[idx] * multiplier).min(1.0);
125        max_so_far = max_so_far.max(adjusted);
126        corrected[idx] = max_so_far;
127    }
128
129    let reject: Vec<bool> = corrected.iter().map(|&p| p <= alpha).collect();
130
131    Ok(MultipleCorrectionResult {
132        pvalues_corrected: corrected,
133        reject,
134        method: "Holm-Bonferroni".to_string(),
135        alpha,
136    })
137}
138
139// ========================================================================
140// Hochberg (step-up)
141// ========================================================================
142
143/// Applies the Hochberg step-up correction.
144///
145/// The Hochberg procedure is a step-up version of the Holm method. It is valid
146/// under certain positive dependence conditions (PRDS). The i-th largest p-value
147/// is multiplied by (m - i + 1), enforcing monotonicity from the top.
148///
149/// # Arguments
150///
151/// * `pvalues` - Slice of raw p-values
152/// * `alpha` - Significance level
153///
154/// # Returns
155///
156/// A `MultipleCorrectionResult` with adjusted p-values and rejection decisions.
157///
158/// # Examples
159///
160/// ```
161/// use scirs2_stats::multiple_testing::hochberg;
162///
163/// let pvals = vec![0.01, 0.04, 0.03, 0.005];
164/// let result = hochberg(&pvals, 0.05).expect("Hochberg correction failed");
165/// assert!(result.reject[3]); // smallest p-value should be rejected
166/// ```
167pub fn hochberg(pvalues: &[f64], alpha: f64) -> StatsResult<MultipleCorrectionResult> {
168    validate_inputs(pvalues, alpha)?;
169
170    let m = pvalues.len();
171
172    // Sort indices by p-value (ascending)
173    let mut order: Vec<usize> = (0..m).collect();
174    order.sort_by(|&a, &b| {
175        pvalues[a]
176            .partial_cmp(&pvalues[b])
177            .unwrap_or(std::cmp::Ordering::Equal)
178    });
179
180    let mut corrected = vec![0.0_f64; m];
181
182    // Step-up: go from largest to smallest, enforce monotonicity downward
183    let mut min_so_far = 1.0_f64;
184    for rank in (0..m).rev() {
185        let idx = order[rank];
186        let multiplier = (m - rank) as f64;
187        let adjusted = (pvalues[idx] * multiplier).min(1.0);
188        min_so_far = min_so_far.min(adjusted);
189        corrected[idx] = min_so_far;
190    }
191
192    let reject: Vec<bool> = corrected.iter().map(|&p| p <= alpha).collect();
193
194    Ok(MultipleCorrectionResult {
195        pvalues_corrected: corrected,
196        reject,
197        method: "Hochberg".to_string(),
198        alpha,
199    })
200}
201
202// ========================================================================
203// Benjamini-Hochberg (FDR control)
204// ========================================================================
205
206/// Applies the Benjamini-Hochberg procedure for FDR control.
207///
208/// The Benjamini-Hochberg (BH) procedure controls the false discovery rate (FDR)
209/// rather than the family-wise error rate. FDR is the expected proportion of false
210/// discoveries among all discoveries. The BH procedure is less conservative than
211/// FWER-controlling methods and is widely used in genomics, neuroimaging, etc.
212///
213/// # Arguments
214///
215/// * `pvalues` - Slice of raw p-values
216/// * `alpha` - Target FDR level (e.g., 0.05)
217///
218/// # Returns
219///
220/// A `MultipleCorrectionResult` with adjusted p-values (q-values) and rejection decisions.
221///
222/// # Examples
223///
224/// ```
225/// use scirs2_stats::multiple_testing::benjamini_hochberg;
226///
227/// let pvals = vec![0.001, 0.008, 0.039, 0.041, 0.042, 0.06, 0.074, 0.205, 0.212, 0.216];
228/// let result = benjamini_hochberg(&pvals, 0.05).expect("BH correction failed");
229/// // The first few p-values should still be significant after BH correction
230/// assert!(result.reject[0]);
231/// ```
232pub fn benjamini_hochberg(pvalues: &[f64], alpha: f64) -> StatsResult<MultipleCorrectionResult> {
233    validate_inputs(pvalues, alpha)?;
234
235    let m = pvalues.len();
236    let mf = m as f64;
237
238    // Sort indices by p-value (ascending)
239    let mut order: Vec<usize> = (0..m).collect();
240    order.sort_by(|&a, &b| {
241        pvalues[a]
242            .partial_cmp(&pvalues[b])
243            .unwrap_or(std::cmp::Ordering::Equal)
244    });
245
246    let mut corrected = vec![0.0_f64; m];
247
248    // Step-up: start from the largest p-value
249    let mut min_so_far = 1.0_f64;
250    for rank in (0..m).rev() {
251        let idx = order[rank];
252        let rank_1based = (rank + 1) as f64;
253        let adjusted = (pvalues[idx] * mf / rank_1based).min(1.0);
254        min_so_far = min_so_far.min(adjusted);
255        corrected[idx] = min_so_far;
256    }
257
258    let reject: Vec<bool> = corrected.iter().map(|&p| p <= alpha).collect();
259
260    Ok(MultipleCorrectionResult {
261        pvalues_corrected: corrected,
262        reject,
263        method: "Benjamini-Hochberg".to_string(),
264        alpha,
265    })
266}
267
268// ========================================================================
269// Benjamini-Yekutieli
270// ========================================================================
271
272/// Applies the Benjamini-Yekutieli procedure for FDR control under arbitrary dependence.
273///
274/// The BY procedure is a modification of BH that controls the FDR under arbitrary
275/// dependence between the test statistics. It is more conservative than BH because
276/// it divides by an additional factor c(m) = sum(1/i, i=1..m).
277///
278/// # Arguments
279///
280/// * `pvalues` - Slice of raw p-values
281/// * `alpha` - Target FDR level
282///
283/// # Returns
284///
285/// A `MultipleCorrectionResult` with adjusted p-values and rejection decisions.
286///
287/// # Examples
288///
289/// ```
290/// use scirs2_stats::multiple_testing::benjamini_yekutieli;
291///
292/// let pvals = vec![0.001, 0.008, 0.039, 0.041, 0.06];
293/// let result = benjamini_yekutieli(&pvals, 0.05).expect("BY correction failed");
294/// assert!(result.reject[0]); // Smallest p-value should still be significant
295/// ```
296pub fn benjamini_yekutieli(pvalues: &[f64], alpha: f64) -> StatsResult<MultipleCorrectionResult> {
297    validate_inputs(pvalues, alpha)?;
298
299    let m = pvalues.len();
300    let mf = m as f64;
301
302    // c(m) = sum(1/i for i=1..m)
303    let cm: f64 = (1..=m).map(|i| 1.0 / i as f64).sum();
304
305    // Sort indices by p-value (ascending)
306    let mut order: Vec<usize> = (0..m).collect();
307    order.sort_by(|&a, &b| {
308        pvalues[a]
309            .partial_cmp(&pvalues[b])
310            .unwrap_or(std::cmp::Ordering::Equal)
311    });
312
313    let mut corrected = vec![0.0_f64; m];
314
315    // Step-up with c(m) correction
316    let mut min_so_far = 1.0_f64;
317    for rank in (0..m).rev() {
318        let idx = order[rank];
319        let rank_1based = (rank + 1) as f64;
320        let adjusted = (pvalues[idx] * mf * cm / rank_1based).min(1.0);
321        min_so_far = min_so_far.min(adjusted);
322        corrected[idx] = min_so_far;
323    }
324
325    let reject: Vec<bool> = corrected.iter().map(|&p| p <= alpha).collect();
326
327    Ok(MultipleCorrectionResult {
328        pvalues_corrected: corrected,
329        reject,
330        method: "Benjamini-Yekutieli".to_string(),
331        alpha,
332    })
333}
334
335// ========================================================================
336// Sidak correction
337// ========================================================================
338
339/// Applies the Sidak correction.
340///
341/// The Sidak correction is based on the formula: adjusted_p = 1 - (1 - p)^m.
342/// For independent tests, it is slightly less conservative than Bonferroni.
343///
344/// # Arguments
345///
346/// * `pvalues` - Slice of raw p-values
347/// * `alpha` - Significance level
348///
349/// # Returns
350///
351/// A `MultipleCorrectionResult` with adjusted p-values and rejection decisions.
352///
353/// # Examples
354///
355/// ```
356/// use scirs2_stats::multiple_testing::sidak;
357///
358/// let pvals = vec![0.01, 0.04, 0.03, 0.005];
359/// let result = sidak(&pvals, 0.05).expect("Sidak correction failed");
360/// // Sidak is slightly less conservative than Bonferroni
361/// assert!(result.pvalues_corrected[3] < 0.02); // 1 - (1 - 0.005)^4
362/// ```
363pub fn sidak(pvalues: &[f64], alpha: f64) -> StatsResult<MultipleCorrectionResult> {
364    validate_inputs(pvalues, alpha)?;
365
366    let m = pvalues.len() as f64;
367    let corrected: Vec<f64> = pvalues
368        .iter()
369        .map(|&p| {
370            // 1 - (1 - p)^m, but handle edge cases
371            if p >= 1.0 {
372                1.0
373            } else if p <= 0.0 {
374                0.0
375            } else {
376                (1.0 - (1.0 - p).powf(m)).min(1.0)
377            }
378        })
379        .collect();
380
381    let reject: Vec<bool> = corrected.iter().map(|&p| p <= alpha).collect();
382
383    Ok(MultipleCorrectionResult {
384        pvalues_corrected: corrected,
385        reject,
386        method: "Sidak".to_string(),
387        alpha,
388    })
389}
390
391// ========================================================================
392// Convenience: apply correction by name
393// ========================================================================
394
395/// Applies a multiple testing correction by method name.
396///
397/// This is a convenience function that dispatches to the appropriate correction
398/// method based on the given name string.
399///
400/// # Arguments
401///
402/// * `pvalues` - Slice of raw p-values
403/// * `alpha` - Significance level
404/// * `method` - Method name: "bonferroni", "holm", "hochberg", "fdr_bh", "fdr_by", "sidak"
405///
406/// # Returns
407///
408/// A `MultipleCorrectionResult` with adjusted p-values and rejection decisions.
409pub fn multipletests(
410    pvalues: &[f64],
411    alpha: f64,
412    method: &str,
413) -> StatsResult<MultipleCorrectionResult> {
414    match method {
415        "bonferroni" => bonferroni(pvalues, alpha),
416        "holm" | "holm-bonferroni" => holm_bonferroni(pvalues, alpha),
417        "hochberg" => hochberg(pvalues, alpha),
418        "fdr_bh" | "benjamini-hochberg" | "bh" => benjamini_hochberg(pvalues, alpha),
419        "fdr_by" | "benjamini-yekutieli" | "by" => benjamini_yekutieli(pvalues, alpha),
420        "sidak" => sidak(pvalues, alpha),
421        _ => Err(StatsError::InvalidArgument(format!(
422            "Unknown correction method '{}'. Valid methods: bonferroni, holm, hochberg, fdr_bh, fdr_by, sidak",
423            method
424        ))),
425    }
426}
427
428// ========================================================================
429// Input validation helper
430// ========================================================================
431
432fn validate_inputs(pvalues: &[f64], alpha: f64) -> StatsResult<()> {
433    if pvalues.is_empty() {
434        return Err(StatsError::InvalidArgument(
435            "p-values array cannot be empty".to_string(),
436        ));
437    }
438    if alpha <= 0.0 || alpha >= 1.0 {
439        return Err(StatsError::InvalidArgument(format!(
440            "alpha must be in (0, 1), got {}",
441            alpha
442        )));
443    }
444    for (i, &p) in pvalues.iter().enumerate() {
445        if p < 0.0 || p > 1.0 {
446            return Err(StatsError::InvalidArgument(format!(
447                "p-value at index {} is {}, must be in [0, 1]",
448                i, p
449            )));
450        }
451    }
452    Ok(())
453}
454
455// ========================================================================
456// Tests
457// ========================================================================
458
459#[cfg(test)]
460mod tests {
461    use super::*;
462
463    #[test]
464    fn test_bonferroni_basic() {
465        let pvals = vec![0.01, 0.04, 0.03, 0.005];
466        let result = bonferroni(&pvals, 0.05).expect("Bonferroni failed");
467
468        assert_eq!(result.pvalues_corrected.len(), 4);
469        assert!((result.pvalues_corrected[0] - 0.04).abs() < 1e-10);
470        assert!((result.pvalues_corrected[1] - 0.16).abs() < 1e-10);
471        assert!((result.pvalues_corrected[2] - 0.12).abs() < 1e-10);
472        assert!((result.pvalues_corrected[3] - 0.02).abs() < 1e-10);
473
474        assert!(result.reject[0]); // 0.04 <= 0.05
475        assert!(!result.reject[1]); // 0.16 > 0.05
476        assert!(!result.reject[2]); // 0.12 > 0.05
477        assert!(result.reject[3]); // 0.02 <= 0.05
478    }
479
480    #[test]
481    fn test_bonferroni_clamped() {
482        let pvals = vec![0.5, 0.6, 0.7];
483        let result = bonferroni(&pvals, 0.05).expect("Bonferroni failed");
484        // 0.5 * 3 = 1.5 => clamped to 1.0
485        assert!((result.pvalues_corrected[0] - 1.0).abs() < 1e-10);
486    }
487
488    #[test]
489    fn test_holm_bonferroni_basic() {
490        let pvals = vec![0.01, 0.04, 0.03, 0.005];
491        let result = holm_bonferroni(&pvals, 0.05).expect("Holm failed");
492
493        assert_eq!(result.pvalues_corrected.len(), 4);
494        // Sorted order: [0.005(3), 0.01(0), 0.03(2), 0.04(1)]
495        // Rank 0: 0.005 * 4 = 0.02
496        // Rank 1: 0.01 * 3 = 0.03, max(0.02, 0.03) = 0.03
497        // Rank 2: 0.03 * 2 = 0.06, max(0.03, 0.06) = 0.06
498        // Rank 3: 0.04 * 1 = 0.04, max(0.06, 0.04) = 0.06
499        assert!((result.pvalues_corrected[3] - 0.02).abs() < 1e-10);
500        assert!((result.pvalues_corrected[0] - 0.03).abs() < 1e-10);
501        assert!(result.reject[3]); // 0.02 <= 0.05
502        assert!(result.reject[0]); // 0.03 <= 0.05
503    }
504
505    #[test]
506    fn test_hochberg_basic() {
507        let pvals = vec![0.01, 0.04, 0.03, 0.005];
508        let result = hochberg(&pvals, 0.05).expect("Hochberg failed");
509
510        assert_eq!(result.pvalues_corrected.len(), 4);
511        // Should be at least as powerful as Bonferroni
512        let bonf_result = bonferroni(&pvals, 0.05).expect("Bonferroni failed");
513        for i in 0..pvals.len() {
514            assert!(result.pvalues_corrected[i] <= bonf_result.pvalues_corrected[i] + 1e-10);
515        }
516    }
517
518    #[test]
519    fn test_benjamini_hochberg_basic() {
520        let pvals = vec![
521            0.001, 0.008, 0.039, 0.041, 0.042, 0.06, 0.074, 0.205, 0.212, 0.216,
522        ];
523        let result = benjamini_hochberg(&pvals, 0.05).expect("BH failed");
524
525        assert_eq!(result.pvalues_corrected.len(), 10);
526        // All corrected p-values should be in [0, 1]
527        for &p in &result.pvalues_corrected {
528            assert!(p >= 0.0 && p <= 1.0);
529        }
530        // The first p-value should still be significant
531        assert!(result.reject[0]);
532    }
533
534    #[test]
535    fn test_benjamini_hochberg_monotonicity() {
536        let pvals = vec![0.01, 0.02, 0.03, 0.04, 0.05];
537        let result = benjamini_hochberg(&pvals, 0.05).expect("BH failed");
538
539        // Adjusted p-values should be non-decreasing when sorted by original p-value
540        let mut order: Vec<usize> = (0..pvals.len()).collect();
541        order.sort_by(|&a, &b| {
542            pvals[a]
543                .partial_cmp(&pvals[b])
544                .unwrap_or(std::cmp::Ordering::Equal)
545        });
546
547        for i in 1..order.len() {
548            assert!(
549                result.pvalues_corrected[order[i]]
550                    >= result.pvalues_corrected[order[i - 1]] - 1e-10,
551                "Monotonicity violation at rank {}",
552                i
553            );
554        }
555    }
556
557    #[test]
558    fn test_benjamini_yekutieli_basic() {
559        let pvals = vec![0.001, 0.008, 0.039, 0.041, 0.06];
560        let result = benjamini_yekutieli(&pvals, 0.05).expect("BY failed");
561
562        // BY should be more conservative than BH
563        let bh_result = benjamini_hochberg(&pvals, 0.05).expect("BH failed");
564        for i in 0..pvals.len() {
565            assert!(result.pvalues_corrected[i] >= bh_result.pvalues_corrected[i] - 1e-10);
566        }
567    }
568
569    #[test]
570    fn test_sidak_basic() {
571        let pvals = vec![0.01, 0.04, 0.03, 0.005];
572        let result = sidak(&pvals, 0.05).expect("Sidak failed");
573
574        // Sidak should be slightly less conservative than Bonferroni
575        let bonf_result = bonferroni(&pvals, 0.05).expect("Bonferroni failed");
576        for i in 0..pvals.len() {
577            assert!(result.pvalues_corrected[i] <= bonf_result.pvalues_corrected[i] + 1e-10);
578        }
579    }
580
581    #[test]
582    fn test_sidak_values() {
583        let pvals = vec![0.005];
584        let result = sidak(&pvals, 0.05).expect("Sidak failed");
585        // For m=1, Sidak(p) = 1 - (1-p)^1 = p
586        assert!((result.pvalues_corrected[0] - 0.005).abs() < 1e-10);
587
588        let pvals2 = vec![0.01, 0.01];
589        let result2 = sidak(&pvals2, 0.05).expect("Sidak failed");
590        // 1 - (1 - 0.01)^2 = 1 - 0.99^2 = 1 - 0.9801 = 0.0199
591        assert!((result2.pvalues_corrected[0] - 0.0199).abs() < 1e-4);
592    }
593
594    #[test]
595    fn test_multipletests_dispatch() {
596        let pvals = vec![0.01, 0.04, 0.03, 0.005];
597
598        let r1 = multipletests(&pvals, 0.05, "bonferroni").expect("dispatch failed");
599        assert_eq!(r1.method, "Bonferroni");
600
601        let r2 = multipletests(&pvals, 0.05, "holm").expect("dispatch failed");
602        assert_eq!(r2.method, "Holm-Bonferroni");
603
604        let r3 = multipletests(&pvals, 0.05, "hochberg").expect("dispatch failed");
605        assert_eq!(r3.method, "Hochberg");
606
607        let r4 = multipletests(&pvals, 0.05, "fdr_bh").expect("dispatch failed");
608        assert_eq!(r4.method, "Benjamini-Hochberg");
609
610        let r5 = multipletests(&pvals, 0.05, "fdr_by").expect("dispatch failed");
611        assert_eq!(r5.method, "Benjamini-Yekutieli");
612
613        let r6 = multipletests(&pvals, 0.05, "sidak").expect("dispatch failed");
614        assert_eq!(r6.method, "Sidak");
615
616        let r7 = multipletests(&pvals, 0.05, "unknown_method");
617        assert!(r7.is_err());
618    }
619
620    #[test]
621    fn test_empty_input() {
622        let empty: Vec<f64> = vec![];
623        assert!(bonferroni(&empty, 0.05).is_err());
624        assert!(holm_bonferroni(&empty, 0.05).is_err());
625        assert!(hochberg(&empty, 0.05).is_err());
626        assert!(benjamini_hochberg(&empty, 0.05).is_err());
627        assert!(benjamini_yekutieli(&empty, 0.05).is_err());
628        assert!(sidak(&empty, 0.05).is_err());
629    }
630
631    #[test]
632    fn test_invalid_alpha() {
633        let pvals = vec![0.01, 0.05];
634        assert!(bonferroni(&pvals, 0.0).is_err());
635        assert!(bonferroni(&pvals, 1.0).is_err());
636        assert!(bonferroni(&pvals, -0.1).is_err());
637    }
638
639    #[test]
640    fn test_invalid_pvalues() {
641        let pvals = vec![0.01, 1.5];
642        assert!(bonferroni(&pvals, 0.05).is_err());
643
644        let pvals2 = vec![-0.01, 0.05];
645        assert!(bonferroni(&pvals2, 0.05).is_err());
646    }
647
648    #[test]
649    fn test_single_pvalue() {
650        let pvals = vec![0.03];
651        let result = bonferroni(&pvals, 0.05).expect("single pval failed");
652        assert!((result.pvalues_corrected[0] - 0.03).abs() < 1e-10);
653        assert!(result.reject[0]);
654    }
655
656    #[test]
657    fn test_all_significant() {
658        let pvals = vec![0.001, 0.002, 0.003];
659        let result = bonferroni(&pvals, 0.05).expect("all sig failed");
660        // 0.001 * 3 = 0.003, 0.002 * 3 = 0.006, 0.003 * 3 = 0.009
661        assert!(result.reject.iter().all(|&r| r));
662    }
663
664    #[test]
665    fn test_none_significant() {
666        let pvals = vec![0.5, 0.6, 0.7, 0.8, 0.9];
667        let result = bonferroni(&pvals, 0.05).expect("none sig failed");
668        assert!(result.reject.iter().all(|&r| !r));
669    }
670
671    #[test]
672    fn test_holm_vs_bonferroni_power() {
673        // Holm should reject at least as many as Bonferroni
674        let pvals = vec![0.001, 0.01, 0.04, 0.06, 0.10];
675        let bonf = bonferroni(&pvals, 0.05).expect("Bonferroni failed");
676        let holm = holm_bonferroni(&pvals, 0.05).expect("Holm failed");
677
678        let bonf_count: usize = bonf.reject.iter().filter(|&&r| r).count();
679        let holm_count: usize = holm.reject.iter().filter(|&&r| r).count();
680        assert!(holm_count >= bonf_count);
681    }
682
683    #[test]
684    fn test_bh_less_conservative_than_bonferroni() {
685        let pvals = vec![0.001, 0.01, 0.04, 0.06, 0.10];
686        let bonf = bonferroni(&pvals, 0.05).expect("Bonferroni failed");
687        let bh = benjamini_hochberg(&pvals, 0.05).expect("BH failed");
688
689        let bonf_count: usize = bonf.reject.iter().filter(|&&r| r).count();
690        let bh_count: usize = bh.reject.iter().filter(|&&r| r).count();
691        // BH (FDR) should reject at least as many as Bonferroni (FWER)
692        assert!(bh_count >= bonf_count);
693    }
694}
scirs2_stats/multiple_testing.rs

scirs2_stats/
multiple_testing.rs
