fdars_core/scalar_on_function/
mod.rs

1//! Scalar-on-function regression with mixed scalar/functional covariates.
2//!
3//! Implements models of the form:
4//! ```text
5//! y = α + ∫β(t)X(t)dt + γᵀz + ε
6//! ```
7//! where X(t) is a functional predictor, z is a vector of scalar covariates,
8//! β(t) is the functional coefficient, and γ is the vector of scalar coefficients.
9//!
10//! # Methods
11//!
12//! - [`fregre_lm`]: FPC-based functional linear model with optional scalar covariates
13//! - [`fregre_l1`]: L1 (median) robust functional regression via IRLS
14//! - [`fregre_huber`]: Huber M-estimation robust functional regression via IRLS
15//! - [`fregre_np_mixed`]: Nonparametric kernel regression with product kernels
16//! - [`functional_logistic`]: Logistic regression for binary outcomes
17//! - [`fregre_cv`]: Cross-validation for number of FPC components
18
19use crate::error::FdarError;
20use crate::linalg::cholesky_solve as linalg_cholesky_solve;
21use crate::matrix::FdMatrix;
22use crate::regression::{FpcaResult, PlsResult};
23
24mod bootstrap;
25mod cv;
26mod fregre_lm;
27mod logistic;
28mod multi;
29mod nonparametric;
30mod pls;
31mod robust;
32#[cfg(test)]
33mod tests;
34
35// Re-export all public items from submodules
36pub use bootstrap::{bootstrap_ci_fregre_lm, bootstrap_ci_functional_logistic};
37pub use cv::{fregre_basis_cv, fregre_np_cv};
38pub use fregre_lm::{fregre_cv, fregre_lm, model_selection_ncomp, predict_fregre_lm};
39pub use logistic::{functional_logistic, predict_functional_logistic};
40pub use multi::{fregre_lm_multi, fregre_lm_multi_cv, predict_fregre_lm_multi, MultiCvResult};
41pub use nonparametric::{
42    fregre_np_from_distances, fregre_np_mixed, predict_fregre_np, predict_fregre_np_from_distances,
43};
44pub use pls::{fregre_pls, predict_fregre_pls};
45pub use robust::{fregre_huber, fregre_l1, predict_fregre_robust};
46
47// ---------------------------------------------------------------------------
48// Result types
49// ---------------------------------------------------------------------------
50
51/// Result of functional linear regression.
52#[derive(Debug, Clone, PartialEq)]
53#[non_exhaustive]
54pub struct FregreLmResult {
55    /// Intercept α
56    pub intercept: f64,
57    /// Functional coefficient β(t), evaluated on the original grid (length m)
58    pub beta_t: Vec<f64>,
59    /// Pointwise standard errors of β(t) (length m)
60    pub beta_se: Vec<f64>,
61    /// Scalar coefficients γ (one per scalar covariate)
62    pub gamma: Vec<f64>,
63    /// Fitted values ŷ (length n)
64    pub fitted_values: Vec<f64>,
65    /// Residuals y - ŷ (length n)
66    pub residuals: Vec<f64>,
67    /// R² statistic
68    pub r_squared: f64,
69    /// Adjusted R²
70    pub r_squared_adj: f64,
71    /// Standard errors of all coefficients (intercept, FPC scores, scalar covariates)
72    pub std_errors: Vec<f64>,
73    /// Number of FPC components used
74    pub ncomp: usize,
75    /// FPCA result (for projecting new data)
76    pub fpca: FpcaResult,
77    /// Regression coefficients on (FPC scores, scalar covariates) — internal
78    pub coefficients: Vec<f64>,
79    /// Residual standard error
80    pub residual_se: f64,
81    /// GCV criterion value (if computed)
82    pub gcv: f64,
83    /// Akaike Information Criterion
84    pub aic: f64,
85    /// Bayesian Information Criterion
86    pub bic: f64,
87}
88
89/// Result of nonparametric functional regression with mixed predictors.
90#[derive(Debug, Clone, PartialEq)]
91#[non_exhaustive]
92pub struct FregreNpResult {
93    /// Fitted values ŷ (length n)
94    pub fitted_values: Vec<f64>,
95    /// Residuals y - ŷ (length n)
96    pub residuals: Vec<f64>,
97    /// R² statistic
98    pub r_squared: f64,
99    /// Bandwidth for functional distance kernel
100    pub h_func: f64,
101    /// Bandwidth for scalar covariates kernel
102    pub h_scalar: f64,
103    /// Leave-one-out CV error
104    pub cv_error: f64,
105}
106
107/// Result of robust (L1 or Huber) functional regression.
108#[derive(Debug, Clone, PartialEq)]
109#[non_exhaustive]
110pub struct FregreRobustResult {
111    /// Intercept
112    pub intercept: f64,
113    /// Functional coefficient β(t), evaluated on the original grid (length m)
114    pub beta_t: Vec<f64>,
115    /// Fitted values ŷ (length n)
116    pub fitted_values: Vec<f64>,
117    /// Residuals y - ŷ (length n)
118    pub residuals: Vec<f64>,
119    /// Regression coefficients (intercept, FPC scores, scalar covariates)
120    pub coefficients: Vec<f64>,
121    /// Number of FPC components used
122    pub ncomp: usize,
123    /// FPCA result (for projecting new data)
124    pub fpca: FpcaResult,
125    /// Number of IRLS iterations performed
126    pub iterations: usize,
127    /// Whether the IRLS algorithm converged
128    pub converged: bool,
129    /// Final IRLS weights (length n)
130    pub weights: Vec<f64>,
131    /// R² statistic
132    pub r_squared: f64,
133}
134
135/// Result of functional logistic regression.
136#[derive(Debug, Clone, PartialEq)]
137#[non_exhaustive]
138pub struct FunctionalLogisticResult {
139    /// Intercept α
140    pub intercept: f64,
141    /// Functional coefficient β(t), evaluated on the original grid (length m)
142    pub beta_t: Vec<f64>,
143    /// Pointwise standard errors of β(t) (length m)
144    pub beta_se: Vec<f64>,
145    /// Scalar coefficients γ (one per scalar covariate)
146    pub gamma: Vec<f64>,
147    /// Predicted probabilities P(Y=1) (length n)
148    pub probabilities: Vec<f64>,
149    /// Predicted class labels (0 or 1)
150    pub predicted_classes: Vec<usize>,
151    /// Number of FPC components used
152    pub ncomp: usize,
153    /// Classification accuracy on training data
154    pub accuracy: f64,
155    /// Standard errors of all coefficients (intercept, FPC scores, scalar covariates)
156    pub std_errors: Vec<f64>,
157    /// Regression coefficients on (FPC scores, scalar covariates) — internal
158    pub coefficients: Vec<f64>,
159    /// Log-likelihood at convergence
160    pub log_likelihood: f64,
161    /// Number of IRLS iterations
162    pub iterations: usize,
163    /// FPCA result (for projecting new data)
164    pub fpca: FpcaResult,
165    /// Akaike Information Criterion
166    pub aic: f64,
167    /// Bayesian Information Criterion
168    pub bic: f64,
169}
170
171/// Result of cross-validation for K selection.
172#[derive(Debug, Clone, PartialEq)]
173#[non_exhaustive]
174pub struct FregreCvResult {
175    /// Candidate K values tested
176    pub k_values: Vec<usize>,
177    /// CV error for each K
178    pub cv_errors: Vec<f64>,
179    /// Optimal K (minimizing CV error)
180    pub optimal_k: usize,
181    /// Minimum CV error
182    pub min_cv_error: f64,
183    /// Out-of-fold predictions at optimal K (length n, each predicted when held out)
184    pub oof_predictions: Vec<f64>,
185    /// Fold assignment for each observation (0..n_folds)
186    pub fold_assignments: Vec<usize>,
187    /// Per-fold MSE at optimal K
188    pub fold_errors: Vec<f64>,
189}
190
191/// Result of PLS-based scalar-on-function regression.
192#[derive(Debug, Clone, PartialEq)]
193#[non_exhaustive]
194pub struct PlsRegressionResult {
195    /// Intercept α
196    pub intercept: f64,
197    /// Functional coefficient β(t), evaluated on the original grid (length m)
198    pub beta_t: Vec<f64>,
199    /// Scalar coefficients γ (one per scalar covariate)
200    pub gamma: Vec<f64>,
201    /// Fitted values ŷ (length n)
202    pub fitted_values: Vec<f64>,
203    /// Residuals y - ŷ (length n)
204    pub residuals: Vec<f64>,
205    /// R² statistic
206    pub r_squared: f64,
207    /// Adjusted R²
208    pub r_squared_adj: f64,
209    /// Number of PLS components used
210    pub ncomp: usize,
211    /// PLS result (for projecting new data)
212    pub pls: PlsResult,
213    /// Regression coefficients on (intercept, PLS scores, scalar covariates)
214    pub coefficients: Vec<f64>,
215    /// Residual standard error
216    pub residual_se: f64,
217    /// Akaike Information Criterion
218    pub aic: f64,
219    /// Bayesian Information Criterion
220    pub bic: f64,
221}
222
223/// Result of multi-predictor functional linear regression.
224#[derive(Debug, Clone, PartialEq)]
225#[non_exhaustive]
226pub struct MultiFregreLmResult {
227    /// Intercept α
228    pub intercept: f64,
229    /// Functional coefficients beta_k(t) for each predictor, each length m_k.
230    pub beta_t: Vec<Vec<f64>>,
231    /// Scalar coefficients γ (one per scalar covariate)
232    pub gamma: Vec<f64>,
233    /// Fitted values ŷ (length n)
234    pub fitted_values: Vec<f64>,
235    /// Residuals y - ŷ (length n)
236    pub residuals: Vec<f64>,
237    /// R²
238    pub r_squared: f64,
239    /// Adjusted R²
240    pub r_squared_adj: f64,
241    /// Number of FPC components used per functional predictor
242    pub ncomp: Vec<usize>,
243    /// FPCA results for each functional predictor (for projection)
244    pub fpcas: Vec<FpcaResult>,
245    /// Regression coefficients [intercept, scores_1..., scores_2..., ..., scalars...]
246    pub coefficients: Vec<f64>,
247    /// Residual standard error
248    pub residual_se: f64,
249    /// AIC
250    pub aic: f64,
251    /// BIC
252    pub bic: f64,
253}
254
255/// Criterion used for model selection.
256#[derive(Debug, Clone, Copy, PartialEq)]
257pub enum SelectionCriterion {
258    /// Akaike Information Criterion
259    Aic,
260    /// Bayesian Information Criterion
261    Bic,
262    /// Generalized Cross-Validation
263    Gcv,
264}
265
266/// Result of ncomp model selection.
267#[derive(Debug, Clone, PartialEq)]
268#[non_exhaustive]
269pub struct ModelSelectionResult {
270    /// Best number of FPC components by the chosen criterion
271    pub best_ncomp: usize,
272    /// (ncomp, AIC, BIC, GCV) for each candidate
273    pub criteria: Vec<(usize, f64, f64, f64)>,
274}
275
276/// Result of bootstrap confidence intervals for β(t).
277#[derive(Debug, Clone, PartialEq)]
278#[non_exhaustive]
279pub struct BootstrapCiResult {
280    /// Pointwise lower bound (length m).
281    pub lower: Vec<f64>,
282    /// Pointwise upper bound (length m).
283    pub upper: Vec<f64>,
284    /// Original β(t) estimate (length m).
285    pub center: Vec<f64>,
286    /// Simultaneous lower bound (sup-norm adjusted, length m).
287    pub sim_lower: Vec<f64>,
288    /// Simultaneous upper bound (sup-norm adjusted, length m).
289    pub sim_upper: Vec<f64>,
290    /// Number of bootstrap replicates that converged.
291    pub n_boot_success: usize,
292}
293
294/// Result of lambda selection for basis regression via cross-validation.
295#[derive(Debug, Clone, PartialEq)]
296#[non_exhaustive]
297pub struct FregreBasisCvResult {
298    /// Optimal smoothing parameter lambda.
299    pub optimal_lambda: f64,
300    /// Mean CV error for each lambda.
301    pub cv_errors: Vec<f64>,
302    /// SE of CV error across folds for each lambda.
303    pub cv_se: Vec<f64>,
304    /// Lambda values tested.
305    pub lambda_values: Vec<f64>,
306    /// Minimum mean CV error.
307    pub min_cv_error: f64,
308}
309
310/// Result of bandwidth selection for nonparametric regression via CV.
311#[derive(Debug, Clone, PartialEq)]
312#[non_exhaustive]
313pub struct FregreNpCvResult {
314    /// Optimal bandwidth.
315    pub optimal_h: f64,
316    /// Mean CV error for each bandwidth.
317    pub cv_errors: Vec<f64>,
318    /// SE of CV error across folds for each bandwidth.
319    pub cv_se: Vec<f64>,
320    /// Bandwidth values tested.
321    pub h_values: Vec<f64>,
322    /// Minimum mean CV error.
323    pub min_cv_error: f64,
324}
325
326// ---------------------------------------------------------------------------
327// Shared linear algebra helpers (delegated to crate::linalg)
328// ---------------------------------------------------------------------------
329
330// Re-export for use by submodules and explain/ modules that import from
331// `crate::scalar_on_function::{cholesky_factor, cholesky_forward_back, compute_xtx}`.
332pub(crate) use crate::linalg::cholesky_factor;
333pub(crate) use crate::linalg::cholesky_forward_back;
334pub(crate) use crate::linalg::compute_xtx;
335
336/// Compute X'y (length p).
337fn compute_xty(x: &FdMatrix, y: &[f64]) -> Vec<f64> {
338    let (n, p) = x.shape();
339    (0..p)
340        .map(|k| {
341            let mut s = 0.0;
342            for i in 0..n {
343                s += x[(i, k)] * y[i];
344            }
345            s
346        })
347        .collect()
348}
349
350/// Solve Ax = b via Cholesky decomposition (A must be symmetric positive definite).
351pub(super) fn cholesky_solve(a: &[f64], b: &[f64], p: usize) -> Result<Vec<f64>, FdarError> {
352    linalg_cholesky_solve(a, b, p)
353}
354
355/// Compute hat matrix diagonal: H_ii = x_i' (X'X)^{-1} x_i, given Cholesky factor L of X'X.
356pub(crate) fn compute_hat_diagonal(x: &FdMatrix, l: &[f64]) -> Vec<f64> {
357    let (n, p) = x.shape();
358    let mut hat_diag = vec![0.0; n];
359    for i in 0..n {
360        let mut v = vec![0.0; p];
361        for j in 0..p {
362            v[j] = x[(i, j)];
363            for k in 0..j {
364                v[j] -= l[j * p + k] * v[k];
365            }
366            v[j] /= l[j * p + j];
367        }
368        hat_diag[i] = v.iter().map(|vi| vi * vi).sum();
369    }
370    hat_diag
371}
372
373/// Compute diagonal of (X'X)^{-1} given Cholesky factor L, then SE = sqrt(sigma² * diag).
374fn compute_ols_std_errors(l: &[f64], p: usize, sigma2: f64) -> Vec<f64> {
375    let mut se = vec![0.0; p];
376    for j in 0..p {
377        let mut v = vec![0.0; p];
378        v[j] = 1.0;
379        for k in 0..p {
380            for kk in 0..k {
381                v[k] -= l[k * p + kk] * v[kk];
382            }
383            v[k] /= l[k * p + k];
384        }
385        se[j] = (sigma2 * v.iter().map(|vi| vi * vi).sum::<f64>()).sqrt();
386    }
387    se
388}
389
390// ---------------------------------------------------------------------------
391// Design matrix and coefficient recovery
392// ---------------------------------------------------------------------------
393
394/// Build design matrix: \[1, ξ_1, ..., ξ_K, z_1, ..., z_p\].
395/// Validate inputs for fregre_lm / functional_logistic.
396fn validate_fregre_inputs(
397    n: usize,
398    m: usize,
399    y: &[f64],
400    scalar_covariates: Option<&FdMatrix>,
401) -> Result<(), FdarError> {
402    if n < 3 {
403        return Err(FdarError::InvalidDimension {
404            parameter: "data",
405            expected: "at least 3 rows".to_string(),
406            actual: format!("{n}"),
407        });
408    }
409    if m == 0 {
410        return Err(FdarError::InvalidDimension {
411            parameter: "data",
412            expected: "at least 1 column".to_string(),
413            actual: "0".to_string(),
414        });
415    }
416    if y.len() != n {
417        return Err(FdarError::InvalidDimension {
418            parameter: "y",
419            expected: format!("{n}"),
420            actual: format!("{}", y.len()),
421        });
422    }
423    if let Some(sc) = scalar_covariates {
424        if sc.nrows() != n {
425            return Err(FdarError::InvalidDimension {
426                parameter: "scalar_covariates",
427                expected: format!("{n} rows"),
428                actual: format!("{} rows", sc.nrows()),
429            });
430        }
431    }
432    Ok(())
433}
434
435/// Resolve ncomp: auto-select via CV if 0, otherwise clamp.
436fn resolve_ncomp(
437    ncomp: usize,
438    data: &FdMatrix,
439    y: &[f64],
440    scalar_covariates: Option<&FdMatrix>,
441    n: usize,
442    m: usize,
443) -> Result<usize, FdarError> {
444    if ncomp == 0 {
445        let cv = fregre_cv(data, y, scalar_covariates, 1, m.min(n - 1).min(20), 5)?;
446        Ok(cv.optimal_k)
447    } else {
448        Ok(ncomp.min(n - 1).min(m))
449    }
450}
451
452pub(crate) fn build_design_matrix(
453    fpca_scores: &FdMatrix,
454    ncomp: usize,
455    scalar_covariates: Option<&FdMatrix>,
456    n: usize,
457) -> FdMatrix {
458    let p_scalar = scalar_covariates.map_or(0, super::matrix::FdMatrix::ncols);
459    let p_total = 1 + ncomp + p_scalar;
460    let mut design = FdMatrix::zeros(n, p_total);
461    for i in 0..n {
462        design[(i, 0)] = 1.0;
463        for k in 0..ncomp {
464            design[(i, 1 + k)] = fpca_scores[(i, k)];
465        }
466        if let Some(sc) = scalar_covariates {
467            for j in 0..p_scalar {
468                design[(i, 1 + ncomp + j)] = sc[(i, j)];
469            }
470        }
471    }
472    design
473}
474
475/// Recover functional coefficient β(t) = Σ_k γ_k φ_k(t).
476fn recover_beta_t(fpc_coeffs: &[f64], rotation: &FdMatrix, m: usize) -> Vec<f64> {
477    let ncomp = fpc_coeffs.len();
478    let mut beta_t = vec![0.0; m];
479    for k in 0..ncomp {
480        for j in 0..m {
481            beta_t[j] += fpc_coeffs[k] * rotation[(j, k)];
482        }
483    }
484    beta_t
485}
486
487/// Pointwise standard error of β(t) via error propagation through FPCA rotation.
488///
489/// SE[β(t_j)]² = Σ_k φ_k(t_j)² · SE[γ_k]²
490fn compute_beta_se(gamma_se: &[f64], rotation: &FdMatrix, m: usize) -> Vec<f64> {
491    let ncomp = gamma_se.len();
492    let mut beta_se = vec![0.0; m];
493    for j in 0..m {
494        let mut var_j = 0.0;
495        for k in 0..ncomp {
496            var_j += rotation[(j, k)].powi(2) * gamma_se[k].powi(2);
497        }
498        beta_se[j] = var_j.sqrt();
499    }
500    beta_se
501}
502
503/// Compute fitted values ŷ = X β.
504fn compute_fitted(design: &FdMatrix, coeffs: &[f64]) -> Vec<f64> {
505    let (n, p) = design.shape();
506    (0..n)
507        .map(|i| {
508            let mut yhat = 0.0;
509            for j in 0..p {
510                yhat += design[(i, j)] * coeffs[j];
511            }
512            yhat
513        })
514        .collect()
515}
516
517/// Compute R² and adjusted R².
518fn compute_r_squared(y: &[f64], residuals: &[f64], p_total: usize) -> (f64, f64) {
519    let n = y.len();
520    let y_mean = y.iter().sum::<f64>() / n as f64;
521    let ss_tot: f64 = y.iter().map(|&yi| (yi - y_mean).powi(2)).sum();
522    let ss_res: f64 = residuals.iter().map(|r| r * r).sum();
523    let r_squared = if ss_tot > 0.0 {
524        1.0 - ss_res / ss_tot
525    } else {
526        0.0
527    };
528    let df_model = (p_total - 1) as f64;
529    let r_squared_adj = if n as f64 - df_model - 1.0 > 0.0 {
530        1.0 - (1.0 - r_squared) * (n as f64 - 1.0) / (n as f64 - df_model - 1.0)
531    } else {
532        r_squared
533    };
534    (r_squared, r_squared_adj)
535}
536
537// ---------------------------------------------------------------------------
538// OLS solver
539// ---------------------------------------------------------------------------
540
541/// Solve ordinary least squares: min ||Xb - y||² via normal equations with Cholesky.
542/// Returns (coefficients, hat_matrix_diagonal) or error if singular.
543fn ols_solve(x: &FdMatrix, y: &[f64]) -> Result<(Vec<f64>, Vec<f64>), FdarError> {
544    let (n, p) = x.shape();
545    if n < p || p == 0 {
546        return Err(FdarError::InvalidDimension {
547            parameter: "design matrix",
548            expected: format!("n >= p and p > 0 (p={p})"),
549            actual: format!("n={n}, p={p}"),
550        });
551    }
552    let xtx = compute_xtx(x);
553    let xty = compute_xty(x, y);
554    let l = cholesky_factor(&xtx, p)?;
555    let b = cholesky_forward_back(&l, &xty, p);
556    let hat_diag = compute_hat_diagonal(x, &l);
557    Ok((b, hat_diag))
558}
559
560/// Sigmoid function: 1 / (1 + exp(-x))
561pub(crate) fn sigmoid(x: f64) -> f64 {
562    if x >= 0.0 {
563        1.0 / (1.0 + (-x).exp())
564    } else {
565        let ex = x.exp();
566        ex / (1.0 + ex)
567    }
568}
569
570// ---------------------------------------------------------------------------
571// Predict methods on result structs
572// ---------------------------------------------------------------------------
573
574impl FregreLmResult {
575    /// Predict new responses. Delegates to [`predict_fregre_lm`].
576    pub fn predict(&self, new_data: &FdMatrix, new_scalar: Option<&FdMatrix>) -> Vec<f64> {
577        predict_fregre_lm(self, new_data, new_scalar)
578    }
579}
580
581impl FregreRobustResult {
582    /// Predict new responses. Delegates to [`predict_fregre_robust`].
583    pub fn predict(&self, new_data: &FdMatrix, new_scalar: Option<&FdMatrix>) -> Vec<f64> {
584        predict_fregre_robust(self, new_data, new_scalar)
585    }
586}
587
588impl FunctionalLogisticResult {
589    /// Predict P(Y=1) for new data. Delegates to [`predict_functional_logistic`].
590    pub fn predict(&self, new_data: &FdMatrix, new_scalar: Option<&FdMatrix>) -> Vec<f64> {
591        predict_functional_logistic(self, new_data, new_scalar)
592    }
593}
594
595impl MultiFregreLmResult {
596    /// Predict new responses. Delegates to [`predict_fregre_lm_multi`].
597    ///
598    /// # Errors
599    ///
600    /// Returns [`FdarError`] if prediction fails due to dimension mismatches.
601    pub fn predict(
602        &self,
603        new_predictors: &[&FdMatrix],
604        new_scalar: Option<&FdMatrix>,
605    ) -> Result<Vec<f64>, FdarError> {
606        predict_fregre_lm_multi(self, new_predictors, new_scalar)
607    }
608}
fdars_core/scalar_on_function/mod.rs

fdars_core/scalar_on_function/
mod.rs
