edge_rs/ql/weights.rs
1//! Adjusted deviances, degrees of freedom and the average quasi-dispersion.
2//!
3//! Port of edgeR's `ql_glm.c`. [`crate::ql::chebyshev`] supplies the two
4//! moments of a unit deviance; this module spends them. Per gene it rescales
5//! each unit deviance by `alpha`, weights the complementary leverage by
6//! `kappa`, and sums over samples. Their ratio is `s2`, the quasi-likelihood
7//! dispersion that `glmQLFit` hands to `squeezeVar`.
8//!
9//! Each gene needs its own weighted design and QR, so the gene fan-out carries
10//! a per-thread scratch buffer.
11//!
12//! ### Deviations from edgePython
13//!
14//! [`compute_prior`] smooths with Cleveland's `lowess` as edgeR does, not
15//! limma's `weightedLowess`. See `UPSTREAM_DEVIATIONS.md` A11.
16//!
17//! ### Deviation from edgeR
18//!
19//! Leverages come from an unpivoted QR, not LINPACK's rank-revealing `dqrdc2`.
20//! The two differ only for rank-deficient weighted designs, which structural
21//! zeros would produce, but [`compute_weight`] returns `(0, 0)` for a zero mean
22//! so the difference cancels. See
23//! `test_structural_zeros_collapse_the_degrees_of_freedom`.
24//!
25//! ### References
26//!
27//! Lund, Nettleton, McCarthy and Smyth, SAGMB 11(5), 2012
28//! Chen, Lun and Smyth, F1000Research 5:1438, 2016
29
30use rayon::prelude::*;
31
32use crate::limma::lowess::lowess;
33use crate::numeric::stats::quantile_type7;
34use crate::prelude::*;
35use crate::ql::chebyshev::{compute_weight, unit_nb_deviance};
36use crate::utils::design::hat_diagonal;
37
38////////////
39// Consts //
40////////////
41
42/// Complementary leverage below which an observation is dropped outright.
43///
44/// edgeR's `thresholdzero`. A sample with leverage one is fitted exactly, so
45/// its residual carries no information.
46const THRESHOLD_ZERO: f64 = 1e-4;
47
48/// Degrees of freedom below which a gene is excluded from the prior trend.
49///
50/// edgeR's `t`. Genes at zero residual df have `s2 = 0` and would drag the
51/// lowess fit towards zero.
52const MIN_DF: f64 = 1e-8;
53
54/// Span of the lowess fit behind the prior.
55///
56/// edgeR's `f`: half the genes per window.
57const PRIOR_SPAN: f64 = 0.5;
58
59/// Robustness iterations behind the prior.
60///
61/// R's `lowess` counts robustness passes after the initial fit, as
62/// [`crate::limma::lowess::lowess`] does. limma's `weightedLowess` counts total
63/// passes, so its equivalent would be four.
64const PRIOR_ITERATIONS: usize = 3;
65
66/// Quantile of the fitted trend taken as the prior.
67///
68/// The 90th percentile, so the prior sits above the bulk of the genes.
69const PRIOR_QUANTILE: f64 = 0.9;
70
71/// Lower bound on the prior before it is raised to the fourth power.
72///
73/// A prior below one would sharpen the dispersions rather than shrink them.
74/// edgeR clamps on the fourth-root scale, so the floor on the returned value
75/// is also one.
76const PRIOR_FLOOR: f64 = 1.0;
77
78//////////////
79// Frontend //
80//////////////
81
82/// Per-gene adjusted deviance, degrees of freedom and their ratio.
83///
84/// All three vectors have one entry per gene, in input order.
85#[derive(Clone, Debug)]
86pub struct AdjustedDeviance {
87 /// Adjusted residual deviance, summed over samples.
88 pub deviance: Vec<f64>,
89 /// Adjusted residual degrees of freedom, summed over samples. Fractional
90 /// because each sample contributes its own `kappa`.
91 pub df: Vec<f64>,
92 /// Quasi-likelihood dispersion, `deviance / df`. Zero where `df` has
93 /// collapsed below `THRESHOLD_ZERO`.
94 pub s2: Vec<f64>,
95}
96
97/// Adjusted deviance and degrees of freedom for every gene.
98///
99/// Port of edgeR's `compute_adjust_vec`. Per gene and sample:
100///
101/// ```text
102/// deviance += alpha * d(y, mu, phi w / prior) * w
103/// df += kappa * (1 - h)
104/// ```
105///
106/// with `(alpha, kappa)` from [`compute_weight`], `d` the unit deviance and `h`
107/// the leverage under working weights `mu / (1 + mu phi w / prior)`.
108/// Observations with complementary leverage below `THRESHOLD_ZERO` drop out
109/// of both sums.
110///
111/// ### Params
112///
113/// * `counts` - Counts, row-major `n_genes * n_samples`
114/// * `n_genes` - Number of genes
115/// * `n_samples` - Number of samples
116/// * `fitted` - Fitted means, row-major `n_genes * n_samples`
117/// * `design` - Design matrix, row-major `n_samples * n_coef`
118/// * `n_coef` - Number of coefficients
119/// * `dispersion` - Negative binomial dispersions, recycled over genes and
120/// samples
121/// * `prior` - Average quasi-dispersion. Length one for edgeR's scalar prior,
122/// or `n_genes` for a per-gene one. Every entry must be finite and positive.
123/// * `weights` - Optional observation weights, recycled the same way
124///
125/// ### Returns
126///
127/// An [`AdjustedDeviance`], or [`EdgeErrors`] if any shape disagrees, the design
128/// has fewer rows than columns, or a dispersion or prior is out of domain.
129#[allow(clippy::too_many_arguments)]
130pub fn compute_adjust_vec<T: EdgeFloat>(
131 counts: &[T],
132 n_genes: usize,
133 n_samples: usize,
134 fitted: &[f64],
135 design: &[f64],
136 n_coef: usize,
137 dispersion: &Recycled<f64>,
138 prior: &[f64],
139 weights: Option<&Recycled<f64>>,
140) -> Result<AdjustedDeviance, EdgeErrors> {
141 validate(
142 counts.len(),
143 n_genes,
144 n_samples,
145 fitted.len(),
146 design.len(),
147 n_coef,
148 )?;
149 dispersion.validate(n_genes, n_samples)?;
150 if let Some(w) = weights {
151 w.validate(n_genes, n_samples)?;
152 }
153 validate_prior(prior, n_genes)?;
154
155 let per_gene: Vec<(f64, f64, f64)> = (0..n_genes)
156 .into_par_iter()
157 .map_init(
158 || Scratch::new(n_samples, n_coef),
159 |scratch, gene| {
160 let start = gene * n_samples;
161 for (slot, value) in scratch.y.iter_mut().zip(&counts[start..start + n_samples]) {
162 *slot = value.to_f64().unwrap_or(0.0);
163 }
164
165 adjust_one_gene(
166 &scratch.y,
167 &fitted[start..start + n_samples],
168 design,
169 n_coef,
170 dispersion.row(gene, n_samples),
171 weights.map(|w| w.row(gene, n_samples)),
172 if prior.len() == 1 {
173 prior[0]
174 } else {
175 prior[gene]
176 },
177 &mut scratch.weighted_design,
178 )
179 },
180 )
181 .collect();
182
183 let mut out = AdjustedDeviance {
184 deviance: Vec::with_capacity(n_genes),
185 df: Vec::with_capacity(n_genes),
186 s2: Vec::with_capacity(n_genes),
187 };
188 for (deviance, df, s2) in per_gene {
189 out.deviance.push(deviance);
190 out.df.push(df);
191 out.s2.push(s2);
192 }
193
194 Ok(out)
195}
196
197/// Average quasi-dispersion from one round of quasi-likelihood dispersions.
198///
199/// Port of edgeR's `compute_prior`. Genes with usable residual df are smoothed
200/// on the fourth-root scale against average log-CPM. The 90th percentile of
201/// the trend is floored at one and raised back to the fourth power.
202///
203/// ### Params
204///
205/// * `ave_log_cpm` - Average log-CPM, one per gene
206/// * `s2` - Quasi-likelihood dispersions, one per gene
207/// * `df` - Adjusted residual degrees of freedom, one per gene
208///
209/// ### Returns
210///
211/// A single-element vector with the average quasi-dispersion, as
212/// [`compute_adjust_vec`] takes for `prior`. Errors with
213/// [`EdgeErrors::LengthMismatch`] if the inputs disagree.
214pub fn compute_prior(ave_log_cpm: &[f64], s2: &[f64], df: &[f64]) -> Result<Vec<f64>, EdgeErrors> {
215 let n = ave_log_cpm.len();
216 if s2.len() != n {
217 return Err(EdgeErrors::LengthMismatch {
218 name: "s2",
219 expected: n,
220 got: s2.len(),
221 });
222 }
223 if df.len() != n {
224 return Err(EdgeErrors::LengthMismatch {
225 name: "df",
226 expected: n,
227 got: df.len(),
228 });
229 }
230
231 let mut abundance = Vec::with_capacity(n);
232 let mut root_s2 = Vec::with_capacity(n);
233 for gene in 0..n {
234 if df[gene] > MIN_DF {
235 abundance.push(ave_log_cpm[gene]);
236 root_s2.push(s2[gene].max(0.0).sqrt().sqrt());
237 }
238 }
239
240 let trend = match root_s2.len() {
241 0 => return Ok(vec![PRIOR_FLOOR]),
242 1 => root_s2,
243 _ => {
244 // edgeR's `compute_ave_qd` calls `lowess(f = 0.5, iter = 3)`.
245 lowess(&abundance, &root_s2, PRIOR_SPAN, PRIOR_ITERATIONS)?
246 }
247 };
248
249 let p = quantile_type7(&trend, PRIOR_QUANTILE)?.max(PRIOR_FLOOR);
250
251 Ok(vec![p * p * p * p])
252}
253
254/// Average quasi-dispersion, found by two rounds of adjustment.
255///
256/// Port of edgeR's `update_prior`. Starts at one, then twice adjusts and
257/// re-estimates. Two rounds, not iterated to convergence, as in edgeR.
258///
259/// ### Params
260///
261/// * `counts` - Counts, row-major `n_genes * n_samples`
262/// * `n_genes` - Number of genes
263/// * `n_samples` - Number of samples
264/// * `fitted` - Fitted means, row-major `n_genes * n_samples`
265/// * `design` - Design matrix, row-major `n_samples * n_coef`
266/// * `n_coef` - Number of coefficients
267/// * `dispersion` - Negative binomial dispersions, recycled over genes and
268/// samples
269/// * `weights` - Optional observation weights, recycled the same way
270/// * `ave_log_cpm` - Average log-CPM, one per gene
271///
272/// ### Returns
273///
274/// A single-element vector holding the average quasi-dispersion, or
275/// [`EdgeErrors`] if any shape disagrees.
276#[allow(clippy::too_many_arguments)]
277pub fn update_prior<T: EdgeFloat>(
278 counts: &[T],
279 n_genes: usize,
280 n_samples: usize,
281 fitted: &[f64],
282 design: &[f64],
283 n_coef: usize,
284 dispersion: &Recycled<f64>,
285 weights: Option<&Recycled<f64>>,
286 ave_log_cpm: &[f64],
287) -> Result<Vec<f64>, EdgeErrors> {
288 if ave_log_cpm.len() != n_genes {
289 return Err(EdgeErrors::LengthMismatch {
290 name: "ave_log_cpm",
291 expected: n_genes,
292 got: ave_log_cpm.len(),
293 });
294 }
295
296 let mut prior = vec![PRIOR_FLOOR];
297 for _ in 0..2 {
298 let adjusted = compute_adjust_vec(
299 counts, n_genes, n_samples, fitted, design, n_coef, dispersion, &prior, weights,
300 )?;
301 prior = compute_prior(ave_log_cpm, &adjusted.s2, &adjusted.df)?;
302 }
303
304 Ok(prior)
305}
306
307/////////////
308// Scratch //
309/////////////
310
311/// Per-thread buffers for the gene fan-out.
312struct Scratch {
313 /// Counts for the current gene, widened to `f64`.
314 y: Vec<f64>,
315 /// `sqrt(W) X` for the current gene, row-major `n_samples * n_coef`.
316 weighted_design: Vec<f64>,
317}
318
319impl Scratch {
320 /// Allocates scratch for one worker.
321 ///
322 /// ### Params
323 ///
324 /// * `n_samples` - Number of samples
325 /// * `n_coef` - Number of coefficients
326 ///
327 /// ### Returns
328 ///
329 /// Buffers sized for a single gene.
330 fn new(n_samples: usize, n_coef: usize) -> Self {
331 Self {
332 y: vec![0.0; n_samples],
333 weighted_design: vec![0.0; n_samples * n_coef],
334 }
335 }
336}
337
338/// Checks the shapes [`compute_adjust_vec`] cannot recover from.
339///
340/// Rejects everything [`hat_diagonal`] would, so the per-gene kernel can treat
341/// it as infallible.
342///
343/// ### Params
344///
345/// * `n_counts` - Length of the counts slice
346/// * `n_genes` - Number of genes
347/// * `n_samples` - Number of samples
348/// * `n_fitted` - Length of the fitted slice
349/// * `n_design` - Length of the design slice
350/// * `n_coef` - Number of coefficients
351///
352/// ### Returns
353///
354/// `Ok(())` when every shape is consistent, otherwise [`EdgeErrors`].
355fn validate(
356 n_counts: usize,
357 n_genes: usize,
358 n_samples: usize,
359 n_fitted: usize,
360 n_design: usize,
361 n_coef: usize,
362) -> Result<(), EdgeErrors> {
363 if n_genes == 0 || n_samples == 0 {
364 return Err(EdgeErrors::EmptyCounts { n_genes, n_samples });
365 }
366 if n_coef == 0 {
367 return Err(EdgeErrors::MustBePositive("n_coef".to_string()));
368 }
369 if n_counts != n_genes * n_samples {
370 return Err(EdgeErrors::LengthMismatch {
371 name: "counts",
372 expected: n_genes * n_samples,
373 got: n_counts,
374 });
375 }
376 if n_fitted != n_genes * n_samples {
377 return Err(EdgeErrors::LengthMismatch {
378 name: "fitted",
379 expected: n_genes * n_samples,
380 got: n_fitted,
381 });
382 }
383 if n_design != n_samples * n_coef {
384 return Err(EdgeErrors::LengthMismatch {
385 name: "design",
386 expected: n_samples * n_coef,
387 got: n_design,
388 });
389 }
390 if n_samples < n_coef {
391 return Err(EdgeErrors::DesignNotFullRank {
392 n_cols: n_coef,
393 rank: n_samples,
394 });
395 }
396 Ok(())
397}
398
399/// Checks the prior is either scalar or one value per gene, and in domain.
400///
401/// ### Params
402///
403/// * `prior` - Average quasi-dispersion, length one or `n_genes`
404/// * `n_genes` - Number of genes
405///
406/// ### Returns
407///
408/// `Ok(())` when usable, otherwise [`EdgeErrors`]. The prior divides a
409/// dispersion, so zero and negative values are rejected up front.
410fn validate_prior(prior: &[f64], n_genes: usize) -> Result<(), EdgeErrors> {
411 if prior.len() != 1 && prior.len() != n_genes {
412 return Err(EdgeErrors::LengthMismatch {
413 name: "prior",
414 expected: n_genes,
415 got: prior.len(),
416 });
417 }
418 if let Some(&bad) = prior.iter().find(|p| !p.is_finite() || **p <= 0.0) {
419 return Err(EdgeErrors::InvalidArgument(format!(
420 "prior must be finite and positive; got {bad}."
421 )));
422 }
423 Ok(())
424}
425
426/// Adjusted deviance, degrees of freedom and their ratio for one gene.
427///
428/// Leverages come from the QR of `sqrt(W) X` with working weight
429/// `mu / (1 + mu phi w / prior)`. edgeR's `qr_hat` uses LINPACK's `dqrdc2`,
430/// which truncates at the numerical rank; [`hat_diagonal`] does not. The two
431/// differ only for a rank-deficient weighted design with non-negligible column
432/// norms, which needs a tiny non-zero fitted mean. edgeR's GLM returns exact
433/// zeros for structural zeros and [`compute_weight`] gives `(0, 0)` there, so
434/// the difference cancels. See
435/// `test_structural_zeros_collapse_the_degrees_of_freedom`.
436///
437/// ### Params
438///
439/// * `y` - Counts for this gene, one per sample
440/// * `mu` - Fitted means for this gene, same length
441/// * `design` - Design matrix, row-major `n_samples * n_coef`
442/// * `n_coef` - Number of coefficients
443/// * `dispersion` - Dispersion row for this gene
444/// * `weights` - Optional weight row for this gene
445/// * `prior` - Average quasi-dispersion for this gene
446/// * `weighted_design` - Scratch of `n_samples * n_coef`, overwritten
447///
448/// ### Returns
449///
450/// `(deviance, df, s2)`.
451#[allow(clippy::too_many_arguments)]
452fn adjust_one_gene(
453 y: &[f64],
454 mu: &[f64],
455 design: &[f64],
456 n_coef: usize,
457 dispersion: RecycledRow<'_, f64>,
458 weights: Option<RecycledRow<'_, f64>>,
459 prior: f64,
460 weighted_design: &mut [f64],
461) -> (f64, f64, f64) {
462 let n_samples = y.len();
463
464 for (sample, &mu_j) in mu.iter().enumerate() {
465 let phi = dispersion.get(sample);
466 let w = weights.map_or(1.0, |w| w.get(sample));
467 let working = (mu_j / (1.0 + (mu_j * phi * w / prior))).sqrt();
468 let row = &design[sample * n_coef..(sample + 1) * n_coef];
469 let out = &mut weighted_design[sample * n_coef..(sample + 1) * n_coef];
470 for (slot, &x) in out.iter_mut().zip(row) {
471 *slot = x * working;
472 }
473 }
474
475 // Infallible: `validate` has already rejected every shape `hat_diagonal`
476 // checks, and the scratch is sized from those same dimensions.
477 let hat = hat_diagonal(weighted_design, n_samples, n_coef)
478 .expect("weighted design shape was validated on entry");
479
480 let mut deviance = 0.0;
481 let mut df = 0.0;
482 for sample in 0..n_samples {
483 let phi = dispersion.get(sample);
484 let w = weights.map_or(1.0, |w| w.get(sample));
485 let (alpha, kappa) = compute_weight(mu[sample], phi, prior / w);
486 let mut unit = unit_nb_deviance(y[sample], mu[sample], phi * w / prior);
487 let mut complement = 1.0 - hat[sample];
488 if complement < THRESHOLD_ZERO {
489 unit = 0.0;
490 complement = 0.0;
491 }
492 deviance += (unit * alpha) * w;
493 df += complement * kappa;
494 }
495
496 let s2 = if df < THRESHOLD_ZERO {
497 0.0
498 } else {
499 deviance / df
500 };
501
502 (deviance, df, s2)
503}
504
505///////////
506// Tests //
507///////////
508
509#[cfg(test)]
510mod tests {
511 use super::*;
512 use approx::assert_relative_eq;
513
514 /// Relative tolerance against edgeR's C. The port reproduces the arithmetic
515 /// operation for operation, so the only slack is the QR.
516 const TOL: f64 = 1e-12;
517
518 /// Six genes by six samples, all values exactly representable in binary.
519 ///
520 /// Gene 3 is zero in the first group and gene 4 is zero everywhere (exact
521 /// zero fitted values from `glmFit`). Gene 5 sits in the thousands, past
522 /// the Chebyshev panel switch.
523 fn counts() -> Vec<f64> {
524 vec![
525 10.0, 12.0, 11.0, 40.0, 44.0, 38.0, //
526 50.0, 48.0, 52.0, 49.0, 51.0, 50.0, //
527 2.0, 0.0, 5.0, 1.0, 3.0, 0.0, //
528 0.0, 0.0, 0.0, 20.0, 22.0, 18.0, //
529 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //
530 1000.0, 1100.0, 900.0, 3000.0, 3200.0, 2800.0,
531 ]
532 }
533
534 /// Fitted means for [`counts`], all dyadic rationals.
535 fn fitted() -> Vec<f64> {
536 vec![
537 11.0, 11.5, 10.75, 40.5, 41.25, 39.5, //
538 50.0, 49.5, 50.25, 49.75, 50.5, 50.0, //
539 2.5, 2.25, 2.75, 1.5, 1.25, 1.75, //
540 0.0, 0.0, 0.0, 20.0, 20.0, 20.0, //
541 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //
542 1000.0, 1050.0, 950.0, 3000.0, 3100.0, 2900.0,
543 ]
544 }
545
546 /// Two-group design, row-major.
547 fn two_group() -> Vec<f64> {
548 vec![
549 1.0, 0.0, //
550 1.0, 0.0, //
551 1.0, 0.0, //
552 1.0, 1.0, //
553 1.0, 1.0, //
554 1.0, 1.0,
555 ]
556 }
557
558 /// Continuous covariate design, row-major.
559 fn continuous() -> Vec<f64> {
560 vec![
561 1.0, -1.5, //
562 1.0, -0.5, //
563 1.0, 0.5, //
564 1.0, 1.5, //
565 1.0, 2.5, //
566 1.0, 3.5,
567 ]
568 }
569
570 /// Runs the six-gene fixture and compares against a reference triple.
571 fn check(
572 design: &[f64],
573 dispersion: Recycled<f64>,
574 prior: &[f64],
575 weights: Option<&Recycled<f64>>,
576 deviance: &[f64],
577 df: &[f64],
578 s2: &[f64],
579 ) {
580 let got = compute_adjust_vec(
581 &counts(),
582 6,
583 6,
584 &fitted(),
585 design,
586 2,
587 &dispersion,
588 prior,
589 weights,
590 )
591 .unwrap();
592
593 for gene in 0..6 {
594 assert_relative_eq!(got.deviance[gene], deviance[gene], max_relative = TOL);
595 assert_relative_eq!(got.df[gene], df[gene], max_relative = TOL);
596 assert_relative_eq!(got.s2[gene], s2[gene], max_relative = TOL);
597 }
598 }
599
600 /// Shared preamble of the `Rscript` blocks below.
601 ///
602 /// ```r
603 /// # Rscript, edgeR 4.8.2
604 /// library(edgeR)
605 /// y <- matrix(c(10,12,11,40,44,38, 50,48,52,49,51,50, 2,0,5,1,3,0,
606 /// 0,0,0,20,22,18, 0,0,0,0,0,0,
607 /// 1000,1100,900,3000,3200,2800), nrow = 6, byrow = TRUE)
608 /// mu <- matrix(c(11,11.5,10.75,40.5,41.25,39.5, 50,49.5,50.25,49.75,50.5,50,
609 /// 2.5,2.25,2.75,1.5,1.25,1.75, 0,0,0,20,20,20, 0,0,0,0,0,0,
610 /// 1000,1050,950,3000,3100,2900), nrow = 6, byrow = TRUE)
611 /// storage.mode(y) <- "double"
612 /// X <- cbind(1, c(0,0,0,1,1,1)); storage.mode(X) <- "double"
613 /// Xc <- cbind(1, c(-1.5,-0.5,0.5,1.5,2.5,3.5)); storage.mode(Xc) <- "double"
614 /// dm <- edgeR:::.compressDispersions(y, 0.125)
615 /// wm <- edgeR:::.compressWeights(y, NULL)
616 /// ```
617 ///
618 /// `.Call(edgeR:::.cxx_compute_adj_vec, y, mu, X, dm, 1.0, wm)`
619 #[test]
620 fn test_two_group_matches_edger() {
621 check(
622 &two_group(),
623 Recycled::scalar(0.125),
624 &[1.0],
625 None,
626 &[
627 0.087_678_574_155_993_14,
628 0.016_535_319_784_658_36,
629 11.401_625_661_199_109,
630 0.111_548_751_889_321_45,
631 0.0,
632 0.056_972_440_150_977_24,
633 ],
634 &[
635 3.969_634_764_628_127,
636 3.999_429_434_054_224,
637 5.340_746_228_772_16,
638 1.992_229_288_606_397_7,
639 0.0,
640 4.001_942_003_921_18,
641 ],
642 &[
643 0.022_087_315_169_965_472,
644 0.004_134_419_685_934_175,
645 2.134_837_562_544_204,
646 0.055_991_924_487_442_87,
647 0.0,
648 0.014_236_198_349_489_956,
649 ],
650 );
651 }
652
653 /// Prior above one, the branch `update_prior` lands on. The prior enters the
654 /// working weights, moments and unit deviance separately.
655 ///
656 /// `.Call(edgeR:::.cxx_compute_adj_vec, y, mu, X, dm, 2.5, wm)`
657 #[test]
658 fn test_prior_above_one_matches_edger() {
659 check(
660 &two_group(),
661 Recycled::scalar(0.125),
662 &[2.5],
663 None,
664 &[
665 0.148_404_728_650_759_82,
666 0.034_023_627_382_698_46,
667 23.944_790_215_814_272,
668 0.187_876_878_253_122_5,
669 0.0,
670 0.141_051_250_719_609_08,
671 ],
672 &[
673 3.944_608_502_302_763_5,
674 3.984_625_332_836_51,
675 9.508_568_001_899_926,
676 1.949_129_177_816_912_6,
677 0.0,
678 4.001_937_312_693_502,
679 ],
680 &[
681 0.037_622_169_237_866_03,
682 0.008_538_726_866_567_972,
683 2.518_233_051_604_596_5,
684 0.096_390_162_535_841_08,
685 0.0,
686 0.035_245_742_173_975_886,
687 ],
688 );
689 }
690
691 /// Per-gene dispersions spanning the Chebyshev cases: gene 5 is at two, past
692 /// case 1, and gene 4 is zero (inert, no counts).
693 ///
694 /// ```r
695 /// dg <- edgeR:::.compressDispersions(y, c(0.05,0.125,0.5,0.25,0,2))
696 /// .Call(edgeR:::.cxx_compute_adj_vec, y, mu, X, dg, 1.5, wm)
697 /// ```
698 #[test]
699 fn test_per_gene_dispersion_matches_edger() {
700 check(
701 &two_group(),
702 Recycled::by_gene(vec![0.05, 0.125, 0.5, 0.25, 0.0, 2.0]),
703 &[1.5],
704 None,
705 &[
706 0.185_285_564_462_572_72,
707 0.023_166_881_640_206_673,
708 14.989_239_815_832_821,
709 0.086_845_649_266_648_48,
710 0.0,
711 0.005_175_385_846_847_329_5,
712 ],
713 &[
714 3.954_646_135_693_23,
715 3.996_255_259_784_583_3,
716 8.051_884_164_459_716,
717 1.968_403_109_458_039_6,
718 0.0,
719 4.722_405_574_548_103_5,
720 ],
721 &[
722 0.046_852_628_049_385_024,
723 0.005_797_147_613_001_947,
724 1.861_581_650_912_709_4,
725 0.044_119_849_663_598_47,
726 0.0,
727 0.001_095_921_509_736_607_7,
728 ],
729 );
730 }
731
732 /// Observation weights, which enter the working weight, the dispersion and
733 /// the multiplier on the adjusted unit deviance.
734 ///
735 /// ```r
736 /// w <- matrix(c(1,1,1,1,1,1, 0.5,0.5,2,2,1,1, 1,0.25,1,1,4,1,
737 /// 2,2,2,0.5,0.5,0.5, 1,1,1,1,1,1, 1,2,1,0.5,1,1),
738 /// nrow = 6, byrow = TRUE)
739 /// wm2 <- edgeR:::.compressWeights(y, w)
740 /// .Call(edgeR:::.cxx_compute_adj_vec, y, mu, X, dm, 1.25, wm2)
741 /// ```
742 #[test]
743 fn test_weights_match_edger() {
744 let weights = Recycled::full(
745 vec![
746 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, //
747 0.5, 0.5, 2.0, 2.0, 1.0, 1.0, //
748 1.0, 0.25, 1.0, 1.0, 4.0, 1.0, //
749 2.0, 2.0, 2.0, 0.5, 0.5, 0.5, //
750 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, //
751 1.0, 2.0, 1.0, 0.5, 1.0, 1.0,
752 ],
753 6,
754 6,
755 )
756 .unwrap();
757
758 check(
759 &two_group(),
760 Recycled::scalar(0.125),
761 &[1.25],
762 Some(&weights),
763 &[
764 0.101_079_410_203_347_93,
765 0.019_878_961_611_844_97,
766 12.984_719_805_109_156,
767 0.093_938_439_126_561_25,
768 0.0,
769 0.071_201_367_293_380_21,
770 ],
771 &[
772 3.954_058_331_932_422_7,
773 3.995_405_261_291_514,
774 6.298_973_859_851_691,
775 1.949_129_177_816_912_6,
776 0.0,
777 4.001_941_748_985_114,
778 ],
779 &[
780 0.025_563_459_544_095_425,
781 0.004_975_455_632_608_117,
782 2.061_402_395_693_523_5,
783 0.048_195_081_267_920_54,
784 0.0,
785 0.017_791_705_067_030_712,
786 ],
787 );
788 }
789
790 /// Continuous covariate. Only the leverages change, so deviances match the
791 /// two-group case bit for bit while the df do not.
792 ///
793 /// `.Call(edgeR:::.cxx_compute_adj_vec, y, mu, Xc, dm, 1.0, wm)`
794 #[test]
795 fn test_continuous_design_matches_edger() {
796 check(
797 &continuous(),
798 Recycled::scalar(0.125),
799 &[1.0],
800 None,
801 &[
802 0.087_678_574_155_993_14,
803 0.016_535_319_784_658_36,
804 11.401_625_661_199_109,
805 0.111_548_751_889_321_45,
806 0.0,
807 0.056_972_440_150_977_24,
808 ],
809 &[
810 3.968_981_768_078_912_3,
811 3.999_429_410_347_023_6,
812 5.390_514_507_550_15,
813 0.996_114_644_303_199,
814 0.0,
815 4.001_941_999_230_005,
816 ],
817 &[
818 0.022_090_949_084_513_38,
819 0.004_134_419_710_441_55,
820 2.115_127_534_714_836,
821 0.111_983_848_974_885_72,
822 0.0,
823 0.014_236_198_366_177_98,
824 ],
825 );
826 }
827
828 /// One gene alone must match the same gene inside a batch.
829 ///
830 /// ```r
831 /// y1 <- y[2,,drop=FALSE]; mu1 <- mu[2,,drop=FALSE]
832 /// dm1 <- edgeR:::.compressDispersions(y1, 0.125)
833 /// wm1 <- edgeR:::.compressWeights(y1, NULL)
834 /// .Call(edgeR:::.cxx_compute_adj_vec, y1, mu1, X, dm1, 1.0, wm1)
835 /// ```
836 #[test]
837 fn test_single_gene_matches_edger() {
838 let out = compute_adjust_vec(
839 &counts()[6..12],
840 1,
841 6,
842 &fitted()[6..12],
843 &two_group(),
844 2,
845 &Recycled::scalar(0.125),
846 &[1.0],
847 None,
848 )
849 .unwrap();
850
851 assert_relative_eq!(
852 out.deviance[0],
853 0.016_535_319_784_658_36,
854 max_relative = TOL
855 );
856 assert_relative_eq!(out.df[0], 3.999_429_434_054_224, max_relative = TOL);
857 assert_relative_eq!(out.s2[0], 0.004_134_419_685_934_175, max_relative = TOL);
858 }
859
860 /// Genes 3 and 4 are the structural-zero cases that make the unpivoted QR
861 /// safe.
862 ///
863 /// edgeR's `dqrdc2` finds rank 1 for gene 3 and rank 0 for gene 4;
864 /// [`hat_diagonal`] reports different leverages. The result is unchanged
865 /// because [`compute_weight`] returns `(0, 0)` at a zero mean. Gene 4 is
866 /// inert and gene 3 keeps exactly the three samples that carry counts.
867 #[test]
868 fn test_structural_zeros_collapse_the_degrees_of_freedom() {
869 let out = compute_adjust_vec(
870 &counts(),
871 6,
872 6,
873 &fitted(),
874 &two_group(),
875 2,
876 &Recycled::scalar(0.125),
877 &[1.0],
878 None,
879 )
880 .unwrap();
881
882 assert_eq!(out.deviance[4], 0.0);
883 assert_eq!(out.df[4], 0.0);
884 assert_eq!(out.s2[4], 0.0);
885 // Three live samples against one fitted coefficient.
886 assert!(out.df[3] > 1.9 && out.df[3] < 2.1);
887 // Genes with no structural zeros keep close to four.
888 assert!(out.df[0] > 3.9);
889 }
890
891 /// The three `Recycled` forms of one dispersion must agree.
892 #[test]
893 fn test_recycled_dispersion_forms_agree() {
894 let run = |dispersion: Recycled<f64>| {
895 compute_adjust_vec(
896 &counts(),
897 6,
898 6,
899 &fitted(),
900 &two_group(),
901 2,
902 &dispersion,
903 &[1.0],
904 None,
905 )
906 .unwrap()
907 };
908
909 let scalar = run(Recycled::scalar(0.125));
910 let per_gene = run(Recycled::by_gene(vec![0.125; 6]));
911 let by_sample = run(Recycled::by_sample(vec![0.125; 6]));
912 let full = run(Recycled::full(vec![0.125; 36], 6, 6).unwrap());
913
914 assert_eq!(scalar.deviance, per_gene.deviance);
915 assert_eq!(scalar.deviance, by_sample.deviance);
916 assert_eq!(scalar.deviance, full.deviance);
917 assert_eq!(scalar.df, per_gene.df);
918 assert_eq!(scalar.df, full.df);
919 }
920
921 /// A per-gene prior of a repeated value must equal the scalar prior.
922 #[test]
923 fn test_per_gene_prior_matches_scalar_prior() {
924 let run = |prior: &[f64]| {
925 compute_adjust_vec(
926 &counts(),
927 6,
928 6,
929 &fitted(),
930 &two_group(),
931 2,
932 &Recycled::scalar(0.125),
933 prior,
934 None,
935 )
936 .unwrap()
937 };
938
939 let scalar = run(&[2.5]);
940 let vector = run(&[2.5; 6]);
941 assert_eq!(scalar.deviance, vector.deviance);
942 assert_eq!(scalar.df, vector.df);
943 }
944
945 /// `f32` counts give the same answer (all fixture counts are exact in `f32`).
946 #[test]
947 fn test_f32_counts_agree_with_f64() {
948 let wide = compute_adjust_vec(
949 &counts(),
950 6,
951 6,
952 &fitted(),
953 &two_group(),
954 2,
955 &Recycled::scalar(0.125),
956 &[1.0],
957 None,
958 )
959 .unwrap();
960
961 let narrow: Vec<f32> = counts().iter().map(|&v| v as f32).collect();
962 let got = compute_adjust_vec(
963 &narrow,
964 6,
965 6,
966 &fitted(),
967 &two_group(),
968 2,
969 &Recycled::scalar(0.125),
970 &[1.0],
971 None,
972 )
973 .unwrap();
974
975 assert_eq!(wide.deviance, got.deviance);
976 assert_eq!(wide.df, got.df);
977 }
978
979 ////////////////////
980 // Prior fixtures //
981 ////////////////////
982
983 /// Relative tolerance for anything that passes through the smoother.
984 /// Rounding only, since `compute_prior` uses the same `lowess` as edgeR.
985 const PRIOR_TOL: f64 = 1e-12;
986
987 /// Twenty-five genes on a dyadic grid of average log-CPM, with dyadic `s2`
988 /// and three df each.
989 fn prior_fixture() -> (Vec<f64>, Vec<f64>, Vec<f64>) {
990 let ave: Vec<f64> = (0..25).map(|i| i as f64 / 4.0 - 2.0).collect();
991 let s2: Vec<f64> = (0..25)
992 .map(|i| ((i * 7 % 13) + 1) as f64 / 8.0 + i as f64 / 16.0)
993 .collect();
994 (ave, s2, vec![3.0; 25])
995 }
996
997 /// The prior against edgeR's `compute_prior`, replicated in R (the C
998 /// function is not exported).
999 ///
1000 /// ```r
1001 /// # Rscript, R 4.5.1 / edgeR 4.8.2
1002 /// ave <- (0:24)/4 - 2
1003 /// s2 <- ((0:24 * 7) %% 13 + 1)/8 + (0:24)/16
1004 /// ans <- lowess(ave, sqrt(sqrt(s2)), f = 0.5, iter = 3,
1005 /// delta = 0.01 * diff(range(ave)))$y
1006 /// max(quantile(ans, 0.9, type = 7, names = FALSE), 1)^4 # 2.423348828030909
1007 /// ```
1008 #[test]
1009 fn test_compute_prior_matches_edger() {
1010 let (ave, s2, df) = prior_fixture();
1011 let got = compute_prior(&ave, &s2, &df).unwrap();
1012
1013 assert_eq!(got.len(), 1);
1014 assert_relative_eq!(got[0], 2.423_348_828_030_909, max_relative = PRIOR_TOL);
1015 }
1016
1017 /// Genes below the df threshold never reach the smoother: dropping the two
1018 /// end genes must equal fitting the survivors alone.
1019 ///
1020 /// ```r
1021 /// a2 <- lowess(ave[2:24], sqrt(sqrt(s2[2:24])), f = 0.5, iter = 3,
1022 /// delta = 0.01 * diff(range(ave[2:24])))$y
1023 /// max(quantile(a2, 0.9, type = 7, names = FALSE), 1)^4 # 2.2385492590160405
1024 /// ```
1025 #[test]
1026 fn test_compute_prior_drops_low_df_genes() {
1027 let (ave, s2, mut df) = prior_fixture();
1028 df[0] = 0.0;
1029 df[24] = MIN_DF;
1030 let filtered = compute_prior(&ave, &s2, &df).unwrap();
1031
1032 let direct = compute_prior(&ave[1..24], &s2[1..24], &[3.0; 23]).unwrap();
1033 assert_eq!(filtered, direct);
1034 assert_relative_eq!(
1035 filtered[0],
1036 2.238_549_259_016_040_5,
1037 max_relative = PRIOR_TOL
1038 );
1039 }
1040
1041 /// No gene left gives the floor; one gene gives its own value. Both are
1042 /// below the smoother's two-point minimum.
1043 #[test]
1044 fn test_compute_prior_degenerate_gene_counts() {
1045 let none = compute_prior(&[1.0, 2.0], &[4.0, 9.0], &[0.0, 0.0]).unwrap();
1046 assert_eq!(none, vec![1.0]);
1047
1048 let one = compute_prior(&[1.0, 2.0], &[16.0, 9.0], &[3.0, 0.0]).unwrap();
1049 assert_relative_eq!(one[0], 16.0, max_relative = 1e-12);
1050 }
1051
1052 /// A trend that sits entirely below one is floored, not inverted.
1053 #[test]
1054 fn test_compute_prior_floors_at_one() {
1055 let ave: Vec<f64> = (0..10).map(|i| i as f64).collect();
1056 let got = compute_prior(&ave, &[0.0625; 10], &[2.0; 10]).unwrap();
1057 assert_eq!(got, vec![1.0]);
1058 }
1059
1060 /// Twenty-four overdispersed genes by six samples, so the prior lands well
1061 /// above its floor. Fitted means are group means rounded to a quarter.
1062 ///
1063 /// Generated by the script quoted on [`test_update_prior_matches_edger`].
1064 #[allow(clippy::type_complexity)]
1065 fn prior_counts() -> (Vec<f64>, Vec<f64>, Vec<f64>) {
1066 let counts = vec![
1067 0.0, 11.0, 1.0, 2.0, 2.0, 2.0, //
1068 1.0, 5.0, 5.0, 8.0, 3.0, 3.0, //
1069 1.0, 3.0, 5.0, 7.0, 3.0, 3.0, //
1070 7.0, 6.0, 0.0, 6.0, 3.0, 6.0, //
1071 3.0, 9.0, 1.0, 4.0, 8.0, 0.0, //
1072 18.0, 11.0, 9.0, 36.0, 11.0, 5.0, //
1073 19.0, 9.0, 7.0, 1.0, 34.0, 17.0, //
1074 20.0, 3.0, 7.0, 22.0, 3.0, 4.0, //
1075 19.0, 14.0, 40.0, 64.0, 37.0, 26.0, //
1076 12.0, 20.0, 27.0, 9.0, 28.0, 54.0, //
1077 41.0, 8.0, 2.0, 4.0, 22.0, 68.0, //
1078 32.0, 38.0, 43.0, 22.0, 13.0, 46.0, //
1079 46.0, 92.0, 8.0, 53.0, 49.0, 66.0, //
1080 13.0, 7.0, 88.0, 28.0, 156.0, 51.0, //
1081 294.0, 31.0, 19.0, 17.0, 158.0, 20.0, //
1082 22.0, 73.0, 14.0, 74.0, 160.0, 42.0, //
1083 76.0, 154.0, 67.0, 123.0, 28.0, 265.0, //
1084 196.0, 111.0, 29.0, 150.0, 18.0, 96.0, //
1085 7.0, 133.0, 102.0, 174.0, 310.0, 73.0, //
1086 133.0, 492.0, 362.0, 14.0, 71.0, 58.0, //
1087 250.0, 551.0, 40.0, 46.0, 54.0, 56.0, //
1088 45.0, 288.0, 209.0, 376.0, 182.0, 152.0, //
1089 458.0, 157.0, 192.0, 524.0, 480.0, 574.0, //
1090 172.0, 271.0, 340.0, 874.0, 456.0, 210.0,
1091 ];
1092
1093 // Group means, three samples each.
1094 let group_means: [(f64, f64); 24] = [
1095 (4.0, 2.0),
1096 (3.75, 4.75),
1097 (3.0, 4.25),
1098 (4.25, 5.0),
1099 (4.25, 4.0),
1100 (12.75, 17.25),
1101 (11.75, 17.25),
1102 (10.0, 9.75),
1103 (24.25, 42.25),
1104 (19.75, 30.25),
1105 (17.0, 31.25),
1106 (37.75, 27.0),
1107 (48.75, 56.0),
1108 (36.0, 78.25),
1109 (114.75, 65.0),
1110 (36.25, 92.0),
1111 (99.0, 138.75),
1112 (112.0, 88.0),
1113 (80.75, 185.75),
1114 (329.0, 47.75),
1115 (280.25, 52.0),
1116 (180.75, 236.75),
1117 (269.0, 526.0),
1118 (261.0, 513.25),
1119 ];
1120 let mut fitted = Vec::with_capacity(144);
1121 for (first, second) in group_means {
1122 fitted.extend_from_slice(&[first, first, first, second, second, second]);
1123 }
1124
1125 let ave_log_cpm = vec![
1126 11.171875, 11.484375, 11.375, 11.59375, 11.46875, 12.921875, 12.9375, 12.4375,
1127 14.015625, 13.6875, 13.65625, 14.046875, 14.59375, 14.828125, 15.453125, 14.84375,
1128 15.8125, 15.546875, 15.890625, 16.515625, 16.203125, 16.546875, 17.53125, 17.421875,
1129 ];
1130
1131 (counts, fitted, ave_log_cpm)
1132 }
1133
1134 /// One adjustment round plus one `compute_prior`: the first half of
1135 /// [`update_prior`].
1136 ///
1137 /// ```r
1138 /// o1 <- .Call(edgeR:::.cxx_compute_adj_vec, y, mu, X, dm, 1.0, wm)
1139 /// cp(ave, o1$s2, o1$df) # 5.2505569919690807
1140 /// ```
1141 #[test]
1142 fn test_prior_after_one_round_matches_edger() {
1143 let (counts, fitted, ave) = prior_counts();
1144 let adjusted = compute_adjust_vec(
1145 &counts,
1146 24,
1147 6,
1148 &fitted,
1149 &two_group(),
1150 2,
1151 &Recycled::scalar(0.125),
1152 &[1.0],
1153 None,
1154 )
1155 .unwrap();
1156
1157 let prior = compute_prior(&ave, &adjusted.s2, &adjusted.df).unwrap();
1158 assert_relative_eq!(prior[0], 5.250_556_991_969_081, max_relative = PRIOR_TOL);
1159 }
1160
1161 /// The two-round loop against edgeR's `compute_ave_qd`, the entry point of
1162 /// `glmQLFit(legacy = FALSE)`.
1163 ///
1164 /// ```r
1165 /// # Rscript, R 4.5.1 / edgeR 4.8.2
1166 /// set.seed(11)
1167 /// ng <- 24; ns <- 6
1168 /// X <- cbind(1, rep(c(0,1), each = 3)); storage.mode(X) <- "double"
1169 /// base <- round(exp(seq(log(4), log(400), length.out = ng)))
1170 /// y <- matrix(0, ng, ns)
1171 /// for (g in 1:ng) y[g,] <- rnbinom(ns, size = 1/0.6, mu = base[g])
1172 /// storage.mode(y) <- "double"
1173 /// mu <- t(apply(y, 1, function(r)
1174 /// rep(round(c(mean(r[1:3]), mean(r[4:6])) * 4) / 4, each = 3)))
1175 /// ave <- round(aveLogCPM(y, dispersion = 0.1) * 64) / 64
1176 /// dm <- edgeR:::.compressDispersions(y, 0.125)
1177 /// wm <- edgeR:::.compressWeights(y, NULL)
1178 /// .Call(edgeR:::.cxx_compute_ave_qd, y, mu, X, dm, ave, wm)
1179 /// # 18.590485410598912
1180 /// ```
1181 #[test]
1182 fn test_update_prior_matches_edger() {
1183 let (counts, fitted, ave) = prior_counts();
1184 let prior = update_prior(
1185 &counts,
1186 24,
1187 6,
1188 &fitted,
1189 &two_group(),
1190 2,
1191 &Recycled::scalar(0.125),
1192 None,
1193 &ave,
1194 )
1195 .unwrap();
1196
1197 assert_eq!(prior.len(), 1);
1198 assert_relative_eq!(prior[0], 18.590_485_410_598_912, max_relative = PRIOR_TOL);
1199 }
1200
1201 /// `update_prior` is exactly two rounds: the hand-run loop matches, one round
1202 /// does not.
1203 #[test]
1204 fn test_update_prior_is_exactly_two_rounds() {
1205 let (counts, fitted, ave) = prior_counts();
1206 let run = |prior: &[f64]| {
1207 let adjusted = compute_adjust_vec(
1208 &counts,
1209 24,
1210 6,
1211 &fitted,
1212 &two_group(),
1213 2,
1214 &Recycled::scalar(0.125),
1215 prior,
1216 None,
1217 )
1218 .unwrap();
1219 compute_prior(&ave, &adjusted.s2, &adjusted.df).unwrap()
1220 };
1221
1222 let first = run(&[1.0]);
1223 let second = run(&first);
1224 let got = update_prior(
1225 &counts,
1226 24,
1227 6,
1228 &fitted,
1229 &two_group(),
1230 2,
1231 &Recycled::scalar(0.125),
1232 None,
1233 &ave,
1234 )
1235 .unwrap();
1236
1237 assert_eq!(got, second);
1238 assert!(first[0] != second[0]);
1239 }
1240
1241 /// Dispersion `2^-20`, the near-Poisson corner: [`unit_nb_deviance`] takes
1242 /// its Poisson expansion while [`compute_weight`] uses the case 1 fit, as in
1243 /// edgeR. Exactly zero is covered by `test_zero_dispersion_matches_edger`.
1244 ///
1245 /// ```r
1246 /// d0 <- edgeR:::.compressDispersions(y, 2^-20)
1247 /// .Call(edgeR:::.cxx_compute_adj_vec, y, mu, X, d0, 1.0, wm)
1248 /// ```
1249 #[test]
1250 fn test_near_poisson_dispersion_matches_edger() {
1251 check(
1252 &two_group(),
1253 Recycled::scalar(2.0_f64.powi(-20)),
1254 &[1.0],
1255 None,
1256 &[
1257 0.360_893_151_979_798_77,
1258 0.122_110_248_576_896_45,
1259 12.329_306_966_107_474,
1260 0.396_987_455_932_764_3,
1261 0.0,
1262 11.686_733_868_487_11,
1263 ],
1264 &[
1265 3.991_092_303_143_825_7,
1266 3.999_478_990_707_440_4,
1267 4.976_310_706_206_089,
1268 1.997_955_076_334_263_3,
1269 0.0,
1270 4.000_035_625_925_281,
1271 ],
1272 &[
1273 0.090_424_656_852_841_87,
1274 0.030_531_538_948_101_14,
1275 2.477_599_911_663_728_6,
1276 0.198_696_887_950_621_38,
1277 0.0,
1278 2.921_657_445_434_316,
1279 ],
1280 );
1281 }
1282
1283 ////////////////////
1284 // Error branches //
1285 ////////////////////
1286
1287 /// Every rejection [`compute_adjust_vec`] can produce, in one table.
1288 #[test]
1289 fn test_rejects_bad_shapes() {
1290 let run = |n_genes: usize,
1291 n_samples: usize,
1292 counts: &[f64],
1293 fitted: &[f64],
1294 design: &[f64],
1295 n_coef: usize| {
1296 compute_adjust_vec(
1297 counts,
1298 n_genes,
1299 n_samples,
1300 fitted,
1301 design,
1302 n_coef,
1303 &Recycled::scalar(0.125),
1304 &[1.0],
1305 None,
1306 )
1307 .unwrap_err()
1308 };
1309
1310 assert!(matches!(
1311 run(0, 6, &[], &[], &two_group(), 2),
1312 EdgeErrors::EmptyCounts { .. }
1313 ));
1314 assert!(matches!(
1315 run(6, 0, &[], &[], &[], 2),
1316 EdgeErrors::EmptyCounts { .. }
1317 ));
1318 assert!(matches!(
1319 run(6, 6, &counts(), &fitted(), &[], 0),
1320 EdgeErrors::MustBePositive(_)
1321 ));
1322 assert!(matches!(
1323 run(5, 6, &counts(), &fitted(), &two_group(), 2),
1324 EdgeErrors::LengthMismatch { name: "counts", .. }
1325 ));
1326 assert!(matches!(
1327 run(6, 6, &counts(), &fitted()[..30], &two_group(), 2),
1328 EdgeErrors::LengthMismatch { name: "fitted", .. }
1329 ));
1330 assert!(matches!(
1331 run(6, 6, &counts(), &fitted(), &two_group(), 3),
1332 EdgeErrors::LengthMismatch { name: "design", .. }
1333 ));
1334 // Two samples, three coefficients: no residual df, and the QR would be
1335 // of a wide matrix.
1336 assert!(matches!(
1337 run(
1338 1,
1339 2,
1340 &[1.0, 2.0],
1341 &[1.0, 2.0],
1342 &[1.0, 0.0, 0.0, 1.0, 1.0, 1.0],
1343 3
1344 ),
1345 EdgeErrors::DesignNotFullRank { .. }
1346 ));
1347 }
1348
1349 #[test]
1350 fn test_rejects_mismatched_dispersion() {
1351 let err = compute_adjust_vec(
1352 &counts(),
1353 6,
1354 6,
1355 &fitted(),
1356 &two_group(),
1357 2,
1358 &Recycled::by_gene(vec![0.125; 3]),
1359 &[1.0],
1360 None,
1361 )
1362 .unwrap_err();
1363 assert!(matches!(
1364 err,
1365 EdgeErrors::LengthMismatch {
1366 name: "by_gene",
1367 ..
1368 }
1369 ));
1370 }
1371
1372 #[test]
1373 fn test_rejects_mismatched_weights() {
1374 let weights = Recycled::by_sample(vec![1.0; 3]);
1375 let err = compute_adjust_vec(
1376 &counts(),
1377 6,
1378 6,
1379 &fitted(),
1380 &two_group(),
1381 2,
1382 &Recycled::scalar(0.125),
1383 &[1.0],
1384 Some(&weights),
1385 )
1386 .unwrap_err();
1387 assert!(matches!(
1388 err,
1389 EdgeErrors::LengthMismatch {
1390 name: "by_sample",
1391 ..
1392 }
1393 ));
1394 }
1395
1396 #[test]
1397 fn test_rejects_wrong_prior_length() {
1398 let err = compute_adjust_vec(
1399 &counts(),
1400 6,
1401 6,
1402 &fitted(),
1403 &two_group(),
1404 2,
1405 &Recycled::scalar(0.125),
1406 &[1.0, 2.0],
1407 None,
1408 )
1409 .unwrap_err();
1410 assert!(matches!(
1411 err,
1412 EdgeErrors::LengthMismatch { name: "prior", .. }
1413 ));
1414 }
1415
1416 #[test]
1417 fn test_rejects_non_positive_prior() {
1418 for bad in [0.0, -1.0, f64::NAN, f64::INFINITY] {
1419 let err = compute_adjust_vec(
1420 &counts(),
1421 6,
1422 6,
1423 &fitted(),
1424 &two_group(),
1425 2,
1426 &Recycled::scalar(0.125),
1427 &[bad],
1428 None,
1429 )
1430 .unwrap_err();
1431 assert!(matches!(err, EdgeErrors::InvalidArgument(_)));
1432 }
1433 }
1434
1435 #[test]
1436 fn test_compute_prior_rejects_mismatched_lengths() {
1437 let err = compute_prior(&[1.0, 2.0], &[1.0], &[1.0, 1.0]).unwrap_err();
1438 assert!(matches!(err, EdgeErrors::LengthMismatch { name: "s2", .. }));
1439
1440 let err = compute_prior(&[1.0, 2.0], &[1.0, 1.0], &[1.0]).unwrap_err();
1441 assert!(matches!(err, EdgeErrors::LengthMismatch { name: "df", .. }));
1442 }
1443
1444 #[test]
1445 fn test_update_prior_rejects_mismatched_ave_log_cpm() {
1446 let err = update_prior(
1447 &counts(),
1448 6,
1449 6,
1450 &fitted(),
1451 &two_group(),
1452 2,
1453 &Recycled::scalar(0.125),
1454 None,
1455 &[0.0, 1.0],
1456 )
1457 .unwrap_err();
1458 assert!(matches!(
1459 err,
1460 EdgeErrors::LengthMismatch {
1461 name: "ave_log_cpm",
1462 ..
1463 }
1464 ));
1465 }
1466
1467 /// `phi = 0` must reach the case 1 Chebyshev fits, not the Poisson ones.
1468 /// edgeR hits this through `glmQLFit(dispersion = 0)` and any zero tagwise
1469 /// dispersion.
1470 ///
1471 /// ```r
1472 /// y <- matrix(c(8,16,32,64,128,256, 4,4,8,8,16,16, 64,64,64,128,128,128),
1473 /// nrow = 3, byrow = TRUE)
1474 /// X <- cbind(1, c(0,0,0,1,1,1))
1475 /// f <- glmQLFit(y, design = X, dispersion = 0,
1476 /// offset = matrix(log(1024), 3, 6), legacy = FALSE)
1477 /// f$deviance.adj # 147.92726615926262 11.90464907383647 0
1478 /// f$df.residual.adj # 5.3435958730008828 9.1582038968715125 3.9057308640743664
1479 /// ```
1480 #[test]
1481 fn test_zero_dispersion_matches_edger() {
1482 let counts: Vec<f64> = vec![
1483 8.0, 16.0, 32.0, 64.0, 128.0, 256.0, //
1484 4.0, 4.0, 8.0, 8.0, 16.0, 16.0, //
1485 64.0, 64.0, 64.0, 128.0, 128.0, 128.0, //
1486 ];
1487 let fitted: Vec<f64> = vec![
1488 18.666_666_666_666_668,
1489 18.666_666_666_666_668,
1490 18.666_666_666_666_668,
1491 149.333_333_333_333_4,
1492 149.333_333_333_333_4,
1493 149.333_333_333_333_4,
1494 5.333_333_333_333_334,
1495 5.333_333_333_333_334,
1496 5.333_333_333_333_334,
1497 13.333_333_333_333_334,
1498 13.333_333_333_333_334,
1499 13.333_333_333_333_334,
1500 64.000_000_000_000_03,
1501 64.000_000_000_000_03,
1502 64.000_000_000_000_03,
1503 128.000_000_000_000_09,
1504 128.000_000_000_000_09,
1505 128.000_000_000_000_09,
1506 ];
1507 let design = vec![1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0];
1508
1509 let out = compute_adjust_vec(
1510 &counts,
1511 3,
1512 6,
1513 &fitted,
1514 &design,
1515 2,
1516 &Recycled::scalar(0.0),
1517 &[15.741_184_694_250_476],
1518 None,
1519 )
1520 .unwrap();
1521
1522 let expected_deviance = [147.927_266_159_262_62, 11.904_649_073_836_47, 0.0];
1523 let expected_df = [
1524 5.343_595_873_000_883,
1525 9.158_203_896_871_512,
1526 3.905_730_864_074_366_4,
1527 ];
1528 for (got, want) in out.deviance.iter().zip(expected_deviance.iter()) {
1529 assert_relative_eq!(got, want, epsilon = 1e-12, max_relative = 1e-12);
1530 }
1531 for (got, want) in out.df.iter().zip(expected_df.iter()) {
1532 assert_relative_eq!(got, want, max_relative = 1e-12);
1533 }
1534 }
1535}