Type:	Package
Title:	Multivariate Functional Principal Component Analysis for Data Observed on Different Dimensional Domains
Version:	1.3-10
Date:	2022-09-15
Maintainer:	Clara Happ-Kurz <chk_R@gmx.de>
Description:	Calculate a multivariate functional principal component analysis for data observed on different dimensional domains. The estimation algorithm relies on univariate basis expansions for each element of the multivariate functional data (Happ & Greven, 2018) <doi:10.1080/01621459.2016.1273115>. Multivariate and univariate functional data objects are represented by S4 classes for this type of data implemented in the package 'funData'. For more details on the general concepts of both packages and a case study, see Happ-Kurz (2020) <doi:10.18637/jss.v093.i05>.
URL:	https://github.com/ClaraHapp/MFPCA
License:	GPL-2
Imports:	abind, foreach, irlba, Matrix(≥ 1.5-0), methods, mgcv (≥ 1.8-33), plyr, stats
Depends:	R (≥ 3.2.0), funData (≥ 1.3-4)
Suggests:	covr, fda, testthat (≥ 2.0.0)
NeedsCompilation:	yes
SystemRequirements:	libfftw3 (>= 3.3.4)
RoxygenNote:	7.2.1
Packaged:	2022-09-15 07:48:19 UTC; clara
Author:	Clara Happ-Kurz [aut, cre]
Repository:	CRAN
Date/Publication:	2022-09-15 08:30:02 UTC

Calculate univariate functional PCA

Description

This function is a slightly adapted version of the fpca.sc function in the refund package for calculating univariate functional principal components based on a smoothed covariance function. The smoothing basis functions are penalized splines.

Usage

.PACE(
  X,
  Y,
  Y.pred = NULL,
  nbasis = 10,
  pve = 0.99,
  npc = NULL,
  makePD = FALSE,
  cov.weight.type = "none"
)

Arguments

Value

References

Di, C., Crainiceanu, C., Caffo, B., and Punjabi, N. (2009). Multilevel functional principal component analysis. Annals of Applied Statistics, 3, 458–488. Yao, F., Mueller, H.-G., and Wang, J.-L. (2005). Functional data analysis for sparse longitudinal data. Journal of the American Statistical Association, 100, 577–590.

See Also

PACE

The functional CP-TPA algorithm

Description

This function implements the functional CP-TPA (FCP-TPA) algorithm, that calculates a smooth PCA for 3D tensor data (i.e. N observations of 2D images with dimension S1 x S2). The results are given in a CANDECOMP/PARAFRAC (CP) model format

X = \sum_{k = 1}^K d_k \cdot u_k \circ v_k \circ w_k

where \circ stands for the outer product, d_k is a scalar and u_k, v_k, w_k are eigenvectors for each direction of the tensor. In this representation, the outer product v_k \circ w_k can be regarded as the k-th eigenimage, while d_k \cdot u_k represents the vector of individual scores for this eigenimage and each observation.

Usage

FCP_TPA(
  X,
  K,
  penMat,
  alphaRange,
  verbose = FALSE,
  tol = 1e-04,
  maxIter = 15,
  adaptTol = TRUE
)

Arguments

Details

The smoothness of the eigenvectors v_k, w_k is induced by penalty matrices for both image directions, that are weighted by smoothing parameters \alpha_{vk}, \alpha_{wk}. The eigenvectors u_k are not smoothed, hence the algorithm does not induce smoothness along observations.

Optimal smoothing parameters are found via a nested generalized cross validation. In each iteration of the TPA (tensor power algorithm), the GCV criterion is optimized via optimize on the interval specified via alphaRange$v (or alphaRange$w, respectively).

The FCP_TPA algorithm is an iterative algorithm. Convergence is assumed if the relative difference between the actual and the previous values are all below the tolerance level tol. The tolerance level is increased automatically, if the algorithm has not converged after maxIter steps and if adaptTol = TRUE. If the algorithm did not converge after maxIter steps (or 2 * maxIter) steps, the function throws a warning.

Value

References

G. I. Allen, "Multi-way Functional Principal Components Analysis", IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing, 2013.

See Also

fcptpaBasis

Examples

 # set.seed(1234)
 
 N <- 100
 S1 <- 75
 S2 <- 75

 # define "true" components
 v <- sin(seq(-pi, pi, length.out = S1))
 w <- exp(seq(-0.5, 1, length.out = S2))
 
 # simulate tensor data with dimensions N x S1 x S2
 X <- rnorm(N, sd = 0.5) %o% v %o% w
 
 # create penalty matrices (penalize first differences for each dimension)
 Pv <- crossprod(diff(diag(S1)))
 Pw <- crossprod(diff(diag(S2)))
 
 # estimate one eigentensor
 res <- FCP_TPA(X, K = 1, penMat = list(v = Pv, w = Pw),
             alphaRange = list(v = c(1e-4, 1e4), w = c(1e-4, 1e4)),
             verbose = TRUE)
 
 # plot the results and compare to true values
 plot(res$V)
 points(v/sqrt(sum(v^2)), pch = 20)
 legend("topleft", legend = c("True", "Estimated"), pch = c(20, 1))
 
 plot(res$W)
 points(w/sqrt(sum(w^2)), pch = 20)
 legend("topleft", legend = c("True", "Estimated"), pch = c(20, 1))

Multivariate functional principal component analysis for functions on different (dimensional) domains

Description

This function calculates a multivariate functional principal component analysis (MFPCA) based on i.i.d. observations x_1, \ldots, x_N of a multivariate functional data-generating process X = (X^{(1)}, \ldots X^{(p)}) with elements X^{(j)} \in L^2(\mathcal{T}_j) defined on a domain \mathcal{T}_j \subset IR^{d_j}. In particular, the elements can be defined on different (dimensional) domains. The results contain the mean function, the estimated multivariate functional principal components \hat \psi_1, \ldots, \hat \psi_M (having the same structure as x_i), the associated eigenvalues \hat \nu_1 \geq \ldots \geq \hat \nu_M > 0 and the individual scores \hat \rho_{im} = \widehat{<x_i, \psi_m>}. Moreover, estimated trajectories for each observation based on the truncated Karhunen-Loeve representation

\hat x_i = \sum_{m = 1}^M \hat \rho_{im} \hat \psi_m

