R/fpca_utils.R

In FPScausal: Functional Propensity Score for Causal Inference

#' Internal FPCA helper
#'
#' Performs Functional Principal Component Analysis on a matrix of functional
#' observations using a B-spline basis expansion, then retains enough components
#' to explain at least `pve` of the total variance.  Accepts either an n x T
#' numeric matrix or an \code{fd} object from the \pkg{fda} package.
#'
#' @param X n x T numeric matrix of observed functional data (rows = subjects),
#'   or an \code{fd} object.
#' @param pve Proportion of variance explained threshold (default 0.95).
#' @param t_grid Numeric vector of observation time points (length T).
#'   Inferred from \code{X} when \code{X} is an \code{fd} object.
#' @param domain Numeric vector c(a, b) giving the domain of the functions.
#'   Inferred from \code{X} when \code{X} is an \code{fd} object.
#' @param nbasis Number of B-spline basis functions. If NULL, set automatically.
#' @return A named list with components:
#'   \item{scr}{n x L matrix of FPC scores.}
#'   \item{efn}{T x L matrix of eigenfunction values evaluated on t_grid.}
#'   \item{mean}{Numeric vector of length T: mean function evaluated on t_grid.}
#'   \item{eval}{Numeric vector of eigenvalues (all, not just L).}
#'   \item{varprop}{Numeric vector of variance proportions for retained
#'     components.}
#'   \item{perc}{Cumulative variance proportions (all components).}
#'   \item{pca_fd}{The raw \code{pca.fd} object from the \pkg{fda} package.}
#'   \item{L}{Number of retained components.}
#'   \item{t_grid}{The t_grid used.}
#'   \item{domain}{The domain used.}
#' @keywords internal
.fps_fpca <- function(X, pve = 0.95, t_grid = NULL, domain = NULL,
                       nbasis = NULL) {

  # ---- If X is already an fd object, use it directly ----
  if (inherits(X, "fd")) {
    rng    <- X$basis$rangeval
    if (is.null(domain))  domain  <- rng
    if (is.null(t_grid))  t_grid  <- seq(rng[1], rng[2], length.out = 100)
    fd_obj <- X

  } else {
    # ---- Build fd from raw matrix ----
    if (!is.matrix(X)) X <- as.matrix(X)
    if (is.null(nbasis)) nbasis <- .auto_nbasis(t_grid)
    basis  <- fda::create.bspline.basis(rangeval = domain, nbasis = nbasis)
    fd_obj <- fda::Data2fd(argvals = t_grid, y = t(X), basisobj = basis)
  }

  pca0  <- fda::pca.fd(fd_obj, centerfns = TRUE)
  pves  <- cumsum(pca0$values) / sum(pca0$values)
  ncomp <- which(pves >= pve)[1]
  if (is.na(ncomp)) {
    ncomp <- length(pca0$values)
  }

  pca <- fda::pca.fd(fd_obj, nharm = ncomp, centerfns = TRUE)

  efn <- fda::eval.fd(t_grid, pca$harmonics)
  colnames(efn) <- NULL

  list(
    scr     = pca$scores,
    efn     = efn,
    mean    = as.vector(fda::eval.fd(t_grid, pca$meanfd)),
    eval    = pca$values,
    varprop = pca$varprop,
    perc    = pves,
    pca_fd  = pca,
    L       = ncomp,
    t_grid  = t_grid,
    domain  = domain
  )
}