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

gam_inference/
effects.rs
