FPCA.FEM: Smooth Functional Principal Component Analysis

Description Usage Arguments Value References Examples

View source: R/smoothFPCA.R

Description

This function implements a smooth functional principal component analysis for data defined over a planar mesh, or a smooth manifold. For details on the model see Lila et al. 2016.

Usage

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

Arguments

locations

A #observations-by-2 matrix in the 2D case and #observations-by-3 matrix in the 2.5D case, where each row specifies the spatial coordinates x and y (and z in 2.5D) 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.

Value

A list with the following variables:

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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
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(mesh$nnodes)
  for (i in 0:(mesh$nnodes-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(mesh$nnodes, 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)

fdaPDE documentation built on July 2, 2020, 2:22 a.m.