1use serde::{Deserialize, Serialize};
33
34use super::join::{Ratio, RatioMethod};
35
36pub const MIN_REPLICATES: usize = 5;
38
39#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
42#[serde(rename_all = "snake_case")]
43pub enum ArmOrder {
44 SubjectFirst,
46 ComparatorFirst,
48}
49
50#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
53#[serde(deny_unknown_fields)]
54pub struct ReplicatePair {
55 pub subject: f64,
57 pub comparator: f64,
59 pub order: ArmOrder,
61}
62
63#[must_use]
70pub fn t_lower_one_sided_95(df: usize) -> f64 {
71 const TABLE: [f64; 30] = [
72 6.314, 2.920, 2.353, 2.132, 2.015, 1.943, 1.895, 1.860, 1.833, 1.812, 1.796, 1.782, 1.771,
73 1.761, 1.753, 1.746, 1.740, 1.734, 1.729, 1.725, 1.721, 1.717, 1.714, 1.711, 1.708, 1.706,
74 1.703, 1.701, 1.699, 1.697,
75 ];
76 match df {
77 0 => f64::INFINITY,
78 d if d <= TABLE.len() => TABLE[d - 1],
79 _ => 1.645,
80 }
81}
82
83#[must_use]
95pub fn log_ratio_lcb(pairs: &[ReplicatePair]) -> Option<Ratio> {
96 if pairs.len() < MIN_REPLICATES || !is_strictly_alternating(pairs) {
97 return None;
98 }
99 let logs = log_ratios(pairs)?;
100 let n = logs.len();
101 let mean = logs.iter().sum::<f64>() / n as f64;
102 let var = logs.iter().map(|l| (l - mean).powi(2)).sum::<f64>() / (n as f64 - 1.0);
103 let se = (var / n as f64).sqrt();
104 let t = t_lower_one_sided_95(n - 1);
105 Some(Ratio {
106 point: mean.exp(),
107 lcb95: Some(t.mul_add(-se, mean).exp()),
108 method: RatioMethod::ReplicateTLower,
109 n,
110 })
111}
112
113#[must_use]
120pub fn log_ratio_point(pairs: &[ReplicatePair]) -> Option<Ratio> {
121 let logs = log_ratios(pairs)?;
122 let n = logs.len();
123 let mean = logs.iter().sum::<f64>() / n as f64;
124 Some(Ratio::reporting_only(
125 mean.exp(),
126 RatioMethod::ReplicateTLower,
127 n,
128 ))
129}
130
131#[must_use]
134pub fn log_ratio_bound_or_point(pairs: &[ReplicatePair]) -> Option<Ratio> {
135 log_ratio_lcb(pairs).or_else(|| log_ratio_point(pairs))
136}
137
138fn log_ratios(pairs: &[ReplicatePair]) -> Option<Vec<f64>> {
139 if pairs.len() < 2 {
140 return None;
141 }
142 pairs
143 .iter()
144 .map(|p| {
145 if p.subject > 0.0 && p.comparator > 0.0 {
146 Some((p.subject / p.comparator).ln())
147 } else {
148 None
149 }
150 })
151 .collect()
152}
153
154fn is_strictly_alternating(pairs: &[ReplicatePair]) -> bool {
156 pairs.windows(2).all(|w| w[0].order != w[1].order)
157}
158
159#[cfg(test)]
160mod tests {
161 use super::*;
162
163 fn alternating(values: &[(f64, f64)]) -> Vec<ReplicatePair> {
164 values
165 .iter()
166 .enumerate()
167 .map(|(i, &(subject, comparator))| ReplicatePair {
168 subject,
169 comparator,
170 order: if i % 2 == 0 {
171 ArmOrder::SubjectFirst
172 } else {
173 ArmOrder::ComparatorFirst
174 },
175 })
176 .collect()
177 }
178
179 #[test]
181 fn one_sided_t_table_matches_published_values() {
182 for (df, want) in [
183 (1_usize, 6.314_f64),
184 (2, 2.920),
185 (3, 2.353),
186 (4, 2.132),
187 (5, 2.015),
188 (6, 1.943),
189 (7, 1.895),
190 (8, 1.860),
191 (9, 1.833),
192 (10, 1.812),
193 (11, 1.796),
194 (12, 1.782),
195 (13, 1.771),
196 (14, 1.761),
197 (15, 1.753),
198 (16, 1.746),
199 (17, 1.740),
200 (18, 1.734),
201 (19, 1.729),
202 (20, 1.725),
203 (21, 1.721),
204 (22, 1.717),
205 (23, 1.714),
206 (24, 1.711),
207 (25, 1.708),
208 (26, 1.706),
209 (27, 1.703),
210 (28, 1.701),
211 (29, 1.699),
212 (30, 1.697),
213 ] {
214 assert_eq!(t_lower_one_sided_95(df), want, "df={df}");
215 }
216 assert_eq!(t_lower_one_sided_95(31), 1.645, "beyond 30, the normal");
217 assert_eq!(t_lower_one_sided_95(1_000), 1.645);
218 assert!(
219 t_lower_one_sided_95(0).is_infinite(),
220 "df=0 supports no bound"
221 );
222 }
223
224 #[test]
227 fn the_table_is_one_sided_not_two_tailed() {
228 assert_eq!(t_lower_one_sided_95(4), 2.132);
229 assert_ne!(t_lower_one_sided_95(4), 2.776);
230 }
231
232 #[test]
234 fn fewer_than_five_replicates_give_no_bound() {
235 let three = alternating(&[(100.0, 90.0), (101.0, 91.0), (99.0, 89.0)]);
236 assert!(log_ratio_lcb(&three).is_none());
237 assert_eq!(MIN_REPLICATES, 5);
238
239 let reporting = log_ratio_point(&three).expect("point estimate exists");
241 assert!(reporting.lcb95.is_none());
242 assert_eq!(reporting.n, 3);
243 assert!(!reporting.passes(0.0), "no bound is not a pass");
244
245 let five = alternating(&[
247 (100.0, 90.0),
248 (101.0, 91.0),
249 (99.0, 89.0),
250 (100.5, 90.5),
251 (100.2, 90.1),
252 ]);
253 assert!(log_ratio_lcb(&five).is_some());
254 }
255
256 #[test]
259 fn non_alternating_order_is_refused() {
260 let mut pairs = alternating(&[
261 (100.0, 90.0),
262 (101.0, 91.0),
263 (99.0, 89.0),
264 (100.5, 90.5),
265 (100.2, 90.1),
266 ]);
267 assert!(log_ratio_lcb(&pairs).is_some(), "control: alternating");
268 pairs[3].order = pairs[2].order;
269 assert!(
270 log_ratio_lcb(&pairs).is_none(),
271 "two consecutive replicates led with the same arm"
272 );
273 assert!(log_ratio_point(&pairs).is_some());
275 }
276
277 #[test]
281 fn log_ratio_bound_is_exponentiated() {
282 let flat = alternating(&[
285 (110.0, 100.0),
286 (220.0, 200.0),
287 (55.0, 50.0),
288 (11.0, 10.0),
289 (1100.0, 1000.0),
290 ]);
291 let r = log_ratio_lcb(&flat).expect("n = 5, alternating");
292 assert!((r.point - 1.10).abs() < 1e-12, "{r:?}");
293 assert!(
294 (r.lcb95.expect("bounded") - 1.10).abs() < 1e-12,
295 "zero variance leaves the bound at the point: {r:?}"
296 );
297 assert_eq!(r.method, RatioMethod::ReplicateTLower);
298 assert_eq!(r.n, 5);
299
300 let skewed = alternating(&[
303 (50.0, 100.0),
304 (200.0, 100.0),
305 (50.0, 100.0),
306 (200.0, 100.0),
307 (100.0, 100.0),
308 ]);
309 let g = log_ratio_lcb(&skewed).expect("n = 5");
310 assert!((g.point - 1.0).abs() < 1e-12, "geometric mean: {g:?}");
311 assert!(g.lcb95.expect("bounded") < g.point, "{g:?}");
312 }
313
314 #[test]
316 fn more_dispersion_lowers_the_bound() {
317 let tight = alternating(&[
318 (110.0, 100.0),
319 (109.0, 100.0),
320 (111.0, 100.0),
321 (110.5, 100.0),
322 (109.5, 100.0),
323 ]);
324 let loose = alternating(&[
325 (60.0, 100.0),
326 (160.0, 100.0),
327 (70.0, 100.0),
328 (150.0, 100.0),
329 (110.0, 100.0),
330 ]);
331 let a = log_ratio_lcb(&tight).expect("n = 5");
332 let b = log_ratio_lcb(&loose).expect("n = 5");
333 assert!(
334 b.lcb95.expect("bounded") < a.lcb95.expect("bounded"),
335 "dispersed {b:?} must bound lower than tight {a:?}"
336 );
337 }
338
339 #[test]
343 fn a_single_replicate_has_no_log_ratio() {
344 let one = alternating(&[(110.0, 100.0)]);
345 assert!(log_ratio_point(&one).is_none());
346 assert!(log_ratio_lcb(&one).is_none());
347 assert!(log_ratio_bound_or_point(&one).is_none());
348 assert!(log_ratio_point(&[]).is_none());
349 let two = alternating(&[(110.0, 100.0), (90.0, 100.0)]);
351 let r = log_ratio_point(&two).expect("two pairs give a point");
352 assert_eq!(r.n, 2);
353 assert!(r.lcb95.is_none());
354 }
355
356 #[test]
359 fn a_zero_lane_has_no_log_ratio() {
360 let zeroed = alternating(&[
361 (110.0, 100.0),
362 (0.0, 100.0),
363 (111.0, 100.0),
364 (110.5, 100.0),
365 (109.5, 100.0),
366 ]);
367 assert!(log_ratio_lcb(&zeroed).is_none());
368 assert!(log_ratio_point(&zeroed).is_none());
369 }
370
371 #[test]
375 fn the_wrapper_falls_back_to_reporting_only() {
376 let three = alternating(&[(100.0, 90.0), (101.0, 91.0), (99.0, 89.0)]);
377 let r = log_ratio_bound_or_point(&three).expect("point estimate");
378 assert!(r.lcb95.is_none());
379 let five = alternating(&[
380 (100.0, 90.0),
381 (101.0, 91.0),
382 (99.0, 89.0),
383 (100.5, 90.5),
384 (100.2, 90.1),
385 ]);
386 assert!(log_ratio_bound_or_point(&five)
387 .expect("bounded")
388 .lcb95
389 .is_some());
390 }
391}