are given if desired (fit = TRUE). The implementation of the observations x_i = (x_i^{(1)}, \ldots , x_i^{(p)}),~ i = 1 , \ldots, N, the mean function and multivariate functional principal components \hat \psi_1, \ldots, \hat \psi_M uses the multiFunData class, which is defined in the package funData.

Usage

MFPCA(
  mFData,
  M,
  uniExpansions,
  weights = rep(1, length(mFData)),
  fit = FALSE,
  approx.eigen = FALSE,
  bootstrap = FALSE,
  nBootstrap = NULL,
  bootstrapAlpha = 0.05,
  bootstrapStrat = NULL,
  verbose = options()$verbose
)

Arguments

Details

Weighted MFPCA

If the elements vary considerably in domain, range or variation, a weight vector w_1 , \ldots, w_p can be supplied and the MFPCA is based on the weighted scalar product

<<f,g>>_w = \sum_{j = 1}^p w_j \int_{\mathcal{T}_j} f^{(j)}(t) g^{(j)}(t) \mathrm{d} t

and the corresponding weighted covariance operator \Gamma_w.

Bootstrap

If bootstrap = TRUE, pointwise bootstrap confidence bands are generated for the multivariate eigenvalues \hat \nu_1, \ldots, \hat \nu_M as well as for multivariate functional principal components \hat \psi_1, \ldots, \hat \psi_M. The parameter nBootstrap gives the number of bootstrap iterations. In each iteration, the observations are resampled on the level of (multivariate) functions and the whole MFPCA is recalculated. In particular, if the univariate basis depends on the data (FPCA approaches), basis functions and scores are both re-estimated. If the basis functions are fixed (e.g. splines), the scores from the original estimate are used to speed up the calculations. The confidence bands for the eigenfunctions are calculated separately for each element as pointwise percentile bootstrap confidence intervals. Analogously, the confidence bands for the eigenvalues are also percentile bootstrap confidence bands. The significance level(s) can be defined by the bootstrapAlpha parameter, which defaults to 5%. As a result, the MFPCA function returns a list CI of the same length as bootstrapAlpha, containing the lower and upper bounds of the confidence bands for the principal components as multiFunData objects of the same structure as mFData. The confidence bands for the eigenvalues are returned in a list CIvalues, containing the upper and lower bounds for each significance level.

Univariate Expansions

The multivariate functional principal component analysis relies on a univariate basis expansion for each element X^{(j)}. The univariate basis representation is calculated using the univDecomp function, that passes the univariate functional observations and optional parameters to the specific function. The univariate decompositions are specified via the uniExpansions argument in the MFPCA function. It is a list of the same length as the mFData object, i.e. having one entry for each element of the multivariate functional data. For each element, uniExpansion must specify at least the type of basis functions to use. Additionally, one may add further parameters. The following basis representations are supported:

• Given basis functions. Then uniExpansions[[j]] = list(type = "given", functions, scores, ortho), where functions is a funData object on the same domain as mFData, containing the given basis functions. The parameters scores and ortho are optional. scores is an N x K matrix containing the scores (or coefficients) of the observed functions for the given basis functions, where N is the number of observed functions and K is the number of basis functions. Note that the scores need to be demeaned to give meaningful results. If scores are not supplied, they are calculated using the given basis functions. The parameter ortho specifies whether the given basis functions are orthonormal orhto = TRUE or not ortho = FALSE. If ortho is not supplied, the functions are treated as non-orthogonal. scores and ortho are not checked for plausibility, use them at your own risk!

• Univariate functional principal component analysis. Then uniExpansions[[j]] = list(type = "uFPCA", nbasis, pve, npc, makePD), where nbasis,pve,npc,makePD are parameters passed to the PACE function for calculating the univariate functional principal component analysis.

• Basis functions expansions from the package fda. Then uniExpansions[[j]] = list(type = "fda", ...), where ... are passed to funData2fd, which heavily builds on eval.fd. If fda is not available, a warning is thrown.

• Spline basis functions (not penalized). Then uniExpansions[[j]] = list(type = "splines1D", bs, m, k), where bs,m,k are passed to the functions univDecomp and univExpansion. For two-dimensional tensor product splines, use type = "splines2D".

• Spline basis functions (with smoothness penalty). Then uniExpansions[[j]] = list(type = "splines1Dpen", bs, m, k), where bs,m,k are passed to the functions univDecomp and univExpansion. Analogously to the unpenalized case, use type = "splines2Dpen" for 2D penalized tensor product splines.

• Cosine basis functions. Use uniExpansions[[j]] = list(type = "DCT2D", qThresh, parallel) for functions one two-dimensional domains (images) and type = "DCT3D" for 3D images. The calculation is based on the discrete cosine transform (DCT) implemented in the C-library fftw3. If this library is not available, the function will throw a warning. qThresh gives the quantile for hard thresholding the basis coefficients based on their absolute value. If parallel = TRUE, the coefficients for different images are calculated in parallel.

See univDecomp and univExpansion for details.

Value

An object of class MFPCAfit containing the following components:

References

C. Happ, S. Greven (2018): Multivariate Functional Principal Component Analysis for Data Observed on Different (Dimensional) Domains. Journal of the American Statistical Association, 113(522): 649-659. DOI: doi:10.1080/01621459.2016.1273115

C. Happ-Kurz (2020): Object-Oriented Software for Functional Data. Journal of Statistical Software, 93(5): 1-38. DOI: doi:10.18637/jss.v093.i05

See Also

See Happ-Kurz (2020. doi:10.18637/jss.v093.i05) for a general introduction to the funData package and it's interplay with MFPCA. This file also includes a case study on how to use MFPCA. Useful functions: multiFunData, PACE, univDecomp, univExpansion, summary, plot, scoreplot

Examples

oldPar <- par(no.readonly = TRUE)

set.seed(1)

### simulate data (one-dimensional domains)
sim <-  simMultiFunData(type = "split", argvals = list(seq(0,1,0.01), seq(-0.5,0.5,0.02)),
                        M = 5, eFunType = "Poly", eValType = "linear", N = 100)

# MFPCA based on univariate FPCA
uFPCA <- MFPCA(sim$simData, M = 5, uniExpansions = list(list(type = "uFPCA"),
                                                                  list(type = "uFPCA")))
