FPCA.FEM {fdaPDE}

Smooth Functional Principal Component Analysis

Description

This function implements a smooth functional principal component analysis over a planar mesh, a smooth manifold or a volume.

Usage

FPCA.FEM(locations = NULL, datamatrix, FEMbasis, lambda, nPC = 1, validation = NULL,
                NFolds = 5,GCVmethod = "Stochastic", nrealizations = 100, search = "tree",
                bary.locations = NULL)

Arguments

locations

A #observations-by-2 matrix in the 2D case and #observations-by-3 matrix in the 2.5D and 3D case, where each row specifies the spatial coordinates x and y (and z in 2.5D and 3D) of the corresponding observation in the datamatrix. If the locations of the observations coincide with (or are a subset of) the nodes of the mesh in the FEMbasis, leave the parameter locations = NULL for a faster implementation.

datamatrix

A matrix of dimensions #samples-by-#locations with the observed data values over the domain for each sample. The datamatrix needs to have zero mean. If the locations argument is left NULL the datamatrix has to be dimensions #samples-by-#nodes where #nodes is the number of nodes of the mesh in the FEMbasis. In this case, each observation is associated to the corresponding node in the mesh. If the data are observed only on a subset of the mesh nodes, fill with NA the values of the datamatrix in correspondence of unobserved data.

FEMbasis

A FEMbasis object describing the Finite Element basis, as created by create.FEM.basis.

lambda

A scalar or vector of smoothing parameters.

nPC

An integer specifying the number of Principal Components to compute.

validation

A string specifying the type of validation to perform. If lambda is a vector, it has to be specified as "GCV" or "KFold". This parameter specify which method of cross-validation is used to select the best parameter lambda among those values of the smoothing parameter specified in lambda for each Principal Component.

NFolds

This parameter is used only in case validation = "KFold". It is an integer specifying the number of folds to use if the KFold cross-validation method for the selection of the best parameter lambda is chosen. Default value is 5.

GCVmethod

This parameter is considered only when validation = "GCV". It can be either "Exact" or "Stochastic". If set to "Exact" the algoritm performs an exact (but possibly slow) computation of the GCV index. If set to "Stochastic" the GCV is approximated by a stochastic algorithm.

nrealizations

The number of realizations to be used in the stochastic algorithm for the estimation of GCV.

search

a flag to decide the search algorithm type (tree or naive or walking search algorithm).

bary.locations

A list with three vectors: locations, location points which are same as the given locations options. (checks whether both locations are the same); element ids, a vector of element id of the points from the mesh where they are located; barycenters, a vector of barycenter of points from the located element.

Value

A list with the following variables:

loadings.FEM

A FEM object that represents the L^2-normalized functional loadings for each Principal Component computed.

scores

A #samples-by-#PrincipalComponents matrix that represents the unnormalized scores or PC vectors.

lambda

A vector of length #PrincipalComponents with the values of the smoothing parameter lambda chosen for that Principal Component.

variance_explained

A vector of length #PrincipalComponents where each value represent the variance explained by that component.

cumsum_percentage

A vector of length #PrincipalComponents containing the cumulative percentage of the variance explained by the first components.

bary.locations

A barycenter information of the given locations if the locations are not mesh nodes.

References

Lila, E., Aston, J.A.D., Sangalli, L.M., 2016a. Smooth Principal Component Analysis over two-dimensional manifolds with an application to neuroimaging. Ann. Appl. Stat., 10(4), pp. 1854-1879.

Examples

library(fdaPDE)

## Load the hub data
data(hub2.5D)
hub2.5D.nodes = hub2.5D$hub2.5D.nodes
hub2.5D.triangles = hub2.5D$hub2.5D.triangles

mesh = create.mesh.2.5D(nodes = hub2.5D.nodes, triangles = hub2.5D.triangles)
## Create the Finite Element basis
FEMbasis = create.FEM.basis(mesh)
## Create a datamatrix
datamatrix = NULL
for(ii in 1:50){
  a1 = rnorm(1, mean = 1, sd = 1)
  a2 = rnorm(1, mean = 1, sd = 1)
  a3 = rnorm(1, mean = 1, sd = 1)

  func_evaluation = numeric(nrow(mesh$nodes))
  for (i in 0:(nrow(mesh$nodes)-1)){
    func_evaluation[i+1] = a1* sin(2*pi*mesh$nodes[i+1,1]) +
                           a2* sin(2*pi*mesh$nodes[i+1,2]) +
                           a3*sin(2*pi*mesh$nodes[i+1,3]) + 1
  }
  data = func_evaluation + rnorm(nrow(mesh$nodes), mean = 0, sd = 0.5)
  datamatrix = rbind(datamatrix, data)
}
## Compute the mean of the datamatrix and subtract it to the data
data_bar = colMeans(datamatrix)
data_demean = matrix(rep(data_bar,50), nrow=50, byrow=TRUE)

datamatrix_demeaned = datamatrix - data_demean
## Set the smoothing parameter lambda
lambda = 0.00375
## Estimate the first 2 Principal Components
FPCA_solution = FPCA.FEM(datamatrix = datamatrix_demeaned,
                      FEMbasis = FEMbasis, lambda = lambda, nPC = 2)

## Plot the functional loadings of the estimated Principal Components
plot(FPCA_solution$loadings.FEM)