fdars_core/explain/
importance.rs

1//! Permutation importance, pointwise importance, and conditional permutation importance.
2
3use super::helpers::{
4    clone_scores_matrix, compute_conditioning_bins, compute_score_variance,
5    logistic_accuracy_from_scores, permute_component, project_scores, shuffle_global,
6};
7use crate::error::FdarError;
8use crate::matrix::FdMatrix;
9use crate::scalar_on_function::{sigmoid, FregreLmResult, FunctionalLogisticResult};
10use rand::prelude::*;
11
12// ===========================================================================
13// FPC Permutation Importance
14// ===========================================================================
15
16/// Result of FPC permutation importance.
17#[derive(Debug, Clone, PartialEq)]
18pub struct FpcPermutationImportance {
19    /// R^2 (or accuracy) drop per component (length ncomp).
20    pub importance: Vec<f64>,
21    /// Baseline metric (R^2 or accuracy).
22    pub baseline_metric: f64,
23    /// Mean metric after permuting each component.
24    pub permuted_metric: Vec<f64>,
25}
26
27/// Permutation importance for a linear functional regression (metric = R^2).
28///
29/// # Errors
30///
31/// Returns [`FdarError::InvalidDimension`] if `data` has zero rows, its column
32/// count does not match `fit.fpca.mean`, or `y.len()` does not match the row
33/// count.
34/// Returns [`FdarError::InvalidParameter`] if `n_perm` is zero.
35/// Returns [`FdarError::ComputationFailed`] if the total sum of squares is zero.
36///
37/// # Examples
38///
39/// ```
40/// use fdars_core::matrix::FdMatrix;
41/// use fdars_core::scalar_on_function::fregre_lm;
42/// use fdars_core::explain::fpc_permutation_importance;
43///
44/// let (n, m) = (20, 30);
45/// let data = FdMatrix::from_column_major(
46///     (0..n * m).map(|k| {
47///         let i = (k % n) as f64;
48///         let j = (k / n) as f64;
49///         ((i + 1.0) * j * 0.2).sin()
50///     }).collect(),
51///     n, m,
52/// ).unwrap();
53/// let y: Vec<f64> = (0..n).map(|i| (i as f64 * 0.5).sin()).collect();
54/// let fit = fregre_lm(&data, &y, None, 3).unwrap();
55/// let imp = fpc_permutation_importance(&fit, &data, &y, 10, 42).unwrap();
56/// assert_eq!(imp.importance.len(), 3);
57/// ```
58#[must_use = "expensive computation whose result should not be discarded"]
59pub fn fpc_permutation_importance(
60    fit: &FregreLmResult,
61    data: &FdMatrix,
62    y: &[f64],
63    n_perm: usize,
64    seed: u64,
65) -> Result<FpcPermutationImportance, FdarError> {
66    let (n, m) = data.shape();
67    if n == 0 {
68        return Err(FdarError::InvalidDimension {
69            parameter: "data",
70            expected: ">0 rows".into(),
71            actual: "0".into(),
72        });
73    }
74    if n != y.len() {
75        return Err(FdarError::InvalidDimension {
76            parameter: "y",
77            expected: format!("{n} (matching data rows)"),
78            actual: format!("{}", y.len()),
79        });
80    }
81    if m != fit.fpca.mean.len() {
82        return Err(FdarError::InvalidDimension {
83            parameter: "data",
84            expected: format!("{} columns", fit.fpca.mean.len()),
85            actual: format!("{m}"),
86        });
87    }
88    if n_perm == 0 {
89        return Err(FdarError::InvalidParameter {
90            parameter: "n_perm",
91            message: "must be > 0".into(),
92        });
93    }
94    let ncomp = fit.ncomp;
95    let scores = project_scores(
96        data,
97        &fit.fpca.mean,
98        &fit.fpca.rotation,
99        ncomp,
100        &fit.fpca.weights,
101    );
102
103    // Baseline R^2 -- compute from same FPC-only prediction used in permuted path
104    // to ensure consistent comparison (gamma terms are constant across permutations)
105    let y_mean: f64 = y.iter().sum::<f64>() / n as f64;
106    let ss_tot: f64 = y.iter().map(|&yi| (yi - y_mean).powi(2)).sum();
107    if ss_tot == 0.0 {
108        return Err(FdarError::ComputationFailed {
109            operation: "fpc_permutation_importance",
110            detail: "total sum of squares is zero; all response values may be identical — check your data".into(),
111        });
112    }
113    let identity_idx: Vec<usize> = (0..n).collect();
114    let ss_res_base = permuted_ss_res_linear(
115        &scores,
116        &fit.coefficients,
117        y,
118        n,
119        ncomp,
120        ncomp,
121        &identity_idx,
122    );
123    let baseline = 1.0 - ss_res_base / ss_tot;
124
125    // NOT migrated to permutation_test::permutation_pvalue — uses a single ADVANCING StdRng across ALL
126    // components; per-perm reseed would change the importances (Phase-49 CONS-02 Plan A). The Phase-48
127    // hand-off "fold into the already-parallel generic path" is behavior-CHANGING (generic path reseeds
128    // per-component seed+k; this advances one RNG) → DEFERRED to a future behavior-changing phase with a
129    // re-baselined golden, NOT done here.
130    let mut rng = StdRng::seed_from_u64(seed);
131    let mut importance = vec![0.0; ncomp];
132    let mut permuted_metric = vec![0.0; ncomp];
133
134    for k in 0..ncomp {
135        let mut sum_r2 = 0.0;
136        for _ in 0..n_perm {
137            let mut idx: Vec<usize> = (0..n).collect();
138            idx.shuffle(&mut rng);
139            let ss_res_perm =
140                permuted_ss_res_linear(&scores, &fit.coefficients, y, n, ncomp, k, &idx);
141            sum_r2 += 1.0 - ss_res_perm / ss_tot;
142        }
143        let mean_perm = sum_r2 / n_perm as f64;
144        permuted_metric[k] = mean_perm;
145        importance[k] = baseline - mean_perm;
146    }
147
148    Ok(FpcPermutationImportance {
149        importance,
150        baseline_metric: baseline,
151        permuted_metric,
152    })
153}
154
155/// Permutation importance for functional logistic regression (metric = accuracy).
156///
157/// # Errors
158///
159/// Returns [`FdarError::InvalidDimension`] if `data` has zero rows, its column
160/// count does not match `fit.fpca.mean`, or `y.len()` does not match the row
161/// count.
162/// Returns [`FdarError::InvalidParameter`] if `n_perm` is zero.
163#[must_use = "expensive computation whose result should not be discarded"]
164pub fn fpc_permutation_importance_logistic(
165    fit: &FunctionalLogisticResult,
166    data: &FdMatrix,
167    y: &[f64],
168    n_perm: usize,
169    seed: u64,
170) -> Result<FpcPermutationImportance, FdarError> {
171    let (n, m) = data.shape();
172    if n == 0 {
173        return Err(FdarError::InvalidDimension {
174            parameter: "data",
175            expected: ">0 rows".into(),
176            actual: "0".into(),
177        });
178    }
179    if n != y.len() {
180        return Err(FdarError::InvalidDimension {
181            parameter: "y",
182            expected: format!("{n} (matching data rows)"),
183            actual: format!("{}", y.len()),
184        });
185    }
186    if m != fit.fpca.mean.len() {
187        return Err(FdarError::InvalidDimension {
188            parameter: "data",
189            expected: format!("{} columns", fit.fpca.mean.len()),
190            actual: format!("{m}"),
191        });
192    }
193    if n_perm == 0 {
194        return Err(FdarError::InvalidParameter {
195            parameter: "n_perm",
196            message: "must be > 0".into(),
197        });
198    }
199    let ncomp = fit.ncomp;
200    let scores = project_scores(
201        data,
202        &fit.fpca.mean,
203        &fit.fpca.rotation,
204        ncomp,
205        &fit.fpca.weights,
206    );
207
208    let baseline: f64 = (0..n)
209        .filter(|&i| {
210            let pred = if fit.probabilities[i] >= 0.5 {
211                1.0
212            } else {
213                0.0
214            };
215            (pred - y[i]).abs() < 1e-10
216        })
217        .count() as f64
218        / n as f64;
219
220    // NOT migrated to permutation_test::permutation_pvalue — uses a single ADVANCING StdRng across ALL
221    // components; per-perm reseed would change the importances (Phase-49 CONS-02 Plan A). The Phase-48
222    // hand-off "fold into the already-parallel generic path" is behavior-CHANGING (generic path reseeds
223    // per-component seed+k; this advances one RNG) → DEFERRED to a future behavior-changing phase with a
224    // re-baselined golden, NOT done here.
225    let mut rng = StdRng::seed_from_u64(seed);
226    let mut importance = vec![0.0; ncomp];
227    let mut permuted_metric = vec![0.0; ncomp];
228
229    for k in 0..ncomp {
230        let mut sum_acc = 0.0;
231        for _ in 0..n_perm {
232            let mut perm_scores = clone_scores_matrix(&scores, n, ncomp);
233            shuffle_global(&mut perm_scores, &scores, k, n, &mut rng);
234            sum_acc += logistic_accuracy_from_scores(
235                &perm_scores,
236                fit.intercept,
237                &fit.coefficients,
238                y,
239                n,
240                ncomp,
241            );
242        }
243        let mean_acc = sum_acc / n_perm as f64;
244        permuted_metric[k] = mean_acc;
245        importance[k] = baseline - mean_acc;
246    }
247
248    Ok(FpcPermutationImportance {
249        importance,
250        baseline_metric: baseline,
251        permuted_metric,
252    })
253}
254
255/// Compute SS_res with component k shuffled by given index permutation.
256fn permuted_ss_res_linear(
257    scores: &FdMatrix,
258    coefficients: &[f64],
259    y: &[f64],
260    n: usize,
261    ncomp: usize,
262    k: usize,
263    perm_idx: &[usize],
264) -> f64 {
265    (0..n)
266        .map(|i| {
267            let mut yhat = coefficients[0];
268            for c in 0..ncomp {
269                let s = if c == k {
270                    scores[(perm_idx[i], c)]
271                } else {
272                    scores[(i, c)]
273                };
274                yhat += coefficients[1 + c] * s;
275            }
276            (y[i] - yhat).powi(2)
277        })
278        .sum()
279}
280
281// ===========================================================================
282// Pointwise Variable Importance
283// ===========================================================================
284
285/// Result of pointwise variable importance analysis.
286#[derive(Debug, Clone, PartialEq)]
287#[non_exhaustive]
288pub struct PointwiseImportanceResult {
289    /// Importance at each grid point (length m).
290    pub importance: Vec<f64>,
291    /// Normalized importance summing to 1 (length m).
292    pub importance_normalized: Vec<f64>,
293    /// Per-component importance (ncomp x m).
294    pub component_importance: FdMatrix,
295    /// Variance of each FPC score (length ncomp).
296    pub score_variance: Vec<f64>,
297}
298
299/// Pointwise variable importance for a linear functional regression model.
300///
301/// Measures how much X(t_j) contributes to prediction variance via the FPC decomposition.
302///
303/// # Errors
304///
305/// Returns [`FdarError::InvalidParameter`] if `fit.ncomp` is zero.
306/// Returns [`FdarError::InvalidDimension`] if the rotation matrix has zero rows
307/// or the scores matrix has fewer than 2 rows.
308#[must_use = "expensive computation whose result should not be discarded"]
309pub fn pointwise_importance(fit: &FregreLmResult) -> Result<PointwiseImportanceResult, FdarError> {
310    let ncomp = fit.ncomp;
311    let m = fit.fpca.rotation.nrows();
312    let n = fit.fpca.scores.nrows();
313    if ncomp == 0 {
314        return Err(FdarError::InvalidParameter {
315            parameter: "ncomp",
316            message: "must be > 0".into(),
317        });
318    }
319    if m == 0 {
320        return Err(FdarError::InvalidDimension {
321            parameter: "rotation",
322            expected: ">0 rows".into(),
323            actual: "0".into(),
324        });
325    }
326    if n < 2 {
327        return Err(FdarError::InvalidDimension {
328            parameter: "scores",
329            expected: ">=2 rows".into(),
330            actual: format!("{n}"),
331        });
332    }
333
334    let score_variance = compute_score_variance(&fit.fpca.scores, n, ncomp);
335    let (component_importance, importance, importance_normalized) =
336        compute_pointwise_importance_core(
337            &fit.coefficients,
338            &fit.fpca.rotation,
339            &score_variance,
340            ncomp,
341            m,
342        );
343
344    Ok(PointwiseImportanceResult {
345        importance,
346        importance_normalized,
347        component_importance,
348        score_variance,
349    })
350}
351
352/// Pointwise variable importance for a functional logistic regression model.
353///
354/// # Errors
355///
356/// Returns [`FdarError::InvalidParameter`] if `fit.ncomp` is zero.
357/// Returns [`FdarError::InvalidDimension`] if the rotation matrix has zero rows
358/// or the scores matrix has fewer than 2 rows.
359#[must_use = "expensive computation whose result should not be discarded"]
360pub fn pointwise_importance_logistic(
361    fit: &FunctionalLogisticResult,
362) -> Result<PointwiseImportanceResult, FdarError> {
363    let ncomp = fit.ncomp;
364    let m = fit.fpca.rotation.nrows();
365    let n = fit.fpca.scores.nrows();
366    if ncomp == 0 {
367        return Err(FdarError::InvalidParameter {
368            parameter: "ncomp",
369            message: "must be > 0".into(),
370        });
371    }
372    if m == 0 {
373        return Err(FdarError::InvalidDimension {
374            parameter: "rotation",
375            expected: ">0 rows".into(),
376            actual: "0".into(),
377        });
378    }
379    if n < 2 {
380        return Err(FdarError::InvalidDimension {
381            parameter: "scores",
382            expected: ">=2 rows".into(),
383            actual: format!("{n}"),
384        });
385    }
386
387    let score_variance = compute_score_variance(&fit.fpca.scores, n, ncomp);
388    let (component_importance, importance, importance_normalized) =
389        compute_pointwise_importance_core(
390            &fit.coefficients,
391            &fit.fpca.rotation,
392            &score_variance,
393            ncomp,
394            m,
395        );
396
397    Ok(PointwiseImportanceResult {
398        importance,
399        importance_normalized,
400        component_importance,
401        score_variance,
402    })
403}
404
405/// Compute component importance matrix and aggregated importance.
406fn compute_pointwise_importance_core(
407    coefficients: &[f64],
408    rotation: &FdMatrix,
409    score_variance: &[f64],
410    ncomp: usize,
411    m: usize,
412) -> (FdMatrix, Vec<f64>, Vec<f64>) {
413    let mut component_importance = FdMatrix::zeros(ncomp, m);
414    for k in 0..ncomp {
415        let ck = coefficients[1 + k];
416        for j in 0..m {
417            component_importance[(k, j)] = (ck * rotation[(j, k)]).powi(2) * score_variance[k];
418        }
419    }
420
421    let mut importance = vec![0.0; m];
422    for j in 0..m {
423        for k in 0..ncomp {
424            importance[j] += component_importance[(k, j)];
425        }
426    }
427
428    let total: f64 = importance.iter().sum();
429    let importance_normalized = if total > 0.0 {
430        importance.iter().map(|&v| v / total).collect()
431    } else {
432        vec![0.0; m]
433    };
434
435    (component_importance, importance, importance_normalized)
436}
437
438// ===========================================================================
439// Conditional Permutation Importance
440// ===========================================================================
441
442/// Result of conditional permutation importance.
443#[derive(Debug, Clone, PartialEq)]
444#[non_exhaustive]
445pub struct ConditionalPermutationImportanceResult {
446    /// Conditional importance per FPC component, length ncomp.
447    pub importance: Vec<f64>,
448    /// Baseline metric (R^2 or accuracy).
449    pub baseline_metric: f64,
450    /// Mean metric after conditional permutation, length ncomp.
451    pub permuted_metric: Vec<f64>,
452    /// Unconditional (standard) permutation importance for comparison, length ncomp.
453    pub unconditional_importance: Vec<f64>,
454}
455
456/// Conditional permutation importance for a linear functional regression model.
457///
458/// # Errors
459///
460/// Returns [`FdarError::InvalidDimension`] if `data` has zero rows, its column
461/// count does not match `fit.fpca.mean`, or `y.len()` does not match the row
462/// count.
463/// Returns [`FdarError::InvalidParameter`] if `n_perm` or `n_bins` is zero.
464/// Returns [`FdarError::ComputationFailed`] if the total sum of squares is zero.
465#[must_use = "expensive computation whose result should not be discarded"]
466pub fn conditional_permutation_importance(
467    fit: &FregreLmResult,
468    data: &FdMatrix,
469    y: &[f64],
470    scalar_covariates: Option<&FdMatrix>,
471    n_bins: usize,
472    n_perm: usize,
473    seed: u64,
474) -> Result<ConditionalPermutationImportanceResult, FdarError> {
475    let (n, m) = data.shape();
476    if n == 0 {
477        return Err(FdarError::InvalidDimension {
478            parameter: "data",
479            expected: ">0 rows".into(),
480            actual: "0".into(),
481        });
482    }
483    if n != y.len() {
484        return Err(FdarError::InvalidDimension {
485            parameter: "y",
486            expected: format!("{n} (matching data rows)"),
487            actual: format!("{}", y.len()),
488        });
489    }
490    if m != fit.fpca.mean.len() {
491        return Err(FdarError::InvalidDimension {
492            parameter: "data",
493            expected: format!("{} columns", fit.fpca.mean.len()),
494            actual: format!("{m}"),
495        });
496    }
497    if n_perm == 0 {
498        return Err(FdarError::InvalidParameter {
499            parameter: "n_perm",
500            message: "must be > 0".into(),
501        });
502    }
503    if n_bins == 0 {
504        return Err(FdarError::InvalidParameter {
505            parameter: "n_bins",
506            message: "must be > 0".into(),
507        });
508    }
509    let _ = scalar_covariates;
510    let ncomp = fit.ncomp;
511    let scores = project_scores(
512        data,
513        &fit.fpca.mean,
514        &fit.fpca.rotation,
515        ncomp,
516        &fit.fpca.weights,
517    );
518
519    let y_mean: f64 = y.iter().sum::<f64>() / n as f64;
520    let ss_tot: f64 = y.iter().map(|&yi| (yi - y_mean).powi(2)).sum();
521    if ss_tot == 0.0 {
522        return Err(FdarError::ComputationFailed {
523            operation: "conditional_permutation_importance",
524            detail: "total sum of squares is zero; all response values may be identical — check your data".into(),
525        });
526    }
527    let ss_res_base: f64 = fit.residuals.iter().map(|r| r * r).sum();
528    let baseline = 1.0 - ss_res_base / ss_tot;
529
530    let predict_r2 = |score_mat: &FdMatrix| -> f64 {
531        let ss_res: f64 = (0..n)
532            .map(|i| {
533                let mut yhat = fit.coefficients[0];
534                for c in 0..ncomp {
535                    yhat += fit.coefficients[1 + c] * score_mat[(i, c)];
536                }
537                (y[i] - yhat).powi(2)
538            })
539            .sum();
540        1.0 - ss_res / ss_tot
541    };
542
543    // NOT migrated to permutation_test::permutation_pvalue — uses a single ADVANCING StdRng across ALL
544    // components; per-perm reseed would change the importances (Phase-49 CONS-02 Plan A). The Phase-48
545    // hand-off "fold into the already-parallel generic path" is behavior-CHANGING (generic path reseeds
546    // per-component seed+k; this advances one RNG) → DEFERRED to a future behavior-changing phase with a
547    // re-baselined golden, NOT done here.
548    let mut rng = StdRng::seed_from_u64(seed);
549    let mut importance = vec![0.0; ncomp];
550    let mut permuted_metric = vec![0.0; ncomp];
551    let mut unconditional_importance = vec![0.0; ncomp];
552
553    for k in 0..ncomp {
554        let bins = compute_conditioning_bins(&scores, ncomp, k, n, n_bins);
555        let (mean_cond, mean_uncond) =
556            permute_component(&scores, &bins, k, n, ncomp, n_perm, &mut rng, &predict_r2);
557        permuted_metric[k] = mean_cond;
558        importance[k] = baseline - mean_cond;
559        unconditional_importance[k] = baseline - mean_uncond;
560    }
561
562    Ok(ConditionalPermutationImportanceResult {
563        importance,
564        baseline_metric: baseline,
565        permuted_metric,
566        unconditional_importance,
567    })
568}
569
570/// Conditional permutation importance for a functional logistic regression model.
571///
572/// # Errors
573///
574/// Returns [`FdarError::InvalidDimension`] if `data` has zero rows, its column
575/// count does not match `fit.fpca.mean`, or `y.len()` does not match the row
576/// count.
577/// Returns [`FdarError::InvalidParameter`] if `n_perm` or `n_bins` is zero.
578#[must_use = "expensive computation whose result should not be discarded"]
579pub fn conditional_permutation_importance_logistic(
580    fit: &FunctionalLogisticResult,
581    data: &FdMatrix,
582    y: &[f64],
583    scalar_covariates: Option<&FdMatrix>,
584    n_bins: usize,
585    n_perm: usize,
586    seed: u64,
587) -> Result<ConditionalPermutationImportanceResult, FdarError> {
588    let (n, m) = data.shape();
589    if n == 0 {
590        return Err(FdarError::InvalidDimension {
591            parameter: "data",
592            expected: ">0 rows".into(),
593            actual: "0".into(),
594        });
595    }
596    if n != y.len() {
597        return Err(FdarError::InvalidDimension {
598            parameter: "y",
599            expected: format!("{n} (matching data rows)"),
600            actual: format!("{}", y.len()),
601        });
602    }
603    if m != fit.fpca.mean.len() {
604        return Err(FdarError::InvalidDimension {
605            parameter: "data",
606            expected: format!("{} columns", fit.fpca.mean.len()),
607            actual: format!("{m}"),
608        });
609    }
610    if n_perm == 0 {
611        return Err(FdarError::InvalidParameter {
612            parameter: "n_perm",
613            message: "must be > 0".into(),
614        });
615    }
616    if n_bins == 0 {
617        return Err(FdarError::InvalidParameter {
618            parameter: "n_bins",
619            message: "must be > 0".into(),
620        });
621    }
622    let _ = scalar_covariates;
623    let ncomp = fit.ncomp;
624    let scores = project_scores(
625        data,
626        &fit.fpca.mean,
627        &fit.fpca.rotation,
628        ncomp,
629        &fit.fpca.weights,
630    );
631
632    let baseline: f64 = (0..n)
633        .filter(|&i| {
634            let pred = if fit.probabilities[i] >= 0.5 {
635                1.0
636            } else {
637                0.0
638            };
639            (pred - y[i]).abs() < 1e-10
640        })
641        .count() as f64
642        / n as f64;
643
644    let predict_acc = |score_mat: &FdMatrix| -> f64 {
645        let correct: usize = (0..n)
646            .filter(|&i| {
647                let mut eta = fit.intercept;
648                for c in 0..ncomp {
649                    eta += fit.coefficients[1 + c] * score_mat[(i, c)];
650                }
651                let pred = if sigmoid(eta) >= 0.5 { 1.0 } else { 0.0 };
652                (pred - y[i]).abs() < 1e-10
653            })
654            .count();
655        correct as f64 / n as f64
656    };
657
658    // NOT migrated to permutation_test::permutation_pvalue — uses a single ADVANCING StdRng across ALL
659    // components; per-perm reseed would change the importances (Phase-49 CONS-02 Plan A). The Phase-48
660    // hand-off "fold into the already-parallel generic path" is behavior-CHANGING (generic path reseeds
661    // per-component seed+k; this advances one RNG) → DEFERRED to a future behavior-changing phase with a
662    // re-baselined golden, NOT done here.
663    let mut rng = StdRng::seed_from_u64(seed);
664    let mut importance = vec![0.0; ncomp];
665    let mut permuted_metric = vec![0.0; ncomp];
666    let mut unconditional_importance = vec![0.0; ncomp];
667
668    for k in 0..ncomp {
669        let bins = compute_conditioning_bins(&scores, ncomp, k, n, n_bins);
670        let (mean_cond, mean_uncond) =
671            permute_component(&scores, &bins, k, n, ncomp, n_perm, &mut rng, &predict_acc);
672        permuted_metric[k] = mean_cond;
673        importance[k] = baseline - mean_cond;
674        unconditional_importance[k] = baseline - mean_uncond;
675    }
676
677    Ok(ConditionalPermutationImportanceResult {
678        importance,
679        baseline_metric: baseline,
680        permuted_metric,
681        unconditional_importance,
682    })
683}
fdars_core/explain/importance.rs

fdars_core/explain/
importance.rs