summary(uFPCA)
plot(uFPCA) # plot the eigenfunctions as perturbations of the mean
scoreplot(uFPCA) # plot the scores

# MFPCA based on univariate spline expansions
splines <- MFPCA(sim$simData, M = 5, uniExpansions = list(list(type = "splines1D", k = 10),
                                                          list(type = "splines1D", k = 10)),
                 fit = TRUE) # calculate reconstruction, too
summary(splines)
plot(splines) # plot the eigenfunctions as perturbations of the mean
scoreplot(splines) # plot the scores

### Compare estimates to true eigenfunctions
# flip to make results more clear
uFPCA$functions <- flipFuns(sim$trueFuns, uFPCA$functions)
splines$functions <- flipFuns(sim$trueFuns, splines$functions)

par(mfrow = c(1,2))
plot(sim$trueFuns[[1]], main = "Eigenfunctions\n1st Element", lwd = 2)
plot(uFPCA$functions[[1]], lty = 2, add = TRUE)
plot(splines$functions[[1]], lty = 3, add = TRUE)

plot(sim$trueFuns[[2]], main = "Eigenfunctions\n2nd Element", lwd = 2)
plot(uFPCA$functions[[2]], lty = 2, add = TRUE)
plot(splines$functions[[2]], lty = 3, add = TRUE)
legend("bottomleft", c("True", "uFPCA", "splines"), lty = 1:3, lwd = c(2,1,1))

# Test reconstruction for the first 10 observations
plot(sim$simData[[1]], obs = 1:10, main = "Reconstruction\n1st Element", lwd = 2)
plot(splines$fit[[1]], obs = 1:10, lty = 2, col = 1, add = TRUE)

plot(sim$simData[[2]], obs = 1:10, main = "Reconstruction\n2nd Element", lwd = 2)
plot(splines$fit[[2]], obs = 1:10, lty = 2, col = 1, add = TRUE)
legend("bottomleft", c("True", "Reconstruction"), lty = c(1,2), lwd = c(2,1))

# MFPCA with Bootstrap-CI for the first 2 eigenfunctions
### ATTENTION: Takes long

splinesBoot <- MFPCA(sim$simData, M = 2, uniExpansions = list(list(type = "splines1D", k = 10),
                                                          list(type = "splines1D", k = 10)),
                 bootstrap = TRUE, nBootstrap = 100, bootstrapAlpha = c(0.05, 0.1), verbose = TRUE)
summary(splinesBoot)
                                 
plot(splinesBoot$functions[[1]], ylim = c(-2,1.5))
plot(splinesBoot$CI$alpha_0.05$lower[[1]], lty = 2, add = TRUE)
plot(splinesBoot$CI$alpha_0.05$upper[[1]], lty = 2, add = TRUE)
plot(splinesBoot$CI$alpha_0.1$lower[[1]], lty = 3, add = TRUE)
plot(splinesBoot$CI$alpha_0.1$upper[[1]], lty = 3, add = TRUE)
abline(h = 0, col = "gray")
 
plot(splinesBoot$functions[[2]], ylim = c(-1,2.5))
plot(splinesBoot$CI$alpha_0.05$lower[[2]], lty = 2, add = TRUE)
plot(splinesBoot$CI$alpha_0.05$upper[[2]], lty = 2, add = TRUE)
plot(splinesBoot$CI$alpha_0.1$lower[[2]], lty = 3, add = TRUE)
plot(splinesBoot$CI$alpha_0.1$upper[[2]], lty = 3, add = TRUE)
abline(h = 0, col = "gray")
legend("topleft", c("Estimate", "95% CI", "90% CI"), lty = 1:3, lwd = c(2,1,1))

# Plot 95% confidence bands for eigenvalues
plot(1:2, splinesBoot$values, pch = 20, ylim = c(0, 1.5), 
     main = "Estimated eigenvalues with 95% CI",
     xlab = "Eigenvalue no.", ylab = "")
arrows(1:2, splinesBoot$CIvalues$alpha_0.05$lower,
       1:2, splinesBoot$CIvalues$alpha_0.05$upper,
       length = 0.05, angle = 90, code = 3)
points(1:2, sim$trueVals[1:2], pch = 20, col = 4)
legend("topright", c("Estimate", "True value"), pch = 20, col = c(1,4))


### simulate data (two- and one-dimensional domains)
### ATTENTION: Takes long

set.seed(2)
sim <-  simMultiFunData(type = "weighted",
                 argvals = list(list(seq(0,1,0.01), seq(-1,1,0.02)), list(seq(-0.5,0.5,0.01))),
                 M = list(c(4,5), 20), eFunType = list(c("Fourier", "Fourier"), "Poly"),
                 eValType = "exponential", N = 150)

# MFPCA based on univariate spline expansions (for images) and univariate FPCA (for functions)
pca <- MFPCA(sim$simData, M = 10,
             uniExpansions = list(list(type = "splines2D", k = c(10,12)),
                             list(type = "uFPCA")))
summary(pca)
plot(pca) # plot the eigenfunctions as perturbations of the mean
scoreplot(pca) # plot the scores

### Compare to true eigenfunctions
# flip to make results more clear
pca$functions <- flipFuns(sim$trueFuns[1:10], pca$functions)

par(mfrow = c(5,2), mar = rep(2,4))
for(m in 2:6) # for m = 1, image.plot (used in plot(funData)) produces an error...
{
  plot(sim$trueFuns[[1]], main = paste("True, m = ", m), obs = m)
  plot(pca$functions[[1]], main = paste("Estimate, m = ", m), obs = m)
}

par(mfrow = c(1,1))
plot(sim$trueFuns[[2]], main = "Eigenfunctions (2nd element)", lwd = 2, obs=  1:5)
plot(pca$functions[[2]], lty = 2, add = TRUE, obs=  1:5)
legend("bottomleft", c("True", "MFPCA"), lty = 1:2, lwd = c(2,1))

par(oldPar)

Univariate functional principal component analysis by smoothed covariance

Description

This function calculates a univariate functional principal components analysis by smoothed covariance based on code from fpca.sc in package refund.

Usage

PACE(
  funDataObject,
  predData = NULL,
  nbasis = 10,
  pve = 0.99,
  npc = NULL,
  makePD = FALSE,
  cov.weight.type = "none"
)

