gam_inference/
effects.rs

1//! Matrix-level effect and contrast inference.
2//!
3//! This module owns the statistical kernel shared by difference-smooth,
4//! partial-dependence, and other linear-contrast reports.  Callers supply a
5//! coefficient vector, its covariance, and a contrast design.  Presentation
6//! layers only marshal the resulting typed report.
7
8use faer::Side;
9use gam_linalg::faer_ndarray::FaerEigh;
10use gam_math::probability::standard_normal_quantile;
11use gam_math::quantile::quantile_from_sorted;
12use gam_solve::estimate::UnifiedFitResult;
13use ndarray::{Array1, Array2, ArrayView1, ArrayView2};
14use rand::{RngExt, SeedableRng, rngs::StdRng};
15use std::error::Error;
16use std::fmt;
17
18/// Default confidence level for effect bands.
19pub const DEFAULT_BAND_LEVEL: f64 = 0.95;
20/// Default Monte Carlo draw count for simultaneous bands.
21pub const DEFAULT_SIMULATIONS: usize = 10_000;
22/// Default deterministic random seed for simultaneous bands.
23pub const DEFAULT_SIMULATION_SEED: u64 = 12_345;
24
25/// The coefficient covariance definition used by an effect report.
26#[derive(Clone, Copy, Debug, Eq, PartialEq)]
27pub enum CovarianceSource {
28    /// Conditional Bayesian covariance with smoothing parameters fixed (`Vb`).
29    Conditional,
30    /// Bayesian covariance including smoothing-parameter uncertainty (`Vp`).
31    SmoothingCorrected,
32    /// Frequentist sandwich covariance (`Ve`).
33    Frequentist,
34}
35
36impl fmt::Display for CovarianceSource {
37    fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
38        formatter.write_str(match self {
39            Self::Conditional => "conditional",
40            Self::SmoothingCorrected => "smoothing-corrected",
41            Self::Frequentist => "frequentist",
42        })
43    }
44}
45
46/// A covariance borrowed from a fit together with its exact provenance.
47#[derive(Clone, Copy, Debug)]
48pub struct SelectedCovariance<'a> {
49    pub source: CovarianceSource,
50    pub matrix: ArrayView2<'a, f64>,
51}
52
53/// Select a coefficient covariance from a unified fit under an explicit policy.
54///
55/// The requested source is exact: absence is an error and no other covariance
56/// definition is substituted. [`SelectedCovariance::source`] records the
57/// resolved provenance.
58pub fn select_covariance<'a>(
59    fit: &'a UnifiedFitResult,
60    source: CovarianceSource,
61) -> Result<SelectedCovariance<'a>, EffectError> {
62    let matrix =
63        covariance_by_source(fit, source).ok_or(EffectError::MissingCovariance { source })?;
64
65    Ok(SelectedCovariance {
66        source,
67        matrix: matrix.view(),
68    })
69}
70
71fn covariance_by_source(fit: &UnifiedFitResult, source: CovarianceSource) -> Option<&Array2<f64>> {
72    match source {
73        CovarianceSource::Conditional => fit.beta_covariance(),
74        CovarianceSource::SmoothingCorrected => fit.beta_covariance_corrected(),
75        CovarianceSource::Frequentist => fit.beta_covariance_ve(),
76    }
77}
78
79/// Configuration for a pointwise normal-theory confidence band.
80#[derive(Clone, Copy, Debug, PartialEq)]
81pub struct PointwiseBandOptions {
82    pub level: f64,
83}
84
85impl Default for PointwiseBandOptions {
86    fn default() -> Self {
87        Self {
88            level: DEFAULT_BAND_LEVEL,
89        }
90    }
91}
92
93/// Configuration for a simulated simultaneous confidence band.
94#[derive(Clone, Copy, Debug, PartialEq)]
95pub struct SimultaneousBandOptions {
96    pub level: f64,
97    pub simulations: usize,
98    pub seed: u64,
99}
100
101impl Default for SimultaneousBandOptions {
102    fn default() -> Self {
103        Self {
104            level: DEFAULT_BAND_LEVEL,
105            simulations: DEFAULT_SIMULATIONS,
106            seed: DEFAULT_SIMULATION_SEED,
107        }
108    }
109}
110
111/// Confidence-band procedure for a linear effect curve.
112#[derive(Clone, Copy, Debug, PartialEq)]
113pub enum BandOptions {
114    /// Independent marginal normal intervals at each contrast row.
115    Pointwise(PointwiseBandOptions),
116    /// A common critical value calibrated from the supremum of the standardized
117    /// Gaussian effect curve.
118    Simultaneous(SimultaneousBandOptions),
119}
120
121impl Default for BandOptions {
122    fn default() -> Self {
123        Self::Pointwise(PointwiseBandOptions::default())
124    }
125}
126
127/// A matrix-level effect report, with one entry per contrast-design row.
128#[derive(Clone, Debug, PartialEq)]
129pub struct EffectReport {
130    pub center: Array1<f64>,
131    pub se: Array1<f64>,
132    pub lower: Array1<f64>,
133    pub upper: Array1<f64>,
134    pub critical: f64,
135}
136
137/// Typed failures from covariance selection or effect-band construction.
138#[derive(Clone, Debug, PartialEq)]
139pub enum EffectError {
140    MissingCovariance {
141        source: CovarianceSource,
142    },
143    EmptyCoefficients,
144    EmptyContrastDesign,
145    InvalidLevel {
146        level: f64,
147    },
148    InvalidSimulationCount,
149    CovarianceShape {
150        rows: usize,
151        columns: usize,
152        expected: usize,
153    },
154    ContrastShape {
155        columns: usize,
156        expected: usize,
157    },
158    NonFiniteInput {
159        input: &'static str,
160    },
161    NonSymmetricCovariance {
162        row: usize,
163        column: usize,
164        difference: f64,
165        tolerance: f64,
166    },
167    IndefiniteCovariance {
168        matrix: &'static str,
169        minimum_eigenvalue: f64,
170        tolerance: f64,
171    },
172    Eigendecomposition {
173        matrix: &'static str,
174        detail: String,
175    },
176    NormalQuantile {
177        detail: String,
178    },
179}
180
181impl fmt::Display for EffectError {
182    fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
183        match self {
184            Self::MissingCovariance { source } => {
185                write!(formatter, "fit has no {source} coefficient covariance")
186            }
187            Self::EmptyCoefficients => {
188                formatter.write_str("beta must contain at least one coefficient")
189            }
190            Self::EmptyContrastDesign => {
191                formatter.write_str("contrast design must contain at least one row")
192            }
193            Self::InvalidLevel { level } => {
194                write!(
195                    formatter,
196                    "confidence level must be finite and in (0, 1), got {level}"
197                )
198            }
199            Self::InvalidSimulationCount => {
200                formatter.write_str("simultaneous-band simulation count must be positive")
201            }
202            Self::CovarianceShape {
203                rows,
204                columns,
205                expected,
206            } => write!(
207                formatter,
208                "covariance must have shape {expected}x{expected}, got {rows}x{columns}"
209            ),
210            Self::ContrastShape { columns, expected } => write!(
211                formatter,
212                "contrast design must have {expected} columns, got {columns}"
213            ),
214            Self::NonFiniteInput { input } => {
215                write!(formatter, "{input} contains a non-finite value")
216            }
217            Self::NonSymmetricCovariance {
218                row,
219                column,
220                difference,
221                tolerance,
222            } => write!(
223                formatter,
224                "covariance is not symmetric at ({row}, {column}): absolute difference {difference:e} exceeds tolerance {tolerance:e}"
225            ),
226            Self::IndefiniteCovariance {
227                matrix,
228                minimum_eigenvalue,
229                tolerance,
230            } => write!(
231                formatter,
232                "{matrix} is materially indefinite: minimum eigenvalue {minimum_eigenvalue:e} is below -{tolerance:e}"
233            ),
234            Self::Eigendecomposition { matrix, detail } => {
235                write!(formatter, "{matrix} eigendecomposition failed: {detail}")
236            }
237            Self::NormalQuantile { detail } => {
238                write!(
239                    formatter,
240                    "normal critical-value calculation failed: {detail}"
241                )
242            }
243        }
244    }
245}
246
247impl Error for EffectError {}
248
249/// Construct an effect report directly from a unified fit and an explicit
250/// covariance-selection policy.
251pub fn effect_report_from_fit(
252    fit: &UnifiedFitResult,
253    source: CovarianceSource,
254    contrast_design: ArrayView2<'_, f64>,
255    options: BandOptions,
256) -> Result<EffectReport, EffectError> {
257    let selected = select_covariance(fit, source)?;
258    effect_report(fit.beta.view(), selected.matrix, contrast_design, options)
259}
260
261/// Compute centers, standard errors, bounds, and the common critical value for
262/// a linear effect curve.
263///
264/// For `m x p` contrast design `C`, coefficient vector `beta`, and covariance
265/// `V`, the report center is `C beta` and its covariance is `C V C'`.
266/// Pointwise bands compute only the diagonal of that covariance, with O(p)
267/// working memory. Simultaneous bands calibrate `max_i |Z_i|` for the
268/// standardized Gaussian curve and factor whichever covariance space is
269/// smaller: coefficient space when `p <= m`, projected curve space otherwise.
270/// Positive-semidefinite singular matrices are supported without a ridge.
271pub fn effect_report(
272    beta: ArrayView1<'_, f64>,
273    covariance: ArrayView2<'_, f64>,
274    contrast_design: ArrayView2<'_, f64>,
275    options: BandOptions,
276) -> Result<EffectReport, EffectError> {
277    validate_inputs(beta, covariance, contrast_design, options)?;
278
279    let covariance = validated_symmetric_matrix(covariance)?;
280    let center = contrast_design.dot(&beta);
281    let (se, critical) = match options {
282        BandOptions::Pointwise(pointwise) => {
283            let se = pointwise_standard_errors(contrast_design, covariance.view())?;
284            let critical = standard_normal_quantile(0.5 * (1.0 + pointwise.level))
285                .map_err(|detail| EffectError::NormalQuantile { detail })?;
286            (se, critical)
287        }
288        BandOptions::Simultaneous(simultaneous) => {
289            let curve_factor = simultaneous_curve_factor(contrast_design, covariance.view())?;
290            let se = factor_standard_errors(&curve_factor);
291            let critical = simultaneous_critical(
292                &curve_factor,
293                se.view(),
294                simultaneous.level,
295                simultaneous.simulations,
296                simultaneous.seed,
297            );
298            (se, critical)
299        }
300    };
301
302    let half_width = se.mapv(|value| critical * value);
303    let lower = &center - &half_width;
304    let upper = &center + &half_width;
305    Ok(EffectReport {
306        center,
307        se,
308        lower,
309        upper,
310        critical,
311    })
312}
313
314fn pointwise_standard_errors(
315    contrast_design: ArrayView2<'_, f64>,
316    covariance: ArrayView2<'_, f64>,
317) -> Result<Array1<f64>, EffectError> {
318    let p = covariance.nrows();
319    let mut se = Array1::<f64>::zeros(contrast_design.nrows());
320    let mut product = vec![0.0_f64; p];
321    for (row_index, row) in contrast_design.rows().into_iter().enumerate() {
322        product.fill(0.0);
323        for covariance_row in 0..p {
324            for column in 0..p {
325                product[covariance_row] += covariance[[covariance_row, column]] * row[column];
326            }
327        }
328        let variance = row
329            .iter()
330            .zip(&product)
331            .map(|(&loading, &projected)| loading * projected)
332            .sum::<f64>();
333        let scale = row
334            .iter()
335            .zip(&product)
336            .map(|(&loading, &projected)| (loading * projected).abs())
337            .sum::<f64>();
338        let tolerance = roundoff_tolerance(scale, p);
339        if variance < -tolerance {
340            return Err(EffectError::IndefiniteCovariance {
341                matrix: "projected curve covariance",
342                minimum_eigenvalue: variance,
343                tolerance,
344            });
345        }
346        se[row_index] = variance.max(0.0).sqrt();
347    }
348    Ok(se)
349}
350
351fn simultaneous_curve_factor(
352    contrast_design: ArrayView2<'_, f64>,
353    covariance: ArrayView2<'_, f64>,
354) -> Result<Array2<f64>, EffectError> {
355    if covariance.nrows() <= contrast_design.nrows() {
356        let coefficient_eigen = psd_eigendecomposition(covariance, "coefficient covariance")?;
357        let coefficient_factor = covariance_factor(&coefficient_eigen);
358        let factor = contrast_design.dot(&coefficient_factor);
359        if factor.iter().any(|value| !value.is_finite()) {
360            return Err(EffectError::NonFiniteInput {
361                input: "projected curve factor",
362            });
363        }
364        return Ok(factor);
365    }
366
367    let mut curve_covariance = contrast_design.dot(&covariance).dot(&contrast_design.t());
368    if curve_covariance.iter().any(|value| !value.is_finite()) {
369        return Err(EffectError::NonFiniteInput {
370            input: "projected curve covariance",
371        });
372    }
373    symmetrize_in_place(&mut curve_covariance);
374    let curve_eigen =
375        psd_eigendecomposition(curve_covariance.view(), "projected curve covariance")?;
376    Ok(covariance_factor(&curve_eigen))
377}
378
379fn validate_inputs(
380    beta: ArrayView1<'_, f64>,
381    covariance: ArrayView2<'_, f64>,
382    contrast_design: ArrayView2<'_, f64>,
383    options: BandOptions,
384) -> Result<(), EffectError> {
385    if beta.is_empty() {
386        return Err(EffectError::EmptyCoefficients);
387    }
388    if contrast_design.nrows() == 0 {
389        return Err(EffectError::EmptyContrastDesign);
390    }
391    let p = beta.len();
392    if covariance.dim() != (p, p) {
393        return Err(EffectError::CovarianceShape {
394            rows: covariance.nrows(),
395            columns: covariance.ncols(),
396            expected: p,
397        });
398    }
399    if contrast_design.ncols() != p {
400        return Err(EffectError::ContrastShape {
401            columns: contrast_design.ncols(),
402            expected: p,
403        });
404    }
405    if beta.iter().any(|value| !value.is_finite()) {
406        return Err(EffectError::NonFiniteInput { input: "beta" });
407    }
408    if covariance.iter().any(|value| !value.is_finite()) {
409        return Err(EffectError::NonFiniteInput {
410            input: "covariance",
411        });
412    }
413    if contrast_design.iter().any(|value| !value.is_finite()) {
414        return Err(EffectError::NonFiniteInput {
415            input: "contrast design",
416        });
417    }
418
419    let (level, simulations) = match options {
420        BandOptions::Pointwise(pointwise) => (pointwise.level, None),
421        BandOptions::Simultaneous(simultaneous) => {
422            (simultaneous.level, Some(simultaneous.simulations))
423        }
424    };
425    if !level.is_finite() || !(0.0..1.0).contains(&level) || level == 0.0 {
426        return Err(EffectError::InvalidLevel { level });
427    }
428    if simulations == Some(0) {
429        return Err(EffectError::InvalidSimulationCount);
430    }
431    Ok(())
432}
433
434struct PsdEigen {
435    vectors: Array2<f64>,
436    active: Vec<(usize, f64)>,
437}
438
439fn validated_symmetric_matrix(matrix: ArrayView2<'_, f64>) -> Result<Array2<f64>, EffectError> {
440    let n = matrix.nrows();
441    let scale = matrix
442        .iter()
443        .fold(0.0_f64, |maximum, value| maximum.max(value.abs()));
444    let symmetry_tolerance = roundoff_tolerance(scale, n);
445    let mut symmetric = matrix.to_owned();
446    for row in 0..n {
447        for column in 0..row {
448            let difference = (matrix[[row, column]] - matrix[[column, row]]).abs();
449            if difference > symmetry_tolerance {
450                return Err(EffectError::NonSymmetricCovariance {
451                    row,
452                    column,
453                    difference,
454                    tolerance: symmetry_tolerance,
455                });
456            }
457            let average = 0.5 * (matrix[[row, column]] + matrix[[column, row]]);
458            symmetric[[row, column]] = average;
459            symmetric[[column, row]] = average;
460        }
461    }
462    Ok(symmetric)
463}
464
465fn symmetrize_in_place(matrix: &mut Array2<f64>) {
466    for row in 0..matrix.nrows() {
467        for column in 0..row {
468            let average = 0.5 * (matrix[[row, column]] + matrix[[column, row]]);
469            matrix[[row, column]] = average;
470            matrix[[column, row]] = average;
471        }
472    }
473}
474
475fn psd_eigendecomposition(
476    matrix: ArrayView2<'_, f64>,
477    label: &'static str,
478) -> Result<PsdEigen, EffectError> {
479    let n = matrix.nrows();
480    let symmetric = validated_symmetric_matrix(matrix)?;
481
482    let (values, vectors) =
483        symmetric
484            .eigh(Side::Lower)
485            .map_err(|error| EffectError::Eigendecomposition {
486                matrix: label,
487                detail: error.to_string(),
488            })?;
489    let spectral_scale = values
490        .iter()
491        .fold(0.0_f64, |maximum, value| maximum.max(value.abs()));
492    let tolerance = roundoff_tolerance(spectral_scale, n);
493    let minimum_eigenvalue = values.iter().copied().fold(f64::INFINITY, f64::min);
494    if minimum_eigenvalue < -tolerance {
495        return Err(EffectError::IndefiniteCovariance {
496            matrix: label,
497            minimum_eigenvalue,
498            tolerance,
499        });
500    }
501    let active = values
502        .iter()
503        .copied()
504        .enumerate()
505        .filter_map(|(column, value)| (value > tolerance).then(|| (column, value.sqrt())))
506        .collect();
507    Ok(PsdEigen { vectors, active })
508}
509
510fn roundoff_tolerance(scale: f64, dimension: usize) -> f64 {
511    scale * f64::EPSILON * dimension.max(1) as f64
512}
513
514fn covariance_factor(eigen: &PsdEigen) -> Array2<f64> {
515    let mut factor = Array2::zeros((eigen.vectors.nrows(), eigen.active.len()));
516    for (active_column, &(eigen_column, eigenvalue_sqrt)) in eigen.active.iter().enumerate() {
517        for row in 0..eigen.vectors.nrows() {
518            factor[[row, active_column]] = eigen.vectors[[row, eigen_column]] * eigenvalue_sqrt;
519        }
520    }
521    factor
522}
523
524fn factor_standard_errors(curve_factor: &Array2<f64>) -> Array1<f64> {
525    let variances = Array1::from_iter(
526        curve_factor
527            .rows()
528            .into_iter()
529            .map(|row| row.iter().map(|value| value * value).sum::<f64>()),
530    );
531    let variance_scale = variances.iter().copied().fold(0.0_f64, f64::max);
532    let variance_tolerance = roundoff_tolerance(variance_scale, curve_factor.ncols());
533    variances.mapv(|variance| {
534        if variance > variance_tolerance {
535            variance.sqrt()
536        } else {
537            0.0
538        }
539    })
540}
541
542fn simultaneous_critical(
543    curve_factor: &Array2<f64>,
544    se: ArrayView1<'_, f64>,
545    level: f64,
546    simulations: usize,
547    seed: u64,
548) -> f64 {
549    if curve_factor.ncols() == 0 {
550        return 0.0;
551    }
552
553    let mut standardized_factor = curve_factor.clone();
554    for row in 0..standardized_factor.nrows() {
555        if se[row] == 0.0 {
556            standardized_factor.row_mut(row).fill(0.0);
557        } else {
558            standardized_factor
559                .row_mut(row)
560                .mapv_inplace(|value| value / se[row]);
561        }
562    }
563
564    let mut rng = StdRng::seed_from_u64(seed);
565    let mut normal_coordinates = vec![0.0; curve_factor.ncols()];
566    let mut maxima = Vec::with_capacity(simulations);
567    for _ in 0..simulations {
568        fill_standard_normals(&mut rng, &mut normal_coordinates);
569        let maximum = standardized_factor
570            .rows()
571            .into_iter()
572            .map(|row| {
573                row.iter()
574                    .zip(&normal_coordinates)
575                    .map(|(&loading, &coordinate)| loading * coordinate)
576                    .sum::<f64>()
577                    .abs()
578            })
579            .fold(0.0_f64, f64::max);
580        maxima.push(maximum);
581    }
582    maxima.sort_by(f64::total_cmp);
583    quantile_from_sorted(&maxima, level)
584}
585
586fn fill_standard_normals(rng: &mut StdRng, output: &mut [f64]) {
587    for pair in output.chunks_mut(2) {
588        let uniform_radius = rng.random::<f64>().max(f64::MIN_POSITIVE);
589        let uniform_angle = rng.random::<f64>();
590        let radius = (-2.0 * uniform_radius.ln()).sqrt();
591        let angle = std::f64::consts::TAU * uniform_angle;
592        pair[0] = radius * angle.cos();
593        if pair.len() == 2 {
594            pair[1] = radius * angle.sin();
595        }
596    }
597}
598
599#[cfg(test)]
600mod tests {
601    use super::*;
602    use approx::assert_abs_diff_eq;
603    use ndarray::array;
604
605    #[test]
606    fn closed_form_centers_and_standard_errors() {
607        let beta = array![2.0, -1.0];
608        let covariance = array![[4.0, 1.0], [1.0, 9.0]];
609        let contrast = array![[1.0, 0.0], [1.0, 2.0]];
610
611        let report = effect_report(
612            beta.view(),
613            covariance.view(),
614            contrast.view(),
615            BandOptions::default(),
616        )
617        .unwrap();
618
619        assert_abs_diff_eq!(report.center[0], 2.0, epsilon = 1e-14);
620        assert_abs_diff_eq!(report.center[1], 0.0, epsilon = 1e-14);
621        assert_abs_diff_eq!(report.se[0], 2.0, epsilon = 1e-14);
622        assert_abs_diff_eq!(report.se[1], 44.0_f64.sqrt(), epsilon = 1e-13);
623    }
624
625    #[test]
626    fn singular_psd_simulation_is_reproducible() {
627        let beta = array![0.5, -0.5];
628        let covariance = array![[1.0, 1.0], [1.0, 1.0]];
629        let contrast = array![[1.0, 0.0], [0.0, 1.0], [1.0, -1.0]];
630        let options = BandOptions::Simultaneous(SimultaneousBandOptions {
631            simulations: 2_000,
632            ..SimultaneousBandOptions::default()
633        });
634
635        let first =
636            effect_report(beta.view(), covariance.view(), contrast.view(), options).unwrap();
637        let second =
638            effect_report(beta.view(), covariance.view(), contrast.view(), options).unwrap();
639
640        assert_eq!(first, second);
641        assert_abs_diff_eq!(first.se[0], 1.0, epsilon = 1e-14);
642        assert_abs_diff_eq!(first.se[1], 1.0, epsilon = 1e-14);
643        assert_eq!(first.se[2], 0.0);
644        assert!(first.critical.is_finite());
645    }
646
647    #[test]
648    fn materially_indefinite_covariance_is_rejected() {
649        let error = effect_report(
650            array![0.0, 0.0].view(),
651            array![[1.0, 0.0], [0.0, -0.1]].view(),
652            array![[0.0, 1.0]].view(),
653            BandOptions::default(),
654        )
655        .unwrap_err();
656
657        assert!(matches!(
658            error,
659            EffectError::IndefiniteCovariance {
660                matrix: "projected curve covariance",
661                ..
662            }
663        ));
664    }
665
666    #[test]
667    fn pointwise_band_uses_two_sided_normal_critical_value() {
668        let report = effect_report(
669            array![0.0].view(),
670            array![[1.0]].view(),
671            array![[1.0]].view(),
672            BandOptions::default(),
673        )
674        .unwrap();
675
676        // True two-sided 95% normal critical value z_{0.975} = 1.9599639845400545.
677        // The former golden 1.959963986120195 was the RAW Acklam-approximation
678        // output (absolute error ~1.6e-9); the quantile now carries a two-round
679        // Halley refinement against erfc and returns the true value, so the
680        // golden pins the mathematically correct quantile at a tightened
681        // tolerance.
682        assert_abs_diff_eq!(report.critical, 1.959_963_984_540_054, epsilon = 1e-9);
683    }
684
685    #[test]
686    fn simultaneous_critical_is_beta_and_contrast_sign_invariant() {
687        let covariance = array![[2.0, 0.25], [0.25, 1.0]];
688        let contrast = array![[1.0, 0.5], [-0.25, 1.0]];
689        let options = BandOptions::Simultaneous(SimultaneousBandOptions {
690            simulations: 1_000,
691            ..SimultaneousBandOptions::default()
692        });
693        let first = effect_report(
694            array![1.0, -2.0].view(),
695            covariance.view(),
696            contrast.view(),
697            options,
698        )
699        .unwrap();
700        let shifted = effect_report(
701            array![8.0, 3.0].view(),
702            covariance.view(),
703            contrast.view(),
704            options,
705        )
706        .unwrap();
707        let signed = effect_report(
708            array![1.0, -2.0].view(),
709            covariance.view(),
710            (-&contrast).view(),
711            options,
712        )
713        .unwrap();
714
715        assert_eq!(first.critical, shifted.critical);
716        assert_eq!(first.critical, signed.critical);
717        assert_eq!(first.se, shifted.se);
718        assert_eq!(first.se, signed.se);
719        for row in 0..contrast.nrows() {
720            assert_abs_diff_eq!(signed.center[row], -first.center[row], epsilon = 1e-14);
721            assert_abs_diff_eq!(signed.lower[row], -first.upper[row], epsilon = 1e-14);
722            assert_abs_diff_eq!(signed.upper[row], -first.lower[row], epsilon = 1e-14);
723        }
724    }
725}
gam_inference/effects.rs

gam_inference/
effects.rs
