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
118
119
120
121
122
123
124
125
126
127
//! Tukey's Honestly Significant Difference (HSD) post-hoc test.
//!
//! Mirrors the reference `pairwise_tukeyhsd`. For all pairs of groups it
//! reports the mean differences, the studentized-range based confidence
//! intervals, the adjusted p-values, and the reject decisions at family-wise
//! error rate `alpha`. The critical value and p-values come from the
//! studentized-range distribution (see [`crate::srange`]).
use crate::srange::{srange_ppf, srange_sf};
/// Result of a pairwise Tukey HSD comparison.
#[derive(Debug, Clone)]
pub struct TukeyHsdResult {
/// `(i, j)` group-index pairs (upper triangle, `i < j`).
pub pairs: Vec<(usize, usize)>,
/// Mean difference `mean_j − mean_i` for each pair.
pub meandiffs: Vec<f64>,
/// Standard error of each pairwise difference (studentized-range scaling).
pub std_pairs: Vec<f64>,
/// Lower and upper confidence bounds for each mean difference.
pub confint: Vec<(f64, f64)>,
/// Studentized-range adjusted p-value for each pair.
pub pvalues: Vec<f64>,
/// Reject the null of equal means at level `alpha`.
pub reject: Vec<bool>,
/// Critical value of the studentized range at level `alpha`.
pub q_crit: f64,
/// Total residual degrees of freedom `N − k`.
pub df_total: f64,
/// Pooled within-group variance (mean squared error).
pub variance: f64,
}
/// All pairwise Tukey HSD comparisons.
///
/// `data` holds the response values and `groups` the integer group label of
/// each observation (labels need not be contiguous; distinct labels define the
/// groups, sorted ascending). `alpha` is the family-wise error rate.
pub fn pairwise_tukeyhsd(data: &[f64], groups: &[usize], alpha: f64) -> TukeyHsdResult {
assert_eq!(data.len(), groups.len(), "data and groups length mismatch");
// Collect the sorted unique group labels.
let mut labels: Vec<usize> = groups.to_vec();
labels.sort_unstable();
labels.dedup();
let k = labels.len();
// Per-group sums, counts.
let mut nobs = vec![0.0_f64; k];
let mut sums = vec![0.0_f64; k];
let label_idx = |g: usize| labels.iter().position(|&l| l == g).unwrap();
for (&v, &g) in data.iter().zip(groups) {
let gi = label_idx(g);
nobs[gi] += 1.0;
sums[gi] += v;
}
let means: Vec<f64> = sums.iter().zip(&nobs).map(|(&s, &n)| s / n).collect();
// Pooled within-group variance: SS_within / (N − k).
let total_n: f64 = nobs.iter().sum();
let mut ss_within = 0.0;
for (&v, &g) in data.iter().zip(groups) {
let gi = label_idx(g);
let d = v - means[gi];
ss_within += d * d;
}
let df_total = total_n - k as f64;
let variance = ss_within / df_total;
// All upper-triangle pairs.
let mut pairs = Vec::new();
let mut meandiffs = Vec::new();
let mut std_pairs = Vec::new();
for i in 0..k {
for j in (i + 1)..k {
pairs.push((i, j));
// meandiffs use mean[j] − mean[i] (reference sign convention).
meandiffs.push(means[j] - means[i]);
// var of the studentized-range difference: var * (1/n_i + 1/n_j) / 2.
let var_pair = variance * (1.0 / nobs[i] + 1.0 / nobs[j]) / 2.0;
std_pairs.push(var_pair.sqrt());
}
}
let kf = k as f64;
let q_crit = srange_ppf(1.0 - alpha, kf, df_total);
let mut pvalues = Vec::with_capacity(pairs.len());
let mut reject = Vec::with_capacity(pairs.len());
let mut confint = Vec::with_capacity(pairs.len());
for idx in 0..pairs.len() {
let st_range = meandiffs[idx].abs() / std_pairs[idx];
let pval = srange_sf(st_range, kf, df_total);
pvalues.push(pval);
reject.push(st_range > q_crit);
let crit_int = std_pairs[idx] * q_crit;
confint.push((meandiffs[idx] - crit_int, meandiffs[idx] + crit_int));
}
TukeyHsdResult {
pairs,
meandiffs,
std_pairs,
confint,
pvalues,
reject,
q_crit,
df_total,
variance,
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn three_groups_pairs() {
let data = [1.0, 2.0, 3.0, 5.0, 6.0, 4.0, 8.0, 9.0, 7.0, 10.0];
let groups = [0, 0, 0, 1, 1, 1, 2, 2, 2, 2];
let res = pairwise_tukeyhsd(&data, &groups, 0.05);
assert_eq!(res.pairs.len(), 3);
assert_eq!(res.df_total, 7.0);
// Group means clearly separated -> all rejected.
assert!(res.reject.iter().all(|&r| r));
}
}