Skip to contents

Performs Sparse and Functional PCA on a data matrix, allowing for both sparsity and smoothness in the estimated principal components. Includes heuristic estimation of penalty parameters. The spatial smoothness penalty is constructed based on provided spatial coordinates.

Usage

second_diff_matrix(n)

sfpca(
  X,
  K,
  spat_cds,
  lambda_u = NULL,
  lambda_v = NULL,
  alpha_u = NULL,
  alpha_v = NULL,
  Omega_u = NULL,
  penalty_u = "l1",
  penalty_v = "l1",
  uthresh = 0.9,
  vthresh = 0.9,
  knn = min(6, ncol(X) - 1),
  max_iter = 100,
  tol = 1e-06,
  verbose = FALSE
)

construct_spatial_penalty(spat_cds, method = "distance", k = 6L)

Arguments

X

A numeric data matrix of dimensions n (observations/time points) by p (variables/space).

K

The number of principal components to estimate.

spat_cds

A matrix of spatial coordinates for each column of X (variables). Each row corresponds to a spatial dimension (e.g., x, y, z), and each column corresponds to a variable.

lambda_u

Sparsity penalty parameter for u. If NULL, estimated heuristically.

lambda_v

Sparsity penalty parameter for v. If NULL, estimated heuristically.

alpha_u

Smoothness penalty parameter for u. If NULL, estimated heuristically.

alpha_v

Smoothness penalty parameter for v. If NULL, estimated heuristically.

Omega_u

A positive semi-definite matrix for smoothness penalty on u. If NULL, defaults to second differences penalty (sparse matrix).

penalty_u

The penalty function for u. Can be "l1" (lasso), "scad", or a custom function.

penalty_v

The penalty function for v. Can be "l1" (lasso), "scad", or a custom function.

uthresh

Quantile for selecting `lambda_u` when estimated.

vthresh

Quantile for selecting `lambda_v` when estimated.

knn

Number of nearest neighbours for constructing `Omega_v`.

max_iter

Maximum number of iterations for the alternating optimization.

tol

Tolerance for convergence.

verbose

Logical; if TRUE, prints progress messages.

Value

A list containing the estimated singular values `d`, left singular vectors `u`, right singular vectors `v`, and the penalty parameters used.

Examples

library(Matrix)
set.seed(123)
# Simulate a small example due to resource constraints
n <- 100  # Number of time points
p <- 50   # Number of spatial locations
X <- matrix(rnorm(n * p), n, p)
# Simulate spatial coordinates
spat_cds <- matrix(runif(p * 3), nrow = 3, ncol = p)  # 3D coordinates
result <- sfpca(X, K = 2, spat_cds = spat_cds)