solow_stats/
tukey.rs

1//! Tukey's Honestly Significant Difference (HSD) post-hoc test.
2//!
3//! Mirrors the reference `pairwise_tukeyhsd`. For all pairs of groups it
4//! reports the mean differences, the studentized-range based confidence
5//! intervals, the adjusted p-values, and the reject decisions at family-wise
6//! error rate `alpha`. The critical value and p-values come from the
7//! studentized-range distribution (see [`crate::srange`]).
8
9use crate::srange::{srange_ppf, srange_sf};
10
11/// Result of a pairwise Tukey HSD comparison.
12#[derive(Debug, Clone)]
13pub struct TukeyHsdResult {
14    /// `(i, j)` group-index pairs (upper triangle, `i < j`).
15    pub pairs: Vec<(usize, usize)>,
16    /// Mean difference `mean_j − mean_i` for each pair.
17    pub meandiffs: Vec<f64>,
18    /// Standard error of each pairwise difference (studentized-range scaling).
19    pub std_pairs: Vec<f64>,
20    /// Lower and upper confidence bounds for each mean difference.
21    pub confint: Vec<(f64, f64)>,
22    /// Studentized-range adjusted p-value for each pair.
23    pub pvalues: Vec<f64>,
24    /// Reject the null of equal means at level `alpha`.
25    pub reject: Vec<bool>,
26    /// Critical value of the studentized range at level `alpha`.
27    pub q_crit: f64,
28    /// Total residual degrees of freedom `N − k`.
29    pub df_total: f64,
30    /// Pooled within-group variance (mean squared error).
31    pub variance: f64,
32}
33
34/// All pairwise Tukey HSD comparisons.
35///
36/// `data` holds the response values and `groups` the integer group label of
37/// each observation (labels need not be contiguous; distinct labels define the
38/// groups, sorted ascending). `alpha` is the family-wise error rate.
39pub fn pairwise_tukeyhsd(data: &[f64], groups: &[usize], alpha: f64) -> TukeyHsdResult {
40    assert_eq!(data.len(), groups.len(), "data and groups length mismatch");
41
42    // Collect the sorted unique group labels.
43    let mut labels: Vec<usize> = groups.to_vec();
44    labels.sort_unstable();
45    labels.dedup();
46    let k = labels.len();
47
48    // Per-group sums, counts.
49    let mut nobs = vec![0.0_f64; k];
50    let mut sums = vec![0.0_f64; k];
51    let label_idx = |g: usize| labels.iter().position(|&l| l == g).unwrap();
52    for (&v, &g) in data.iter().zip(groups) {
53        let gi = label_idx(g);
54        nobs[gi] += 1.0;
55        sums[gi] += v;
56    }
57    let means: Vec<f64> = sums.iter().zip(&nobs).map(|(&s, &n)| s / n).collect();
58
59    // Pooled within-group variance: SS_within / (N − k).
60    let total_n: f64 = nobs.iter().sum();
61    let mut ss_within = 0.0;
62    for (&v, &g) in data.iter().zip(groups) {
63        let gi = label_idx(g);
64        let d = v - means[gi];
65        ss_within += d * d;
66    }
67    let df_total = total_n - k as f64;
68    let variance = ss_within / df_total;
69
70    // All upper-triangle pairs.
71    let mut pairs = Vec::new();
72    let mut meandiffs = Vec::new();
73    let mut std_pairs = Vec::new();
74    for i in 0..k {
75        for j in (i + 1)..k {
76            pairs.push((i, j));
77            // meandiffs use mean[j] − mean[i] (reference sign convention).
78            meandiffs.push(means[j] - means[i]);
79            // var of the studentized-range difference: var * (1/n_i + 1/n_j) / 2.
80            let var_pair = variance * (1.0 / nobs[i] + 1.0 / nobs[j]) / 2.0;
81            std_pairs.push(var_pair.sqrt());
82        }
83    }
84
85    let kf = k as f64;
86    let q_crit = srange_ppf(1.0 - alpha, kf, df_total);
87
88    let mut pvalues = Vec::with_capacity(pairs.len());
89    let mut reject = Vec::with_capacity(pairs.len());
90    let mut confint = Vec::with_capacity(pairs.len());
91    for idx in 0..pairs.len() {
92        let st_range = meandiffs[idx].abs() / std_pairs[idx];
93        let pval = srange_sf(st_range, kf, df_total);
94        pvalues.push(pval);
95        reject.push(st_range > q_crit);
96        let crit_int = std_pairs[idx] * q_crit;
97        confint.push((meandiffs[idx] - crit_int, meandiffs[idx] + crit_int));
98    }
99
100    TukeyHsdResult {
101        pairs,
102        meandiffs,
103        std_pairs,
104        confint,
105        pvalues,
106        reject,
107        q_crit,
108        df_total,
109        variance,
110    }
111}
112
113#[cfg(test)]
114mod tests {
115    use super::*;
116
117    #[test]
118    fn three_groups_pairs() {
119        let data = [1.0, 2.0, 3.0, 5.0, 6.0, 4.0, 8.0, 9.0, 7.0, 10.0];
120        let groups = [0, 0, 0, 1, 1, 1, 2, 2, 2, 2];
121        let res = pairwise_tukeyhsd(&data, &groups, 0.05);
122        assert_eq!(res.pairs.len(), 3);
123        assert_eq!(res.df_total, 7.0);
124        // Group means clearly separated -> all rejected.
125        assert!(res.reject.iter().all(|&r| r));
126    }
127}
solow_stats/tukey.rs

solow_stats/
tukey.rs
