fdars_core/explain/
pdp.rs

1//! PDP/ICE and beta decomposition.
2
3use super::helpers::*;
4use crate::matrix::FdMatrix;
5use crate::scalar_on_function::{sigmoid, FregreLmResult, FunctionalLogisticResult};
6
7/// Result of a functional partial dependence plot.
8pub struct FunctionalPdpResult {
9    /// FPC score grid values (length n_grid).
10    pub grid_values: Vec<f64>,
11    /// Average prediction across observations at each grid point (length n_grid).
12    pub pdp_curve: Vec<f64>,
13    /// Individual conditional expectation curves (n x n_grid).
14    pub ice_curves: FdMatrix,
15    /// Which FPC component was varied.
16    pub component: usize,
17}
18
19/// Functional PDP/ICE for a linear functional regression model.
20///
21/// Varies the FPC score for `component` across a grid while keeping other scores
22/// fixed, producing ICE curves and their average (PDP).
23///
24/// For a linear model, ICE curves are parallel lines (same slope, different intercepts).
25///
26/// # Arguments
27/// * `fit` -- A fitted [`FregreLmResult`]
28/// * `data` -- Original functional predictor matrix (n x m)
29/// * `scalar_covariates` -- Optional scalar covariates (n x p)
30/// * `component` -- Which FPC component to vary (0-indexed, must be < fit.ncomp)
31/// * `n_grid` -- Number of grid points (must be >= 2)
32pub fn functional_pdp(
33    fit: &FregreLmResult,
34    data: &FdMatrix,
35    _scalar_covariates: Option<&FdMatrix>,
36    component: usize,
37    n_grid: usize,
38) -> Option<FunctionalPdpResult> {
39    let (n, m) = data.shape();
40    if component >= fit.ncomp
41        || n_grid < 2
42        || n == 0
43        || m != fit.fpca.mean.len()
44        || n != fit.fitted_values.len()
45    {
46        return None;
47    }
48
49    let ncomp = fit.ncomp;
50    let scores = project_scores(data, &fit.fpca.mean, &fit.fpca.rotation, ncomp);
51    let grid_values = make_grid(&scores, component, n_grid);
52
53    let coef_c = fit.coefficients[1 + component];
54    let mut ice_curves = FdMatrix::zeros(n, n_grid);
55    for i in 0..n {
56        let base = fit.fitted_values[i] - coef_c * scores[(i, component)];
57        for g in 0..n_grid {
58            ice_curves[(i, g)] = base + coef_c * grid_values[g];
59        }
60    }
61
62    let pdp_curve = ice_to_pdp(&ice_curves, n, n_grid);
63
64    Some(FunctionalPdpResult {
65        grid_values,
66        pdp_curve,
67        ice_curves,
68        component,
69    })
70}
71
72/// Functional PDP/ICE for a functional logistic regression model.
73///
74/// Predictions pass through sigmoid, so ICE curves are non-parallel.
75///
76/// # Arguments
77/// * `fit` -- A fitted [`FunctionalLogisticResult`]
78/// * `data` -- Original functional predictor matrix (n x m)
79/// * `scalar_covariates` -- Optional scalar covariates (n x p)
80/// * `component` -- Which FPC component to vary (0-indexed, must be < fit.ncomp)
81/// * `n_grid` -- Number of grid points (must be >= 2)
82pub fn functional_pdp_logistic(
83    fit: &FunctionalLogisticResult,
84    data: &FdMatrix,
85    scalar_covariates: Option<&FdMatrix>,
86    component: usize,
87    n_grid: usize,
88) -> Option<FunctionalPdpResult> {
89    let (n, m) = data.shape();
90    if component >= fit.ncomp || n_grid < 2 || n == 0 || m != fit.fpca.mean.len() {
91        return None;
92    }
93
94    let ncomp = fit.ncomp;
95    let p_scalar = fit.gamma.len();
96    if p_scalar > 0 && scalar_covariates.is_none() {
97        return None;
98    }
99
100    let scores = project_scores(data, &fit.fpca.mean, &fit.fpca.rotation, ncomp);
101    let grid_values = make_grid(&scores, component, n_grid);
102
103    let mut ice_curves = FdMatrix::zeros(n, n_grid);
104    let coef_c = fit.coefficients[1 + component];
105    for i in 0..n {
106        let eta_base = logistic_eta_base(
107            fit.intercept,
108            &fit.coefficients,
109            &fit.gamma,
110            &scores,
111            scalar_covariates,
112            i,
113            ncomp,
114            component,
115        );
116        for g in 0..n_grid {
117            ice_curves[(i, g)] = sigmoid(eta_base + coef_c * grid_values[g]);
118        }
119    }
120
121    let pdp_curve = ice_to_pdp(&ice_curves, n, n_grid);
122
123    Some(FunctionalPdpResult {
124        grid_values,
125        pdp_curve,
126        ice_curves,
127        component,
128    })
129}
130
131// ---------------------------------------------------------------------------
132// Beta decomposition
133// ---------------------------------------------------------------------------
134
135/// Per-FPC decomposition of the functional coefficient beta(t).
136pub struct BetaDecomposition {
137    /// `components[k]` = coef_k * phi_k(t), each of length m.
138    pub components: Vec<Vec<f64>>,
139    /// FPC regression coefficients (length ncomp).
140    pub coefficients: Vec<f64>,
141    /// Proportion of ||beta(t)||^2 explained by each component.
142    pub variance_proportion: Vec<f64>,
143}
144
145/// Decompose beta(t) = sum_k coef_k * phi_k(t) for a linear functional regression.
146pub fn beta_decomposition(fit: &FregreLmResult) -> Option<BetaDecomposition> {
147    let ncomp = fit.ncomp;
148    let m = fit.fpca.mean.len();
149    if ncomp == 0 || m == 0 {
150        return None;
151    }
152    decompose_beta(&fit.coefficients, &fit.fpca.rotation, ncomp, m)
153}
154
155/// Decompose beta(t) for a functional logistic regression.
156pub fn beta_decomposition_logistic(fit: &FunctionalLogisticResult) -> Option<BetaDecomposition> {
157    let ncomp = fit.ncomp;
158    let m = fit.fpca.mean.len();
159    if ncomp == 0 || m == 0 {
160        return None;
161    }
162    decompose_beta(&fit.coefficients, &fit.fpca.rotation, ncomp, m)
163}
164
165fn decompose_beta(
166    coefficients: &[f64],
167    rotation: &FdMatrix,
168    ncomp: usize,
169    m: usize,
170) -> Option<BetaDecomposition> {
171    let mut components = Vec::with_capacity(ncomp);
172    let mut coefs = Vec::with_capacity(ncomp);
173    let mut norms_sq = Vec::with_capacity(ncomp);
174
175    for k in 0..ncomp {
176        let ck = coefficients[1 + k];
177        coefs.push(ck);
178        let comp: Vec<f64> = (0..m).map(|j| ck * rotation[(j, k)]).collect();
179        let nsq: f64 = comp.iter().map(|v| v * v).sum();
180        norms_sq.push(nsq);
181        components.push(comp);
182    }
183
184    let total_sq: f64 = norms_sq.iter().sum();
185    let variance_proportion = if total_sq > 0.0 {
186        norms_sq.iter().map(|&s| s / total_sq).collect()
187    } else {
188        vec![0.0; ncomp]
189    };
190
191    Some(BetaDecomposition {
192        components,
193        coefficients: coefs,
194        variance_proportion,
195    })
196}
197
198// ---------------------------------------------------------------------------
199// Significant regions
200// ---------------------------------------------------------------------------
201
202/// Direction of a significant region.
203#[derive(Debug, Clone, Copy, PartialEq, Eq)]
204pub enum SignificanceDirection {
205    Positive,
206    Negative,
207}
208
209/// A contiguous interval where the confidence band excludes zero.
210#[derive(Debug, Clone)]
211pub struct SignificantRegion {
212    /// Start index (inclusive).
213    pub start_idx: usize,
214    /// End index (inclusive).
215    pub end_idx: usize,
216    /// Direction of the effect.
217    pub direction: SignificanceDirection,
218}
219
220/// Identify contiguous regions where the CI `[lower, upper]` excludes zero.
221pub fn significant_regions(lower: &[f64], upper: &[f64]) -> Option<Vec<SignificantRegion>> {
222    if lower.len() != upper.len() || lower.is_empty() {
223        return None;
224    }
225    let n = lower.len();
226    let mut regions = Vec::new();
227    let mut i = 0;
228    while i < n {
229        if let Some(d) = detect_direction(lower[i], upper[i]) {
230            let start = i;
231            i += 1;
232            while i < n && detect_direction(lower[i], upper[i]) == Some(d) {
233                i += 1;
234            }
235            regions.push(SignificantRegion {
236                start_idx: start,
237                end_idx: i - 1,
238                direction: d,
239            });
240        } else {
241            i += 1;
242        }
243    }
244    Some(regions)
245}
246
247/// Build CI from beta(t) +/- z * SE, then find significant regions.
248pub fn significant_regions_from_se(
249    beta_t: &[f64],
250    beta_se: &[f64],
251    z_alpha: f64,
252) -> Option<Vec<SignificantRegion>> {
253    if beta_t.len() != beta_se.len() || beta_t.is_empty() {
254        return None;
255    }
256    let lower: Vec<f64> = beta_t
257        .iter()
258        .zip(beta_se)
259        .map(|(b, s)| b - z_alpha * s)
260        .collect();
261    let upper: Vec<f64> = beta_t
262        .iter()
263        .zip(beta_se)
264        .map(|(b, s)| b + z_alpha * s)
265        .collect();
266    significant_regions(&lower, &upper)
267}
268
269/// Detect significance direction at a single point from CI bounds.
270fn detect_direction(lower: f64, upper: f64) -> Option<SignificanceDirection> {
271    if lower > 0.0 {
272        Some(SignificanceDirection::Positive)
273    } else if upper < 0.0 {
274        Some(SignificanceDirection::Negative)
275    } else {
276        None
277    }
278}
fdars_core/explain/pdp.rs

fdars_core/explain/
pdp.rs