Arguments

Value

Warning

This function works only for univariate functional data observed on one-dimensional domains.

See Also

funData, fpcaBasis, univDecomp

Examples

  oldPar <- par(no.readonly = TRUE)

  # simulate data
  sim <- simFunData(argvals = seq(-1,1,0.01), M = 5, eFunType = "Poly",
                    eValType = "exponential", N = 100)

  # calculate univariate FPCA
  pca <- PACE(sim$simData, npc = 5)

  # Plot the results
  par(mfrow = c(1,2))
  plot(sim$trueFuns, lwd = 2, main = "Eigenfunctions")
  # flip estimated functions for correct signs
  plot(flipFuns(sim$trueFuns,pca$functions), lty = 2, add = TRUE)
  legend("bottomright", c("True", "Estimate"), lwd = c(2,1), lty = c(1,2))

  plot(sim$simData, lwd = 2, main = "Some Observations", obs = 1:7)
  plot(pca$fit, lty = 2, obs = 1:7, add = TRUE) # estimates are almost equal to true values
  legend("bottomright", c("True", "Estimate"), lwd = c(2,1), lty = c(1,2))

  par(oldPar)

UMPCA: Uncorrelated Multilinear Principle Component Analysis

Description

This function implements the uncorrelated multilinear principal component analysis for tensors of dimension 2, 3 or 4. The code is basically the same as in the MATLAB toolbox UMPCA by Haiping Lu (Link: https://www.mathworks.com/matlabcentral/fileexchange/35432-uncorrelated-multilinear-principal-component-analysis-umpca, see also references).

Usage

UMPCA(TX, numP)

Arguments

Value

Warning

As this algorithm aims more at uncorrelated features than at an optimal reconstruction of the data, hence it might give poor results when used for the univariate decomposition of images in MFPCA.

References

Haiping Lu, K.N. Plataniotis, and A.N. Venetsanopoulos, "Uncorrelated Multilinear Principal Component Analysis for Unsupervised Multilinear Subspace Learning", IEEE Transactions on Neural Networks, Vol. 20, No. 11, Page: 1820-1836, Nov. 2009.

Examples

set.seed(12345)

 # define "true" components
 a <- sin(seq(-pi, pi, length.out = 100))
 b <- exp(seq(-0.5, 1, length.out = 150))

 # simulate tensor data
 X <- a %o% b %o% rnorm(80, sd = 0.5)

 # estimate one component
 UMPCAres <- UMPCA(X, numP = 1)

 # plot the results and compare to true values
 plot(UMPCAres$Us[[1]][,1])
 points(a/sqrt(sum(a^2)), pch = 20) # eigenvectors are defined only up to a sign change!
 legend("topright", legend = c("True", "Estimated"), pch = c(20, 1))

 plot(UMPCAres$Us[[2]][,1])
 points(b/sqrt(sum(b^2)), pch = 20)
 legend("topleft", legend = c("True", "Estimated"), pch = c(20, 1))

Utility function that calculates matrix of basis-scalar products (one dimension)

Description

If the element X^{(j)} is expanded in basis functions b_i^{(j)}(t),~ i = 1, \ldots, K_j, this function calculates the K_j \times K_j matrix B^{(jj)} with entries

B^{(jj)}_{mn} = \int_{\mathcal{T_j}} b_m^{(j)}(t) b_n^{(j)}(t) \mathrm{d} t

.

Usage

calcBasisIntegrals(basisFunctions, dimSupp, argvals)

Arguments

Value

A matrix containing the scalar product of all combinations of basis functions (matrix B^{(j)})

Warning

This function is implemented only for functions on one- or two-dimensional domains.

See Also

MFPCA, dimSupp

Internal function that implements the MFPCA algorithm for given univariate decompositions

Description

Internal function that implements the MFPCA algorithm for given univariate decompositions

Usage

calcMFPCA(
  N,
  p,
  Bchol,
  M,
  type,
  weights,
  npc,
  argvals,
  uniBasis,
  fit = FALSE,
  approx.eigen = FALSE
)

Arguments

Value

A list containing the following components:

Check input of FCP TPA function

Description

Check input of FCP TPA function

Usage

check_FCP_TPA_input(X, K, penMat, alphaRange, verbose, tol, maxIter, adaptTol)

Value

No return value, called for side effects

Calculate and threshold DCT for an image

Description

This function calculates the (orthonormal) discrete cosine transformation for an image and returns thresholded DCT coefficients using the C-library fftw3 (see http://www.fftw.org/).

Usage

dct2D(image, qThresh)

Arguments

Value

Warning

If the C-library fftw3 is not available when the package MFPCA is installed, this function is disabled an will throw an error. For full functionality install the C-library fftw3 from http://www.fftw.org/ and reinstall MFPCA. This function has not been tested with ATLAS/MKL/OpenBLAS.

See Also

dctBasis2D

Calculate and threshold DCT for an 3D image

Description

This function calculates the (orthonormal) discrete cosine transformation for a 3D image and returns thresholded DCT coefficients using the C-library fftw3 (see http://www.fftw.org/).

Usage

dct3D(image, qThresh)

Arguments

Value

Warning

If the C-library fftw3 is not available when the package MFPCA is installed, this function is disabled an will throw an error. For full functionality install the C-library fftw3 from http://www.fftw.org/ and reinstall MFPCA. This function has not been tested with ATLAS/MKL/OpenBLAS.

See Also

dctBasis3D

Calculate a cosine basis representation for functional data on two- or three-dimensional domains

Description

These functions calculate a tensor cosine basis representation for functional data on two- or three-dimensional domains based on a discrete cosine transformation (DCT) using the C-library fftw3 (http://www.fftw.org/). Coefficients under a given threshold are set to 0 to reduce complexity and for denoising.

Usage

dctBasis2D(funDataObject, qThresh, parallel = FALSE)

dctBasis3D(funDataObject, qThresh, parallel = FALSE)

Arguments

Details

Given the (discretized) observed functions X_i, the function dctBasis2D calculates a basis representation

X_i(s,t) = \sum_{m = 0}^{K_1-1} \sum_{n = 0}^{K_2-1} \theta_{mn} f_{mn}(s,t)

of a two-dimensional function X_i(s,t) in terms of (orthogonal) tensor cosine basis functions

f_{mn}(s,t) = c_m c_n \cos(ms) \cos(nt), \quad (s,t) \in \mathcal{T}

with c_m = \frac{1}{\sqrt{\pi}} for m=0 and c_m = \sqrt{\frac{2}{\pi}} for m=1,2,\ldots based on a discrete cosine transform (DCT).

If not thresholded (qThresh = 0), the function returns all non-zero coefficients \theta_{mn} in the basis representation in a sparseMatrix (package Matrix) called scores. Otherwise, coefficients with

|\theta_{mn}| <= q

are set to zero, where q is the qThresh-quantile of |\theta_{mn}|.

For functions X_i(s,t,u) on three-dimensional domains, the function dctBasis3D calculates a basis representation

X_i(s,t,u) = \sum_{m = 0}^{K_1-1} \sum_{n = 0}^{K_2-1} \sum_{k = 0}^{K_3-1} \theta_{mnk} f_{mnk}(s,t,u)

in terms of (orthogonal) tensor cosine basis functions

f_{mnk}(s,t,u) = c_m c_n c_k \cos(ms) \cos(nt) \cos(ku), \quad (s,t,u) \in \mathcal{T}

again with c_m = \frac{1}{\sqrt{pi}} for m=0 and c_m = \sqrt{\frac{2}{pi}} for m=1,2,\ldots based on a discrete cosine transform (DCT). The thresholding works analogous as for the two-dimensional case.

Value

Warning

If the C-library fftw3 is not available when the package MFPCA is installed, this function is disabled an will throw an error. For full functionality install the C-library fftw3 from http://www.fftw.org/ and reinstall MFPCA. This function has not been tested with ATLAS/MKL/OpenBLAS.

See Also

univDecomp, dct2D, dct3D

Calculate linear combinations of orthonormal cosine basis functions on two- or three-dimensional domains

Description

Given scores (coefficients), these functions calculate a linear combination of two- or three-dimensional cosine tensor basis functions on two- or three-dimensional domains using the C-library fftw3 (see http://www.fftw.org/).

Usage

dctFunction2D(scores, argvals, parallel = FALSE)

dctFunction3D(scores, argvals, parallel = FALSE)

Arguments

Value

An object of class funData with N observations on the two- or threedimensional domain specified by argvals, corresponding to the linear combination of orthonormal cosine basis functions.

Warning

If the C-library fftw3 is not available when the package MFPCA is installed, the functions are disabled an will throw an error. For full functionality install the C-library fftw3 from http://www.fftw.org/ and reinstall MFPCA. This function has not been tested with ATLAS/MKL/OpenBLAS.

See Also

univExpansion, idct2D, idct3D, dctBasis2D, dctBasis3D

Calculate a linear combination of arbitrary basis function

Description

This function calculates a linear combination of arbitrary basis functions on domains with arbitrary dimension.

Usage

expandBasisFunction(scores, argvals = functions@argvals, functions)

Arguments

Value

An object of class funData with N observations on argvals, corresponding to the linear combination of the basis functions.

See Also

univExpansion

Calculate a smooth PCA representation for functional data on two-dimensional domains

Description

This function calculates a smooth PCA representation based on the FCP_TPA algorithm (see References) for functional data on two-dimensional domains. In this case, the data can be interpreted as images with S1 x S2 pixels (assuming nObsPoints(funDataObject) = (S1, S2)), i.e. the total data for N observations can be represented as third order tensor of dimension N x S1 x S2.

Usage

fcptpaBasis(
  funDataObject,
  npc,
  smoothingDegree = rep(2, 2),
  alphaRange,
  orderValues = TRUE,
  normalize = FALSE
)

Arguments

Details

The smooth PCA of the tensor data is calculated via the FCP_TPA function. Smoothness is induced by difference penalty matrices for both directions of the images, weighted by smoothing parameters \alpha_v, \alpha_w. The resulting eigenvectors can be interpreted in terms of eigenfunctions and individual scores for each observation. See FCP_TPA for details.

Value

References

G. I. Allen, "Multi-way Functional Principal Components Analysis", In IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing, 2013.

See Also

univDecomp, FCP_TPA

Use a basis from package fda for univariate representation

Description

This function allows to use univariate basis representations from the fda package using the funData2fd function from package funData.

Usage

fdaBasis(funDataObject, ...)

Arguments

Value

Warning

The package fda has be be installed to use this functionality.

See Also

funData2fd, eval.fd

Find the optimal smoothing parameters in FCP_TPA using GCV

Description

These functions find the optimal smoothing parameters \alpha_v, \alpha_w for the two image directions (v and w) in the FCP_TPA algorithm based on generalized cross-validation, which is nested in the tensor power algorithm. Given a range of possible values of \alpha_v (or \alpha_w, respectively), the optimum is found by optimizing the GCV criterion using the function optimize.

Usage

findAlphaVopt(alphaRange, data, u, w, alphaW, OmegaW, GammaV, lambdaV)

findAlphaWopt(alphaRange, data, u, v, alphaV, OmegaV, GammaW, lambdaW)

Arguments

Value

The optimal \alpha_v (or \alpha_w, respectively), found by optimizing the GCV criterion within the given range of possible values.

Functions

• findAlphaWopt:

References

G. I. Allen (2013), "Multi-way Functional Principal Components Analysis", IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing.

J. Z. Huang, H. Shen and A. Buja (2009), "The Analysis of Two-Way Functional Data Using Two-Way Regularized Singular Value Decomposition". Journal of the American Statistical Association, Vol. 104, No. 488, 1609 – 1620.

See Also

FCP_TPA, gcv

Calculate a functional principal component basis representation for functional data on one-dimensional domains

Description

This function calculates a functional principal component basis representation for functional data on one-dimensional domains. The FPCA is calculated via the PACE function, which is built on fpca.sc in the refund package.

Usage

fpcaBasis(
  funDataObject,
  nbasis = 10,
  pve = 0.99,
  npc = NULL,
  makePD = FALSE,
  cov.weight.type = "none"
)

Arguments

Value

See Also

univDecomp, PACE

Generalized cross-validation for the FCP-TPA algorithm

Description

These function calculates the generalized cross-validation criterion for the smoothing parameters \alpha_v or \alpha_w that are used in the FCP_TPA algorithm. As the criterion is symmetric in v and w, this function implements a generic criterion, which is called by findAlphaVopt, findAlphaWopt with the correct values.

Usage

gcv(alpha, n, z, eta, lambda)

Arguments

Details

The criterion can be evaluated in a numerically efficient way, adopting the ideas in Huang, Shen and Buja (2008) to three-ways tensors. TODO!

Value

The value of the GCV criterion.

References

G. I. Allen, "Multi-way Functional Principal Components Analysis", IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing, 2013.

See Also

FCP_TPA

Use given basis functions for univariate representation

Description

Use given basis functions for univariate representation

Usage

givenBasis(funDataObject, functions, scores = NULL, ortho = NULL)

Arguments

Value

Calculate an inverse DCT for an image

Description

This function calculates an inverse (orthonormal) discrete cosine transformation for given coefficients in two dimensions using the C-library fftw3 (see http://www.fftw.org/). As many coefficients are expected to be zero, the values are given in compressed format (indices and values only of non-zero coefficients).

Usage

idct2D(scores, ind, dim)

Arguments

Value

A matrix of dimensions dim, which is a linear combination of cosine tensor basis functions with the given coefficients.

Warning

If the C-library fftw3 is not available when the package MFPCA is installed, this function is disabled an will throw an error. For full functionality install the C-library fftw3 from http://www.fftw.org/ and reinstall MFPCA. This function has not been tested with ATLAS/MKL/OpenBLAS.

See Also

dctBasis2D

Calculate an inverse DCT for a 3D image

Description

This function calculates an inverse (orthonormal) discrete cosine transformation for given coefficients in three dimensions using the C-library fftw3 (see http://www.fftw.org/). As many coefficients are expected to be zero, the values are given in compressed format (indices and values only of non-zero coefficients).

Usage

idct3D(scores, ind, dim)

Arguments

Value

A matrix of dimensions dim, which is a linear combination of cosine tensor basis functions with the given coefficients.

Warning

If the C-library fftw3 is not available when the package MFPCA is installed, this function is disabled an will throw an error. For full functionality install the C-library fftw3 from http://www.fftw.org/ and reinstall MFPCA. This function has not been tested with ATLAS/MKL/OpenBLAS.

See Also

dctBasis3D

Compute the largest eigenvalue and associated eigenvector of a matrix A using the power method

Description

MATLAB Function: Copyright 1999 by Todd K. Moon

Usage

maxeig(A)

Arguments

Value

Calculate multivariate basis expansion

Description

Calculate multivariate basis expansion

Usage

multivExpansion(multiFuns, scores)

Arguments

Value

A multiFunData object containing the expanded functions for each observation.

Calculate the euclidean norm of a vector

Description

Calculate the euclidean norm of a vector

Usage

normVec(x)

Arguments

Value

The euclidean norm of x.

Plot MFPCA results

Description

Plots the eigenfunctions as perturbations of the mean (i.e. the mean function plus/minus a constant factor times each eigenfunction separately). If all elements have a one-dimensional domain, the plots can be combined, otherwise the effects of adding and subtracting are shown in two separate rows for each eigenfunction.

Usage

## S3 method for class 'MFPCAfit'
plot(
  x,
  plotPCs = seq_len(nObs(x$functions)),
  stretchFactor = NULL,
  combined = FALSE,
  ...
)

Arguments

Value

A plot of the principal components as perturbations of the mean.

See Also

MFPCA, plot.funData

Examples

# Simulate multivariate functional data on one-dimensonal domains
# and calculate MFPCA (cf. MFPCA help)
set.seed(1)
# simulate data (one-dimensional domains)
sim <-  simMultiFunData(type = "split", argvals = list(seq(0,1,0.01), seq(-0.5,0.5,0.02)),
                       M = 5, eFunType = "Poly", eValType = "linear", N = 100)
# MFPCA based on univariate FPCA
PCA <- MFPCA(sim$simData, M = 5, uniExpansions = list(list(type = "uFPCA"),
                                                     list(type = "uFPCA")))

# Plot the results
plot(PCA, combined = TRUE) # combine addition and subtraction in one plot

Function prediction based on MFPCA results

Description

Predict functions based on a truncated multivariate Karhunen-Loeve representation:

\hat x = \hat mu + \sum_{m = 1}^M \rho_m \hat \psi_m

with estimated mean function \hat \mu and principal components \psi_m. The scores \rho_m can be either estimated (reconstruction of observed functions) or user-defined (construction of new functions).

Usage

## S3 method for class 'MFPCAfit'
predict(object, scores = object$scores, ...)

Arguments

Value

A multiFunData object containing the predicted functions.

See Also

MFPCA

Examples

#' # Simulate multivariate functional data on one-dimensonal domains
# and calculate MFPCA (cf. MFPCA help)
set.seed(1)
# simulate data (one-dimensional domains)
sim <-  simMultiFunData(type = "split", argvals = list(seq(0,1,0.01), seq(-0.5,0.5,0.02)),
                       M = 5, eFunType = "Poly", eValType = "linear", N = 100)
# MFPCA based on univariate FPCA
PCA <- MFPCA(sim$simData, M = 5, uniExpansions = list(list(type = "uFPCA"),
                                                     list(type = "uFPCA")))

# Reconstruct the original data
pred <- predict(PCA) # default reconstructs data used for the MFPCA fit

# plot the results: 1st element
plot(sim$simData[[1]]) # original data
plot(pred[[1]], add = TRUE, lty = 2) # reconstruction

# plot the results: 2nd element
plot(sim$simData[[2]]) # original data
plot(pred[[2]], add = TRUE, lty = 2) # reconstruction

Print the results of a Multivariate Functional Principal Component Analysis

Description

A print function for class MFPCAfit.

Usage

## S3 method for class 'MFPCAfit'
print(x, ...)

Arguments

Value

No return value, called for side effects

Print summary of a Multivariate Functional Principal Component Analysis

Description

A print method for class MFPCAfit.summary

Usage

## S3 method for class 'summary.MFPCAfit'
print(x, ...)

Arguments

Value

No return value, called for side effects

Scoreplot Generic

Description

Redirects to plot.default

Usage

scoreplot(PCAobject, ...)

Arguments

Value

A bivariate plot of scores.

Plot the Scores of a Multivariate Functional Principal Component Analysis

Description

This function plots two scores of a multivariate functional principal component analysis for each observation.

Usage

## S3 method for class 'MFPCAfit'
scoreplot(PCAobject, choices = 1:2, scale = FALSE, ...)

Arguments

Value

A bivariate plot of scores.

See Also

MFPCA

Examples

# and calculate MFPCA (cf. MFPCA help)
set.seed(1)
# simulate data (one-dimensional domains)
sim <-  simMultiFunData(type = "split", argvals = list(seq(0,1,0.01), seq(-0.5,0.5,0.02)),
                       M = 5, eFunType = "Poly", eValType = "linear", N = 100)
# MFPCA based on univariate FPCA
PCA <- MFPCA(sim$simData, M = 5, uniExpansions = list(list(type = "uFPCA"),
                                                     list(type = "uFPCA")))

# Plot the first two scores
scoreplot(PCA) # no scaling (default)
scoreplot(PCA, scale = TRUE) # scale the scores by the first two eigenvalues 

Screeplot for Multivariate Functional Principal Component Analysis

Description

This function plots the proportion of variance explained by the leading eigenvalues in an MFPCA against the number of the principal component.

Usage

## S3 method for class 'MFPCAfit'
screeplot(
  x,
  npcs = min(10, length(x$values)),
  type = "lines",
  ylim = NULL,
  main = deparse(substitute(x)),
  ...
)

Arguments

Value

A screeplot, showing the decrease of the principal component score.

See Also

MFPCA, screeplot

Examples

# Simulate multivariate functional data on one-dimensonal domains
# and calculate MFPCA (cf. MFPCA help)
set.seed(1)
# simulate data (one-dimensional domains)
sim <-  simMultiFunData(type = "split", argvals = list(seq(0,1,0.01), seq(-0.5,0.5,0.02)),
                       M = 5, eFunType = "Poly", eValType = "linear", N = 100)
# MFPCA based on univariate FPCA
PCA <- MFPCA(sim$simData, M = 5, uniExpansions = list(list(type = "uFPCA"),
                                                     list(type = "uFPCA")))

# screeplot
screeplot(PCA) # default options
screeplot(PCA, npcs = 3, type = "barplot", main= "Screeplot")

Calculate a spline basis decomposition for functional data on one-dimensional domains

Description

These functions calculate a penalized or unpenalized spline basis decomposition for functional data on one-dimensional domains based on the gam function in the mgcv package.

Usage

splineBasis1D(funDataObject, bs = "ps", m = NA, k = -1)

splineBasis1Dpen(funDataObject, bs = "ps", m = NA, k = -1, parallel = FALSE)

Arguments

Value

See Also

univDecomp, gam, foreach

Calculate a spline basis representation for functional data on two-dimensional domains

Description

These functions calculate a penalized or unpenalized tensor product spline basis representation for functional data on two-dimensional domains based on the gam/bam functions in the mgcv package. See Details.

Usage

splineBasis2D(funDataObject, bs = "ps", m = NA, k = -1)

splineBasis2Dpen(funDataObject, bs = "ps", m = NA, k = -1, parallel = FALSE)

Arguments

Details

If the basis representation is calculated without penalization (splineBasis2D), the coefficients are computed using the gam function from the mgcv package. In the case of penalization (splineBasis2Dpen), the function bam (for large GAMs) is used instead.

Value

See Also

univDecomp, splineBasis1D, gam, bam, foreach

Calculate linear combinations of spline basis functions on one-dimensional domains

Description

Given scores (coefficients), this function calculates a linear combination of spline basis functions on one-dimensional domains based on the gam function in the mgcv package.

Usage

splineFunction1D(scores, argvals, bs, m, k)

Arguments

Value

An object of class funData with N observations on argvals, corresponding to the linear combination of spline basis functions.

See Also

univExpansion, gam, splineBasis1D

Calculate linear combinations of spline basis functions on two-dimensional domains

Description

Given scores (coefficients), these functions calculate a linear combination of spline tensor basis functions on two-dimensional domains based on the gam/bam functions in the mgcv package. See Details.

Usage

splineFunction2D(scores, argvals, bs, m, k)

splineFunction2Dpen(scores, argvals, bs, m, k)

Arguments

Details

If the scores have been calculated based on an unpenalized tensor spline basis, the linear combination is computed based on the gam functions ((splineFunction2D)). If the scores were obtained using penalization, the expansion is calculated via bam (splineFunction2Dpen).

Value

An object of class funData with N observations on the two-dimensional domain specified by argvals, corresponding to the linear combination of spline basis functions.

Warning

The function splineFunction2Dpen, which relies on bam has not been tested with ATLAS/MKL/OpenBLAS.

See Also

univExpansion, gam, splineBasis2D

Sample stratified indices according to a factor variable

Description

Sample stratified indices according to a factor variable

Usage

stratSample(f)

Arguments

Value

A vector of the same length as f, containing the resampled indices, stratified according to the levels of f

Summarize a Multivariate Functional Principal Component Analysis

Description

A summary method for class MFPCAfit

Usage

## S3 method for class 'MFPCAfit'
summary(object, ...)

Arguments

Value

An object of class summary.MFPCAfit

Tensor times vector calculation

Description

Functionality adapted from the MATLAB tensor toolbox (https://www.tensortoolbox.org/).

Usage

ttv(A, v, dim)

Arguments

Details

Let A be a tensor with dimensions d_1 \times d_2 \times \ldots \times d_p and let v be a vector of length d_i. Then the tensor-vector-product along the i-th dimension is defined as

B_{j_1 \ldots j_{i-1}j_{i+1} \ldots j_d} = \sum_{i=1}^{d_i} A_{j_1 \ldots j_{i-1} i j_{i+1} \ldots j_d} \cdot v_i.

It can hence be seen as a generalization of the matrix-vector product.

The tensor-vector-product along several dimensions between a tensor A and multiple vectors v_1,...,v_k (k \le p) is defined as a series of consecutive tensor-vector-product along the different dimensions. For consistency, the multiplications are calculated from the dimension of the highest order to the lowest.

Value

An array, the result of the multiplication.

References

B. W. Bader and T. G. Kolda. Algorithm 862: MATLAB tensor classes for fast algorithm prototyping, ACM Transactions on Mathematical Software 32(4):635-653, December 2006.

See Also

UMPCA

Examples

# create a three-mode tensor
a1 <- seq(0,1, length.out = 10)
a2 <- seq(-1,1, length.out = 20)
a3 <- seq(-pi, pi, length.out = 15)
A <-a1 %o% a2 %o% a3
dim(A)

# multiply along different dimensions
dim(ttv(A = A, v = list(rnorm(10)), dim = 1))
dim(ttv(A = A, v = list(rnorm(20)), dim = 2))
dim(ttv(A = A, v = list(rnorm(15)), dim = 3))

# multiply along more than one dimension
length(ttv(A = A, v = list(rnorm(10), rnorm(15)), dim = c(1,3)))

Internal function for the Tensor times Vector calculation

Description

Internal function for the Tensor times Vector calculation

Usage

ttvCalculation(A, v, dim)

Arguments

Value

Result of tensor times vector calculation

See Also

ttv

Calculate an uncorrelated multilinear principal component basis representation for functional data on two-dimensional domains

Description

This function calculates an uncorrelated multilinear principal component analysis (UMPCA) representation for functional data on two-dimensional domains. In this case, the data can be interpreted as images with S1 x S2 pixels (assuming nObsPoints(funDataObject) = (S1, S2)), i.e. the total observed data are represented as third order tensor of dimension N x S1 x S2. The UMPCA of a tensor of this kind is calculated via the UMPCA function, which is an R-version of the analogous functions in the UMPCA MATLAB toolbox by Haiping Lu (Link: https://www.mathworks.com/matlabcentral/fileexchange/35432-uncorrelated-multilinear-principal-component-analysis-umpca, see also references).

Usage

umpcaBasis(funDataObject, npc)

Arguments

Value

Warning

As this algorithm aims more at uncorrelated features than at an optimal reconstruction of the data, hence it might give poor results when used for the univariate decomposition of images in MFPCA. The function therefore throws a warning.

References

Haiping Lu, K.N. Plataniotis, and A.N. Venetsanopoulos, "Uncorrelated Multilinear Principal Component Analysis for Unsupervised Multilinear Subspace Learning", IEEE Transactions on Neural Networks, Vol. 20, No. 11, Page: 1820-1836, Nov. 2009.

See Also

univDecomp

Univariate basis decomposition

Description

This function calculates a univariate basis decomposition for a (univariate) functional data object.

Usage

univDecomp(type, funDataObject, ...)

Arguments

Details

Functional data X_i(t) can often be approximated by a linear combination of basis functions b_k(t)

X_i(t) = \sum_{k = 1}^K \theta_{ik} b_k(t), i = 1, \ldots, N.

The basis functions may be prespecified (such as spline basis functions or Fourier bases) or can be estimated from the data (e.g. by functional principal component analysis) and are the same for all observations X_1(t), \ldots, X_n(t). The coefficients (or scores) \theta_{ik} reflect the weight of each basis function b_k(t) for the observed function X_i(t) and can be used to characterize the individual observations.

Value

Warning

The options type = "DCT2D" and type = "DCT3D" have not been tested with ATLAS/MKL/OpenBLAS.

See Also

MFPCA, univExpansion, fpcaBasis, splineBasis1D, splineBasis1Dpen, splineBasis2D, splineBasis2Dpen, umpcaBasis, fcptpaBasis, fdaBasis, dctBasis2D, dctBasis3D

Examples

# generate some data
dat <- simFunData(argvals = seq(0,1,0.01), M = 5,
                  eFunType = "Poly", eValType = "linear", N = 100)$simData

# decompose the data in univariate functional principal components...
decFPCA <- univDecomp(type = "uFPCA", funDataObject = dat, npc = 5)
str(decFPCA)

# or in splines (penalized)
decSplines <- univDecomp(type = "splines1Dpen", funDataObject = dat) # use mgcv's default params
str(decSplines)

Calculate a univariate basis expansion

Description

This function calculates a univariate basis expansion based on given scores (coefficients) and basis functions.

Usage

univExpansion(
  type,
  scores,
  argvals = ifelse(!is.null(functions), functions@argvals, NULL),
  functions,
  params = NULL
)

Arguments

Details

This function calculates functional data X_i(t), i= 1 \ldots N that is represented as a linear combination of basis functions b_k(t)

X_i(t) = \sum_{k = 1}^K \theta_{ik} b_k(t), i = 1, \ldots, N.

The basis functions may be prespecified (such as spline basis functions or Fourier bases) or can be estimated from observed data (e.g. by functional principal component analysis). If type = "default" (i.e. a linear combination of arbitrary basis functions is to be calculated), both scores and basis functions must be supplied.

Value

An object of class funData with N observations on argvals, corresponding to the linear combination of the basis functions.

Warning

The options type = "spline2Dpen", type = "DCT2D" and type = "DCT3D" have not been tested with ATLAS/MKL/OpenBLAS.

See Also

MFPCA, splineFunction1D, splineFunction2D, splineFunction2Dpen, dctFunction2D, dctFunction3D, expandBasisFunction

Examples

oldPar <- par(no.readonly = TRUE)
par(mfrow = c(1,1))

set.seed(1234)

### Spline basis ###
# simulate coefficients (scores) for N = 10 observations and K = 8 basis functions
N <- 10
K <- 8
scores <- t(replicate(n = N, rnorm(K, sd = (K:1)/K)))
dim(scores)

# expand spline basis on [0,1]
funs <- univExpansion(type = "splines1D", scores = scores, argvals = list(seq(0,1,0.01)),
                      functions = NULL, # spline functions are known, need not be given
                      params = list(bs = "ps", m = 2, k = K)) # params for mgcv

plot(funs, main = "Spline reconstruction")

### PCA basis ###
# simulate coefficients (scores) for N = 10 observations and K = 8 basis functions
N <- 10
K <- 8

scores <- t(replicate(n = N, rnorm(K, sd = (K:1)/K)))
dim(scores)

# Fourier basis functions as eigenfunctions
eFuns <- eFun(argvals = seq(0,1,0.01), M = K, type = "Fourier")

# expand eigenfunction basis
funs <-  univExpansion(type = "uFPCA", scores = scores,
                       argvals = NULL, # use argvals of eFuns (default)
                       functions = eFuns)

plot(funs, main = "PCA reconstruction")

par(oldPar)