Fdapace Manual

User Manual:

Open the PDF directly: View PDF PDF.
Page Count: 73

DownloadFdapace-manual
Open PDF In BrowserView PDF
Package ‘fdapace’
April 24, 2019
Type Package
Title Functional Data Analysis and Empirical Dynamics
URL https://github.com/functionaldata/tPACE
BugReports https://github.com/functionaldata/tPACE/issues
Version 0.5.0
Date 2019-03-20,
Maintainer Pantelis Z. Hadjipantelis 
Description Provides implementation of various methods of Functional Data
Analysis (FDA) and Empirical Dynamics. The core of this package is Functional
Principal Component Analysis (FPCA), a key technique for functional data
analysis, for sparsely or densely sampled random trajectories and time courses,
via the Principal Analysis by Conditional Estimation (PACE) algorithm or
numerical integration. PACE is useful for the analysis of data that have been
generated by a sample of underlying (but usually not fully observed) random
trajectories. It does not rely on pre-smoothing of trajectories, which is
problematic if functional data are sparsely sampled. PACE provides options
for functional regression and correlation, for Longitudinal Data Analysis,
the analysis of stochastic processes from samples of realized trajectories,
and for the analysis of underlying dynamics. The core computational algorithms
are implemented using the 'Eigen' C++ library for numerical linear algebra and
'RcppEigen' ``glue''.
License BSD_3_clause + file LICENSE
LazyData false
Imports Rcpp (>= 0.11.5), Hmisc, MASS, Matrix, pracma, numDeriv
LinkingTo Rcpp, RcppEigen
Suggests plot3D, rgl, aplpack, mgcv, ks, gtools, knitr, EMCluster,
minqa, testthat
NeedsCompilation yes
RoxygenNote 6.1.1
VignetteBuilder knitr
Author Xiongtao Dai [aut],
Pantelis Z. Hadjipantelis [aut, cre],
Kynghee Han [aut],
Hao Ji [aut],
1

R topics documented:

2
Shu-Chin Lin [ctb],
Hans-Georg Mueller [cph, ths],
Jane-Ling Wang [cph, ths]

R topics documented:
BwNN . . . . . . . . . .
CheckData . . . . . . . .
CheckOptions . . . . . .
ConvertSupport . . . . .
CreateBasis . . . . . . .
CreateBWPlot . . . . . .
CreateCovPlot . . . . . .
CreateDesignPlot . . . .
CreateDiagnosticsPlot . .
CreateFuncBoxPlot . . .
CreateModeOfVarPlot .
CreateOutliersPlot . . . .
CreatePathPlot . . . . .
CreateScreePlot . . . . .
CreateStringingPlot . . .
cumtrapzRcpp . . . . . .
DynCorr . . . . . . . . .
Dyn_test . . . . . . . . .
FAM . . . . . . . . . . .
FCCor . . . . . . . . . .
FClust . . . . . . . . . .
FCReg . . . . . . . . . .
fdapace . . . . . . . . .
fitted.FPCA . . . . . . .
fitted.FPCAder . . . . .
FOptDes . . . . . . . . .
FPCA . . . . . . . . . .
FPCAder . . . . . . . .
FPCquantile . . . . . . .
FPCReg . . . . . . . . .
FSVD . . . . . . . . . .
FVPA . . . . . . . . . .
GetCrCorYX . . . . . .
GetCrCorYZ . . . . . .
GetCrCovYX . . . . . .
GetCrCovYZ . . . . . .
GetNormalisedSample .
kCFC . . . . . . . . . .
Lwls1D . . . . . . . . .
Lwls2D . . . . . . . . .
Lwls2DDeriv . . . . . .
MakeBWtoZscore02y . .
MakeFPCAInputs . . . .
MakeGPFunctionalData
MakeHCtoZscore02y . .
MakeLNtoZscore02y . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

3
4
4
5
5
6
6
7
8
9
10
11
12
13
14
14
15
16
17
20
21
23
24
25
27
28
29
32
34
35
38
40
41
42
42
43
45
46
47
47
48
49
50
50
51
51

BwNN

3
MakeSparseGP .
medfly25 . . . .
MultiFAM . . . .
NormCurvToArea
predict.FPCA . .
print.FPCA . . .
print.FSVD . . .
print.WFDA . . .
SBFitting . . . .
SelectK . . . . .
SetOptions . . . .
Sparsify . . . . .
str.FPCA . . . . .
Stringing . . . .
trapzRcpp . . . .
TVAM . . . . . .
VCAM . . . . .
WFDA . . . . . .
Wiener . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

Index

BwNN

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

52
52
53
57
57
58
59
59
59
62
63
63
64
64
65
66
68
70
72
73

Minimum bandwidth based on kNN criterion.

Description
Input a list of time points Lt, and the number of unique neighbors k and get the minimum bandwidth
garanteeing k unique neighbours.
Usage
BwNN(Lt, k = 3, onlyMean = FALSE, onlyCov = FALSE)
Arguments
Lt

n-by-1 list of vectors

k

number of unique neighbors for cov and mu (default = 3)

onlyMean

Indicator to return only the minimum bandwidth for the mean

onlyCov

Indicator to return only the minium bandwidth for the covariance

Examples
tinyGrid = list(c(1,7), c(2,3),
BwNN(tinyGrid, k = 2) # c(3,2)

6,

c(2,4), c(4,5))

4

CheckOptions

CheckData

Check data format

Description
Check if there are problems with the form and basic structure of the functional data ’y’ and the
recorded times ’t’.

Usage
CheckData(y, t)

Arguments
y

is a n-by-1 list of vectors

t

is a n-by-1 list of vectors

CheckOptions

Check option format

Description
Check if the options structure is valid and set the ones that are NULL

Usage
CheckOptions(t, optns, n)

Arguments
t

is a n-by-1 list of vectors

optns

is an initialized option list

n

is a total number of sample curves

ConvertSupport

ConvertSupport

5

Convert support of a mu/phi/cov etc. to and from obsGrid and workGrid

Description
Convert the support of a given function 1-D or 2-D function from ’fromGrd’ to ’toGrid’. Both grids
need to be sorted. This is a interpolation/convenience function.
Usage
ConvertSupport(fromGrid, toGrid, mu = NULL, Cov = NULL, phi = NULL,
isCrossCov = FALSE)
Arguments
fromGrid

vector of points with input grid to interpolate from

toGrid

vector of points with the target grid to interpolate on

mu

any vector of function to be interpolated

Cov

a square matrix supported on fromGrid * fromGrid, to be interpolated to toGrid
* toGrid.

phi

any matrix, each column containing a function to be interpolated

isCrossCov

logical, indicating whether the input covariance is a cross-covariance. If so then
the output is not made symmetric.

CreateBasis

Create an orthogonal basis of K functions in [0, 1], with nGrid points.

Description
Create an orthogonal basis of K functions in [0, 1], with nGrid points.
Usage
CreateBasis(K, pts = seq(0, 1, length.out = 50), type = c("cos", "sin",
"fourier", "legendre01", "poly"))
Arguments
K

A positive integer specifying the number of eigenfunctions to generate.

pts

A vector specifying the time points to evaluate the basis functions.

type

A string for the type of orthogonal basis.

Value
A K by nGrid matrix, each column containing an basis function.

6

CreateCovPlot

Examples
basis <- CreateBasis(3, type='fourier')
head(basis)

CreateBWPlot

Functional Principal Component Analysis Bandwidth Diagnostics
plot

Description
This function by default creates the mean and first principal modes of variation plots for 50 If
provided with a derivative options object (?FPCAder) it will return the differentiated mean and first
two principal modes of variations for 50
Usage
CreateBWPlot(fpcaObj, derOptns = NULL, bwMultipliers = NULL)
Arguments
fpcaObj

An FPCA class object returned by FPCA().

derOptns

A list of options to control the derivation parameters; see ?FPCAder. If NULL
standard diagnostics are returned

bwMultipliers

A vector of multipliers that the original ’bwMu’ and ’bwCov’ will be multiplied
by. (default: c(0.50, 0.75, 1.00, 1.25, 1.50)) - default: NULL

Examples
set.seed(1)
n <- 25
pts <- seq(0, 1, by=0.05)
sampWiener <- Wiener(n, pts)
sampWiener <- Sparsify(sampWiener, pts, 10)
res1 <- FPCA(sampWiener$Ly, sampWiener$Lt,
list(dataType='Sparse', error=FALSE, kernel='epan', verbose=FALSE))
CreateBWPlot(res1)

CreateCovPlot

Create the covariance surface plot based on the results from FPCA()
or FPCder().

Description
This function will open a new device if not instructed otherwise.
Usage
CreateCovPlot(fpcaObj, covPlotType = "Fitted", isInteractive = FALSE,
colSpectrum = NULL, ...)

CreateDesignPlot

7

Arguments
fpcaObj

returned object from FPCA().

covPlotType

a string specifying the type of covariance surface to be plotted: ’Smoothed’: plot
the smoothed cov surface ’Fitted’: plot the fitted cov surface

isInteractive

an option for interactive plot: TRUE: interactive plot; FALSE: printable plot

colSpectrum

character vector to be use as input in the ’colorRampPalette’ function defining
the colouring scheme (default: c(’blue’,’red’))

...

other arguments passed into persp3d, persp3D, plot3d or points3D for plotting
options

Examples
set.seed(1)
n <- 20
pts <- seq(0, 1, by=0.05)
sampWiener <- Wiener(n, pts)
sampWiener <- Sparsify(sampWiener, pts, 10)
res <- FPCA(sampWiener$Ly, sampWiener$Lt,
list(dataType='Sparse', error=FALSE, kernel='epan', verbose=TRUE))
CreateCovPlot(res)

CreateDesignPlot

Create the design plot of the functional data.

Description
Display the scatter plot of {(Tij , Til }i≤n, j,l≤ni , which is used to diagnose whether the design time
points are sufficiently dense in the domain of interest so that the 2D smoothed covariance estimate
is feasible. See Figure 2 of Yao et al (2005). This function will open a new device if not instructed
otherwise.
Usage
CreateDesignPlot(Lt, obsGrid = NULL, isColorPlot = TRUE,
noDiagonal = TRUE, addLegend = TRUE, ...)
Arguments
Lt

a list of observed time points for functional data

obsGrid

a vector of sorted observed time points. Default to the unique time points in Lt.

isColorPlot

an option for colorful plot: TRUE: create color plot with color indicating counts
FALSE: create black and white plot with dots indicating observed time pairs

noDiagonal

an option specifying plotting the diagonal design points: TRUE: remove diagonal time pairs FALSE: do not remove diagonal time pairs

addLegend

Logical, default TRUE

...

Other arguments passed into plot().

8

CreateDiagnosticsPlot

References
Yao, Fang, Hans-Georg Mueller, and Jane-Ling Wang. "Functional data analysis for sparse longitudinal data." Journal of the American Statistical Association 100, no. 470 (2005): 577-590.
Examples
set.seed(1)
n <- 20
pts <- seq(0, 1, by=0.05)
sampWiener <- Wiener(n, pts)
sampWiener <- Sparsify(sampWiener, pts, 10)
CreateDesignPlot(sampWiener$Lt, sort(unique(unlist(sampWiener$Lt))))

CreateDiagnosticsPlot Functional Principal Component Analysis Diagnostics plot

Description
Deprecated. Use plot.FPCA instead.
This function plot the results for an FPCA. It prints the design plot, mean function, scree-plot and
the first three eigenfunctions of a sample. If provided with a derivative options object (?FPCAder) it
will return the differentiated mean and first two principal modes of variations for 50%, 75%, 100%,
125% and 150% of the defined bandwidth choice.
Usage
CreateDiagnosticsPlot(...)
## S3 method for class 'FPCA'
plot(x, openNewDev = FALSE, addLegend = TRUE, ...)
Arguments
...

passed into plot.FPCA.

x

An FPCA class object returned by FPCA().

openNewDev

A logical specifying if a new device should be opened - default: FALSE

addLegend

A logical specifying whether to add legend.

Details
The black, red, and green curves stand for the first, second, and third eigenfunctions, respectively.
plot.FPCA is currently implemented only for the original function, but not a derivative FPCA object.

CreateFuncBoxPlot

9

Examples
set.seed(1)
n <- 20
pts <- seq(0, 1, by=0.05)
sampWiener <- Wiener(n, pts)
sampWiener <- Sparsify(sampWiener, pts, 10)
res1 <- FPCA(sampWiener$Ly, sampWiener$Lt,
list(dataType='Sparse', error=FALSE, kernel='epan', verbose=FALSE))
plot(res1)

CreateFuncBoxPlot

Create functional boxplot using ’bagplot’, ’KDE’ or ’pointwise’
methodology

Description
Using an FPCA object create a functional box-plot based on the function scores. The green line
corresponds to the functional median, the dark grey area to the area spanned by the curves within
the 25th and 75-th percentile and the light gret to the area spanned by the curves within the 2.5th
and 97.5-th percentile.
Usage
CreateFuncBoxPlot(fpcaObj, optns = list(), ...)
Arguments
fpcaObj

An object of class FPCA returned by the function FPCA().

optns

A list of options control parameters specified by list(name=value). See ‘Details’.

...

Additional arguments for the ’plot’ function.

Details
Available control options are
ifactor inflation ifactor for the bag-plot defining the loop of bag-plot or multiplying ifactor the
KDE pilot bandwidth matrix. (see ?aplpack::compute.bagplot; ?ks::Hpi respectively; default:
2.58; 2 respectively).
variant string defining the method used (’KDE’, ’pointwise’ or ’bagplot’) (default: ’bagplot’)
unimodal logical specifying if the KDE estimate should be unimodal (default: FALSE, relavant
only for variant=’KDE’)
addIndx vector of indeces corresponding to which samples one should overlay (Default: NULL)
K integer number of the first K components used for the representation. (default: length(fpcaObj$lambda
))
References
P. J. Rousseeuw, I. Ruts, J. W. Tukey (1999): The bagplot: a bivariate boxplot, The American
Statistician, vol. 53, no. 4, 382-387

10

CreateModeOfVarPlot

Examples
set.seed(1)
n <- 20
pts <- seq(0, 1, by=0.05)
sampWiener <- Wiener(n, pts)
sampWiener <- Sparsify(sampWiener, pts, 10)
res <- FPCA(sampWiener$Ly, sampWiener$Lt,
list(dataType='Sparse', error=FALSE, kernel='epan', verbose=TRUE))
CreateFuncBoxPlot(res, list(addIndx=c(1:3)) )

CreateModeOfVarPlot

Functional Principal Component Analysis mode of variation plot

Description
Create the k-th mode of variation plot around the mean. The red-line
√ is the functional mean, the
grey shaded areas show the range of variations around the mean: ±Q λk φk for the dark grey area
Q = 1, and for the light grey are Q = 2. In the case of ’rainbow’ plot the blue edge corresponds to
Q = -3, the green edge to Q = +3 and the red-line to Q = 0 (the mean). In the case of ’minimal’ the
blue line corresponds to Q = -2 and green line to Q = 2.
Usage
CreateModeOfVarPlot(fpcaObj, k = 1, plotType = "standard",
colSpectrum = NULL, ...)
Arguments
fpcaObj

An FPCA class object returned by FPCA().

k

The k-th mode of variation to plot (default k = 1)

plotType

Character indicating the creation of "standard" shaded plot, a "rainbow"-plot or
a "minimal" sketch plot (default: "standard")

colSpectrum

Character vector to be use as input in the ’colorRampPalette’ function defining
the outliers colours (default: c("blue","red", "green"), relavant only for rainbowPlot=TRUE)

...

Additional arguments for the plot function.

Examples
set.seed(1)
n <- 20
pts <- seq(0, 1, by=0.05)
sampWiener <- Wiener(n, pts)
sampWiener <- Sparsify(sampWiener, pts, 10)
res <- FPCA(sampWiener$Ly, sampWiener$Lt,
list(dataType='Sparse', error=FALSE, kernel='epan', verbose=TRUE))
CreateModeOfVarPlot(res)

CreateOutliersPlot

11

CreateOutliersPlot

Functional Principal Component or Functional Singular Value Decomposition Scores Plot using ’bagplot’ or ’KDE’ methodology

Description
This function will create, using the first components scores, a set of convex hulls of the scores based
on ’bagplot’ or ’KDE’ methodology.
Usage
CreateOutliersPlot(fObj, optns = NULL, ...)
Arguments
fObj

A class object returned by FPCA() or FSVD().

optns

A list of options control parameters specified by list(name=value). See ‘Details’.

...

Additional arguments for the ’plot’ function.

Details
Available control options are
ifactor inflation ifactor for the bag-plot defining the loop of bag-plot or multiplying ifactor the
KDE pilot bandwidth matrix. (see ?aplpack::compute.bagplot; ?ks::Hpi respectively; default:
2.58; 2 respectively).
variant string defining the outlier method used (’KDE’, ’NN’ or ’bagplot’) (default: ’KDE’)
unimodal logical specifying if the KDE estimate should be unimodal (default: FALSE, relavant
only for variant=’KDE’)
maxVar logical specifying if during slicing we should used the directions of maximum variance
(default: FALSE for FPCA, TRUE for FSVD)
nSlices integer between 3 and 16, denoting the number of slices to be used (default: 4, relavant
only for groupingType=’slice’)
showSlices logical specifying if during slicing we should show the outline of the slice (default:
FALSE)
colSpectrum character vector to be use as input in the ’colorRampPalette’ function defining the
outliers colours (default: c("red", "yellow", ’blue’), relavant only for groupingType=’slice’)
groupingType string specifying if a slice grouping (’slice’) or a standard percentile/bagplot grouping (’standard’) should be returned (default: ’standard’)
fIndeces a two-component vector with the index of the mode of variation to consider (default:
c(1,2) for FPCA and c(1,1) for FSVD)
Value
An (temporarily) invisible copy of a list containing the labels associated with each of sample curves.

12

CreatePathPlot

References
P. J. Rousseeuw, I. Ruts, J. W. Tukey (1999): The bagplot: a bivariate boxplot, The American
Statistician, vol. 53, no. 4, 382-387 R. J. Hyndman and H. L. Shang. (2010) Rainbow plots,
bagplots, and boxplots for functional data, Journal of Computational and Graphical Statistics, 19(1),
29-45
Examples
## Not run:
set.seed(1)
n <- 420
pts <- seq(0, 1, by=0.05)
sampWiener <- Wiener(n, pts)
sampWiener <- Sparsify(sampWiener, pts, 10)
res <- FPCA(sampWiener$Ly, sampWiener$Lt,
list(dataType='Sparse', error=FALSE, kernel='epan', verbose=TRUE))
CreateOutliersPlot(res)
## End(Not run)

CreatePathPlot

Create the fitted sample path plot based on the results from FPCA().

Description
Create the fitted sample path plot based on the results from FPCA().
Usage
CreatePathPlot(fpcaObj, subset, K = NULL,
inputData = fpcaObj[["inputData"]], showObs = !is.null(inputData),
obsOnly = FALSE, showMean = FALSE, derOptns = list(p = 0), ...)
Arguments
fpcaObj

Returned object from FPCA().

subset

A vector of indices or a logical vector for subsetting the observations.

K

The number of components to reconstruct the fitted sample paths.

inputData

A list of length 2 containing the sparse/dense (unsupported yet) observations.
inputData needs to contain two fields: Lt for a list of time points and Ly for a
list of observations. Default to the ‘inputData‘ field within ‘fpcaObj‘.

showObs

Whether to plot the original observations for each subject.

obsOnly

Whether to show only the original curves.

showMean

Whether to plot the mean function as a bold solid curve.

derOptns

A list of options to control derivation parameters; see ‘fitted.FPCA’. (default =
NULL)

...

other arguments passed into matplot for plotting options

CreateScreePlot

13

Examples
set.seed(1)
n <- 20
pts <- seq(0, 1, by=0.05)
sampWiener <- Wiener(n, pts)
sampWiener <- Sparsify(sampWiener, pts, 10)
res <- FPCA(sampWiener$Ly, sampWiener$Lt,
list(dataType='Sparse', error=FALSE, kernel='epan',
verbose=TRUE))
CreatePathPlot(res, subset=1:5)
# CreatePathPlot has a lot of usages:
## Not run:
CreatePathPlot(res)
CreatePathPlot(res, 1:20)
CreatePathPlot(res, 1:20, showObs=FALSE)
CreatePathPlot(res, 1:20, showMean=TRUE, showObs=FALSE)
CreatePathPlot(res, 1:20, obsOnly=TRUE)
CreatePathPlot(res, 1:20, obsOnly=TRUE, showObs=FALSE)
CreatePathPlot(inputData=sampWiener, subset=1:20, obsOnly=TRUE)
## End(Not run)

CreateScreePlot

Create the scree plot for the fitted eigenvalues

Description
This function will open a new device if not instructed otherwise.
Usage
CreateScreePlot(fpcaObj, ...)
Arguments
fpcaObj

A object of class FPCA returned by the function FPCA().

...

Additional arguments for the ’plot’ function.

Examples
set.seed(1)
n <- 20
pts <- seq(0, 1, by=0.05)
sampWiener <- Wiener(n, pts)
sampWiener <- Sparsify(sampWiener, pts, 10)
res <- FPCA(sampWiener$Ly, sampWiener$Lt,
list(dataType='Sparse', error=FALSE, kernel='epan', verbose=TRUE))
CreateScreePlot(res)

14

cumtrapzRcpp

CreateStringingPlot

Create plots for observed and stringed high dimensional data

Description
The function produces the following three plots: 1) A plot of predictors (standardized if specified
so during stringing) in original order for a subset of observations; 2) A plot of predictors in stringed
order for the same subset of observations; 3) A plot of the stringing function, which is the stringed
order vs. the original order.
Usage
CreateStringingPlot(stringingObj, subset, ...)
Arguments
stringingObj

A stringing object of class "Stringing", returned by the function Stringing.

subset

A vector of indices or a logical vector for subsetting the observations. If missing,
first min(n,50) observations will be plotted where n is the sample size.

...

Other arguments passed into matplot for plotting options

Examples
set.seed(1)
n <- 50
wiener = Wiener(n = n)[,-1]
p = ncol(wiener)
rdmorder = sample(size = p, x=1:p, replace = FALSE)
stringingfit = Stringing(X = wiener[,rdmorder], disOptns = "correlation")
diff_norev = sum(abs(rdmorder[stringingfit$StringingOrder] - 1:p))
diff_rev = sum(abs(rdmorder[stringingfit$StringingOrder] - p:1))
if(diff_rev <= diff_norev){
stringingfit$StringingOrder = rev(stringingfit$StringingOrder)
stringingfit$Ly = lapply(stringingfit$Ly, rev)
}
CreateStringingPlot(stringingfit, 1:20)

cumtrapzRcpp

Cumulative Trapezoid Rule Numerical Integration

Description
Cumulative Trapezoid Rule Numerical Integration using Rcpp
Usage
cumtrapzRcpp(X, Y)

DynCorr

15

Arguments
X
Y

Sorted vector of X values
Vector of Y values.

DynCorr

Dynamical Correlation

Description
Calculate Dynamical Correlation for 2 paired dense regular functional data observed on the same
grid.
Usage
DynCorr(x, y, t)
Arguments
x
y
t

a n by m matrix where rows representing subjects and columns representing
measurements, missings are allowed.
a n by m matrix where rows representing subjects and columns representing
measurements, missings are allowed.
a length m vector of time points where x,y are observed.

Value
A length m vector of individual dynamic correlations
References
Dubin J A, M\"uller H G. Dynamical correlation for multivariate longitudinal data[J]. Journal of
the American Statistical Association, 2005, 100(471): 872-881. Liu S, Zhou Y, Palumbo R, et
al. Dynamical correlation: A new method for quantifying synchrony with multivariate intensive
longitudinal data[J]. Psychological methods, 2016, 21(3): 291.
Examples
set.seed(10)
n=200
# sample size
t=seq(0,1,length.out=100)
# length of data
mu_quad_x=8*t^2-4*t+5
mu_quad_y=8*t^2-12*t+6
fun=rbind(rep(1,length(t)),-t,t^2)
z1=matrix(0,n,3)
z1[,1]=rnorm(n,0,2)
z1[,2]=rnorm(n,0,16/3)
z1[,3]=rnorm(n,0,4)
x1_quad_error=y1_quad_error=matrix(0,nrow=n,ncol=length(t))
for (i in 1:n){
x1_quad_error[i,]=mu_quad_x+z1[i,]%*%fun+rnorm(length(t),0,0.01)
y1_quad_error[i,]=mu_quad_y+2*z1[i,]%*%fun +rnorm(length(t),0,0.01)
}
dyn1_quad=DynCorr(x1_quad_error,y1_quad_error,t)

16

Dyn_test

Dyn_test

Bootstrap test of Dynamic correlation

Description
Perform one sample (H0: Dynamic correlation = 0) or two sample (H0:Dynamic_correlation_1 =
Dynamic_correlation_2) bootstrap test of Dynamical Correlation.
Usage
Dyn_test(x1, y1, t1, x2, y2, t2, B = 1000)
Arguments
x1

a n by m matrix where rows representing subjects and columns representing
measurements, missings are allowed.

y1

a n by m matrix where rows representing subjects and columns representing
measurements, missings are allowed.

t1

a vector of time points where x1,y1 are observed.

x2

(optional if missing will be one sample test) a n by m matrix where rows representing subjects and columns representing measurements, missings are allowed.

y2

(optional if missing will be one sample test) a n by m matrix where rows representing subjects and columns representing measurements, missings are allowed.

t2

(optional if missing will be one sample test) a vector of time points where x2,y2
are observed.

B

number of bootstrap samples.

Value
a list of the following
stats: test statistics.
pval: p-value of the test.

References
Dubin J A, M\"uller H G. Dynamical correlation for multivariate longitudinal data[J]. Journal of the
American Statistical Association, 2005, 100(471): 872-881.
Liu S, Zhou Y, Palumbo R, et al. Dynamical correlation: A new method for quantifying synchrony
with multivariate intensive longitudinal data[J]. Psychological methods, 2016, 21(3): 291.
Examples
n=200
# sample size
t=seq(0,1,length.out=100)
# length of data
mu_quad_x=8*t^2-4*t+5
mu_quad_y=8*t^2-12*t+6
fun=rbind(rep(1,length(t)),-t,t^2)

FAM

17

z1=matrix(0,n,3)
z1[,1]=rnorm(n,0,2)
z1[,2]=rnorm(n,0,16/3)
z1[,3]=rnorm(n,0,4)
# covariance matrix of random effects
x1_quad_error=y1_quad_error=matrix(0,nrow=n,ncol=length(t))
for (i in 1:n){
x1_quad_error[i,]=mu_quad_x+z1[i,]%*%fun+rnorm(length(t),0,0.01)
y1_quad_error[i,]=mu_quad_y+2*z1[i,]%*%fun +rnorm(length(t),0,0.01)
}
bt_DC=Dyn_test(x1_quad_error,y1_quad_error,t,B=1000)

FAM

Functional Additive Models

Description
Functional additive models with a single predictor process
Usage
FAM(Y, Lx, Lt, nEval = 51, newLx = NULL, newLt = NULL,
bwMethod = 0, alpha = 0.7, supp = c(-2, 2), optns = NULL)
Arguments
Y
Lx
Lt
nEval

newLx
newLt
bwMethod
alpha
supp

optns

An n-dimensional vector whose elements consist of scalar responses.
A list of n vectors containing the observed values for each individual. See FPCA
for detail.
A list of n vectors containing the observation time points for each individual.
Each vector should be sorted in ascending order. See FPCA for detail.
The number of evaluation grid points for kernel smoothing (default is 51. If
it is specified as 0, then estimated FPC scores in the training set are used for
evaluation grid instead of equal grid).
A list of the observed values for test set. See predict.FPCA for detail.
A list of the observed time points for test set. See predict.FPCA for detail.
The method of bandwidth selection for kernel smoothing, a positive value for
designating K-fold cross-validtaion and zero for GCV (default is 50)
The shrinkage factor (positive number) for bandwidth selection. See Han et al.
(2016) (default is 0.7).
The lower and upper limits of kernel smoothing domain for studentized FPC
scores, which FPC scores are divided by the square roots of eigenvalues (default
is [-2,2]).
A list of options control parameters specified by list(name=value). See FPCA.

Details
FAM fits functional additive models for a scalar response and single predictor process proposed by
Mueller and Yao (2007) that
K
X
E(Y |X) =
gk (ξk ),
k=1

where ξk stand for the k-th FPC score of the the predictor process.

18

FAM

Value
A list containing the following fields:
mu

Mean estimator of EY

fam

A N by K matrix whose column vectors consist of the component function estimators at the given estimation points.

xi

An N by K matrix whose column vectors consist of N vectors of estimation
points for each component function.

bw

A K-dimensional bandwidth vector.

lambda

A K-dimensional vector containing eigenvalues.

phi

An nWorkGrid by K matrix containing eigenfunctions, supported by WorkGrid.
See FPCA.

workGrid

An nWorkGrid by K_j working grid, the internal regular grid on which the eigen
analysis is carried on. See FPCA.

References
Mueller, H.-G. and Yao, F. (2005), "Functional additive models", JASA, Vol.103, No.484, p.15341544.
Examples
set.seed(1000)
library(MASS)
f1
f2
f3
f4

<<<<-

function(t)
function(t)
function(t)
function(t)

0.5*t
2*cos(2*pi*t/4)
1.5*sin(2*pi*t/4)
2*atan(2*pi*t/4)

n<-250
N<-500
sig <- diag(c(4.0,2.0,1.5,1.2))
scoreX <- mvrnorm(n,mu=rep(0,4),Sigma=sig)
scoreXTest <- mvrnorm(N,mu=rep(0,4),Sigma=sig)
Y <- f1(scoreX[,1]) + f2(scoreX[,2]) + f3(scoreX[,3]) + f4(scoreX[,4]) + rnorm(n,0,0.1)
YTest <- f1(scoreXTest[,1]) + f2(scoreXTest[,2]) +
f3(scoreXTest[,3]) + f4(scoreXTest[,4]) + rnorm(N,0,0.1)
phi1
phi2
phi3
phi4

<<<<-

function(t)
function(t)
function(t)
function(t)

sqrt(2)*sin(2*pi*t)
sqrt(2)*sin(4*pi*t)
sqrt(2)*cos(2*pi*t)
sqrt(2)*cos(4*pi*t)

grid <- seq(0,1,length.out=21)
Lt <- Lx <- list()
for (i in 1:n) {
Lt[[i]] <- grid
Lx[[i]] <- scoreX[i,1]*phi1(grid) + scoreX[i,2]*phi2(grid) +
scoreX[i,3]*phi3(grid) + scoreX[i,4]*phi4(grid) + rnorm(1,0,0.01)

FAM
}
LtTest <- LxTest <- list()
for (i in 1:N) {
LtTest[[i]] <- grid
LxTest[[i]] <- scoreXTest[i,1]*phi1(grid) + scoreXTest[i,2]*phi2(grid) +
scoreXTest[i,3]*phi3(grid) + scoreXTest[i,4]*phi4(grid) + rnorm(1,0,0.01)
}
# estimation
fit <- FAM(Y=Y,Lx=Lx,Lt=Lt)
xi <- fit$xi
par(mfrow=c(2,2))
j <- 1
g1 <- f1(sort(xi[,j]))
tmpSgn <- sign(sum(g1*fit$fam[,j]))
plot(sort(xi[,j]),g1,type='l',col=2,ylim=c(-2.5,2.5),xlab='xi1')
points(sort(xi[,j]),tmpSgn*fit$fam[order(xi[,j]),j],type='l')
j <- 2
g2 <- f2(sort(xi[,j]))
tmpSgn <- sign(sum(g2*fit$fam[,j]))
plot(sort(xi[,j]),g2,type='l',col=2,ylim=c(-2.5,2.5),xlab='xi2')
points(sort(xi[,j]),tmpSgn*fit$fam[order(xi[,j]),j],type='l')
j <- 3
g3 <- f3(sort(xi[,j]))
tmpSgn <- sign(sum(g3*fit$fam[,j]))
plot(sort(xi[,j]),g3,type='l',col=2,ylim=c(-2.5,2.5),xlab='xi3')
points(sort(xi[,j]),tmpSgn*fit$fam[order(xi[,j]),j],type='l')
j <- 4
g4 <- f4(sort(xi[,j]))
tmpSgn <- sign(sum(g4*fit$fam[,j]))
plot(sort(xi[,j]),g4,type='l',col=2,ylim=c(-2.5,2.5),xlab='xi4')
points(sort(xi[,j]),tmpSgn*fit$fam[order(xi[,j]),j],type='l')
# fitting
fit <- FAM(Y=Y,Lx=Lx,Lt=Lt,nEval=0)
yHat <- fit$mu+apply(fit$fam,1,'sum')
par(mfrow=c(1,1))
plot(yHat,Y)
abline(coef=c(0,1),col=2)
# R^2
R2 <- 1-sum((Y-yHat)^2)/sum((Y-mean(Y))^2)
R2
# prediction
fit <- FAM(Y=Y,Lx=Lx,Lt=Lt,newLx=LxTest,newLt=LtTest)
yHat <- fit$mu+apply(fit$fam,1,'sum')

19

20

FCCor
par(mfrow=c(1,1))
plot(yHat,YTest,xlim=c(-10,10))
abline(coef=c(0,1),col=2)

FCCor

Calculate functional correlation between two simultaneously observed
processes.

Description
Calculate functional correlation between two simultaneously observed processes.

Usage
FCCor(x, y, Lt, bw = stop("bw missing"), kern = "epan",
Tout = sort(unique(unlist(Lt))))
Arguments
x

A list of function values corresponding to the first process.

y

A list of function values corresponding to the second process.

Lt

A list of time points for both x and y.

bw

A numeric vector for bandwidth of length either 5 or 1, specifying the bandwidths for E(X), E(Y), var(X), var(Y), and cov(X, Y). If bw is a scalar then all
five bandwidths are chosen to be the same.

kern

Smoothing kernel for mu and covariance; "rect", "gauss", "epan", "gausvar",
"quar" (default: "gauss")

Tout

Output time points. Default to the sorted unique time points.

Details
FCCor calculate only the concurrent correlation corr(X(t), Y(t)) (note that the time points t are the
same). It assumes no measurement error in the observed values.

Value
A list with the following components:
corr

A vector of the correlation corr(X(t), Y(t)) evaluated at Tout.

Tout

Same as the input Tout.

bw

The bandwidths used for E(X), E(Y), var(X), var(Y), and cov(X, Y).

FClust

21

Examples
set.seed(1)
n <- 200
nGridIn <- 50
sparsity <- 1:5 # must have length > 1
bw <- 0.2
kern <- 'epan'
T <- matrix(seq(0.5, 1, length.out=nGridIn))
## Corr(X(t), Y(t)) = 1/2
A <- Wiener(n, T)
B <- Wiener(n, T)
C <- Wiener(n, T) + matrix((1:nGridIn) , n, nGridIn, byrow=TRUE)
X <- A + B
Y <- A + C
indEach <- lapply(1:n, function(x) sort(sample(nGridIn, sample(sparsity, 1))))
tAll <- lapply(1:n, function(i) T[indEach[[i]]])
Xsp <- lapply(1:n, function(i) X[i, indEach[[i]]])
Ysp <- lapply(1:n, function(i) Y[i, indEach[[i]]])
plot(T, FCCor(Xsp, Ysp, tAll, bw)[['corr']], ylim=c(-1, 1))
abline(h=0.5)

FClust

Functional clustering and identifying substructures of longitudinal
data

Description
By default the function will cluster the data using the functional principal component (FPC) scores
from the data’s FPC analysis using EMCluster (Chen and Maitra, 2015) or directly clustering the
functional data using kCFC (Chiou and Li, 2007).
Usage
FClust(Ly, Lt, k = 3, cmethod = "EMCluster", optnsFPCA = NULL,
optnsCS = NULL)
Arguments
Ly

A list of n vectors containing the observed values for each individual. Missing
values specified by NAs are supported for dense case (dataType='dense').

Lt

A list of n vectors containing the observation time points for each individual
corresponding to y.

k

A scalar defining the number of clusters to define; default 3.

cmethod

A string specifying the clusterig method to use (’EMCluster’ or ’kCFC’); default: ’EMCluster’.

optnsFPCA

A list of options control parameters specified by list(name=value) to be used
for by FPCA on the sample y; by default: "list( methodMuCovEst =’smooth’,
FVEthreshold= 0.90, methodBwCov = ’GCV’, methodBwMu = ’GCV’ )". See
‘Details in ?FPCA’.

22

FClust
optnsCS

A list of options control parameters specified by list(name=value) to be used
for cluster-specific FPCA from kCFC; by default: "list( methodMuCovEst =’smooth’,
FVEthreshold= 0.70, methodBwCov = ’GCV’, methodBwMu = ’GCV’ )". See
‘Details in ?FPCA’ and ’?kCFC’.

Details
Within EMCluster we examine the model initiated "EMCluster::em.EM" and return the optimal
model based on ’EMCluster::emcluster’. See ?EMCluster::emcluster for details.

Value
A list containing the following fields:
cluster

A vector of levels 1:k, indicating the cluster to which each curve is allocated.

fpca

An FPCA object derived from the sample used by Rmixmod, otherwise NULL.

clusterObj

Either a EMCluster object or kCFC object.

References
Wei-Chen Chen and Ranjan Maitra, "EMCluster: EM Algorithm for Model-Based Clusttering of
Finite Mixture Gaussian Distribution". (2015)
Julien Jacques and Cristian Preda, "Funclust: A curves clustering method using functional random
variables density approximation". Neurocomputing 112 (2013): 164-171
Jeng-Min Chiou and Pai-Ling Li, "Functional clustering and identifying substructures of longitudinal data". Journal of the Royal Statistical Society B 69 (2007): 679-699

Examples
## Not run:
data(medfly25)
Flies <- MakeFPCAInputs(medfly25$ID, medfly25$Days, medfly25$nEggs)
newClust <- FClust(Flies$Ly, Flies$Lt, k = 2, optnsFPCA =
list(methodMuCovEst = 'smooth', userBwCov = 2, FVEthreshold = 0.90))
# We denote as 'veryLowCount' the group of flies that lay less
# than twenty-five eggs during the 25-day period examined.
veryLowCount = ifelse( sapply( unique(medfly25$ID), function(u)
sum( medfly25$nEggs[medfly25$ID == u] )) < 25, 0, 1)
N <- length(unique(medfly25$ID))
(correctRate <- sum( (1 + veryLowCount) == newClust$cluster) / N) # 99.6%
## End(Not run)

FCReg

FCReg

23

Functional Concurrent Regression by 2D smoothing method.

Description
Functional concurrent regression with dense or sparse functional data for scalar or functional dependent variable.
Usage
FCReg(vars, userBwMu, userBwCov, outGrid, kern = "gauss",
measurementError = TRUE, diag1D = "none", useGAM = FALSE,
returnCov = TRUE)
Arguments
vars

A list of input functional/scalar covariates. Each field corresponds to a functional (a list) or scalar (a vector) covariate. The last entry is assumed to be the
response if no entry is names ’Y’. If a field corresponds to a functional covariate,
it should have two fields: ’Lt’, a list of time points, and ’Ly’, a list of function
values.

userBwMu

A scalar with bandwidth used for smoothing the mean

userBwCov

A scalar with bandwidth used for smoothing the auto- and cross-covariances

outGrid

A vector with the output time points

kern

Smoothing kernel choice, common for mu and covariance; "rect", "gauss", "epan",
"gausvar", "quar" (default: "gauss")

measurementError
Indicator measurement errors on the functional observations should be assumed.
If TRUE the diagonal raw covariance will be removed when smoothing. (default: TRUE)
diag1D

A string specifying whether to use 1D smoothing for the diagonal line of the
covariance. ’none’: don’t use 1D smoothing; ’cross’: use 1D only for crosscovariances; ’all’: use 1D for both auto- and cross-covariances. (default : ’none’)

useGAM

Indicator to use gam smoothing instead of local-linear smoothing (semi-parametric
option) (default: FALSE)

returnCov

Indicator to return the covariance surfaces, which is a four dimensional array.
The first two dimensions correspond to outGrid and the last two correspond to
the covariates and the response, i.e. (i, j, k, l) entry being Cov(X_k(t_i), X_l(t_j))
(default: FALSE)

...

Additional arguments

Details
If measurement error is assumed, the diagonal elements of the raw covariance will be removed. This
could result in highly unstable estimate if the design is very sparse, or strong seasonality presents.

24

fdapace

References
Yao, F., Mueller, H.G., Wang, J.L. "Functional Linear Regression Analysis for Longitudinal Data."
Annals of Statistics 33, (2005): 2873-2903.(Dense data)
Senturk, D., Nguyen, D.V. "Varying Coefficient Models for Sparse Noise-contaminated Longitudinal Data", Statistica Sinica 21(4), (2011): 1831-1856. (Sparse data)
Examples
# Y(t) = \beta_0(t) + \beta_1(t) X_1(t) + \beta_2(t) Z_2 + \epsilon
# Settings
set.seed(1)
n <- 75
nGridIn <- 150
sparsity <- 5:10 # Sparse data sparsity
T <- round(seq(0, 1, length.out=nGridIn), 4) # Functional data support
bw <- 0.1
outGrid <- round(seq(min(T), 1, by=0.05), 2)
# Simulate functional data
mu <- T * 2 # mean function for X_1
sigma <- 1
beta_0 <- 0
beta_1 <- 1
beta_2 <- 1
Z <- MASS::mvrnorm(n, rep(0, 2), diag(2))
X_1 <- Z[, 1, drop=FALSE] %*% matrix(1, 1, nGridIn) + matrix(mu, n, nGridIn, byrow=TRUE)
epsilon <- rnorm(n, sd=sigma)
Y <- matrix(NA, n, nGridIn)
for (i in seq_len(n)) {
Y[i, ] <- beta_0 + beta_1 * X_1[i, ] + beta_2 * Z[i, 2] + epsilon[i]
}
# Sparsify functional data
set.seed(1)
X_1sp <- Sparsify(X_1, T, sparsity)
set.seed(1)
Ysp <- Sparsify(Y, T, sparsity)
vars <- list(X_1=X_1sp, Z_2=Z[, 2], Y=Ysp)
withError2D <- FCReg(vars, bw, bw, outGrid)

fdapace

PACE: Principal Analysis by Conditional Expectation

Description
PACE package for Functional Data Analysis and Empirical Dynamics.

fitted.FPCA

25

Details
PACE is a versatile package that provides implementation of various methods of Functional Data
Analysis (FDA) and Empirical Dynamics. The core of this package is Functional Principal Component Analysis (FPCA), a key technique for functional data analysis, for sparsely or densely sampled
random trajectories and time courses, via the Principal Analysis by Conditional Estimation (PACE)
algorithm. PACE is useful for the analysis of data that have been generated by a sample of underlying (but usually not fully observed) random trajectories. It does not rely on pre-smoothing of
trajectories, which is problematic if functional data are sparsely sampled. PACE provides options
for functional regression and correlation, for Longitudinal Data Analysis, the analysis of stochastic
processes from samples of realized trajectories, and for the analysis of underlying dynamics.
Author(s)

Xiongtao Dai , Pantelis Z. Hadjipantelis , Hao Ji
 Hans-Georg Mueller  Jane-Ling Wang 
Maintainer: Pantelis Z. Hadjipantelis 

fitted.FPCA

Fitted functional sample from FPCA object

Description
Combine the zero-meaned fitted values and the interpolated mean to get the fitted values for the trajectories or the derivatives of these trajectories. Estimates are given on the work-grid, not on the observation grid. Use ConvertSupport to map the estimates to your desired domain. 100*(1-alpha)percentage coverage intervals, or bands, for trajectory estimates (not derivatives) are provided. See
details in example.
Usage
## S3 method for class 'FPCA'
fitted(object, K = NULL, derOptns = list(p = 0),
ciOptns = list(alpha = NULL, cvgMethod = NULL), ...)
Arguments
object

A object of class FPCA returned by the function FPCA().

K

The integer number of the first K components used for the representation. (default: length(fpcaObj$lambda ))

derOptns

A list of options to control the derivation parameters specified by list(name=value).
See ‘Details’. (default = NULL)

ciOptns

A list of options to control the confidence interval/band specified by list(name=value).
See ‘Details’. (default = NULL)

...

Additional arguments

26

fitted.FPCA

Details
Available derivation control options are
p The order of the derivatives returned (default: 0, max: 2)
method The method used to produce the sample of derivatives (’FPC’ (default) or ’QUO’). See
Liu and Mueller (2009) for more details
bw Bandwidth for smoothing the derivatives (default: p * 0.10 * S)
kernelType Smoothing kernel choice; same available types are FPCA(). default(’epan’)
Available confidence interval/band control options are
alpha Significant level for confidence interval/band for trajectory coverage. default=0.05 (currently only work when p=0)
cvgMethod Option for trajectory coverage method beween ’interval’ and ’band’. default=’band’
Value
If alpha is NULL, p>1 or functional observations are dense, an n by length(workGrid) matrix, each
row of which contains a sample. Otherwise, it returns a list which consists of the following items:
workGrid

An evaluation grid for fitted values.

fitted

An n by length(workGrid) matrix, each row of which contains a sample.

cvgUpper

An n by length(workGrid) matrix, each row of which contains the upper alphacoverage limit

cvgLower

An n by length(workGrid) matrix, each row of which contains the lower alphacoverage limit

References
Yao, F., Mueller, H.-G. and Wang, J.-L. "Functional data analysis for sparse longitudinal data",
Journal of the American Statistical Association, vol.100, No. 470 (2005): 577-590.
Liu, Bitao, and Hans-Georg Mueller. "Estimating derivatives for samples of sparsely observed functions, with application to online auction dynamics." Journal of the American Statistical Association
104, no. 486 (2009): 704-717. (Sparse data FPCA)
Examples
set.seed(1)
n <- 100
pts <- seq(0, 1, by=0.05)
sampWiener <- Wiener(n, pts)
sampWiener <- Sparsify(sampWiener, pts, 5:10)
res <- FPCA(sampWiener$Ly, sampWiener$Lt,
list(dataType='Sparse', error=FALSE, kernel='epan', verbose=TRUE))
fittedY <- fitted(res, ciOptns = list(alpha=0.05))
workGrid <- res$workGrid
cvgUpper <- fittedY$cvgUpper
cvgLower <- fittedY$cvgLower
par(mfrow=c(2,3))
ind <- sample(1:n,6)
for (i in 1:6) {

fitted.FPCAder

27

j <- ind[i]
plot(workGrid,cvgUpper[j,],type='l',ylim=c(min(cvgLower[j,]),max(cvgUpper[j,])),col=4,lty=2,
xlab='t', ylab='X(t)', main=paste(j,'-th subject',sep=''))
points(workGrid,cvgLower[j,],type='l',col=4,lty=2)
points(res$inputData$Lt[[j]],res$inputData$Ly[[j]])
}

fitted.FPCAder

Fitted functional sample from FPCAder object

Description
Combine the zero-meaned fitted values and the mean derivative to get the fitted values for the
derivatives trajectories Estimates are given on the work-grid, not on the observation grid. Use
ConvertSupport to map the estimates to your desired domain.
Usage
## S3 method for class 'FPCAder'
fitted(object, K = NULL, ...)
Arguments
object

A object of class FPCA returned by the function FPCA().

K

The integer number of the first K components used for the representation. (default: length(derObj$lambda ))

...

Additional arguments

Value
An n by length(workGrid) matrix, each row of which contains a sample.
References
Liu, Bitao, and Hans-Georg Mueller. "Estimating derivatives for samples of sparsely observed functions, with application to online auction dynamics." Journal of the American Statistical Association
104, no. 486 (2009): 704-717. (Sparse data FPCA)
Examples
set.seed(1)
n <- 20
pts <- seq(0, 1, by=0.05)
sampWiener <- Wiener(n, pts)
sampWiener <- Sparsify(sampWiener, pts, 10)

28

FOptDes

FOptDes

Optimal Designs for Functional and Longitudinal Data for Trajectory
Recovery or Scalar Response Prediction

Description
Optimal Designs for Functional and Longitudinal Data for Trajectory Recovery or Scalar Response
Prediction
Usage
FOptDes(Ly, Lt, Resp, p = 3, optns = list(),
isRegression = !missing(Resp), isSequential = FALSE,
RidgeCand = NULL)
Arguments
Ly

A list of n vectors containing the observed values for each individual. Missing
values specified by NAs are supported for dense case (dataType='dense').

Lt

A list of n vectors containing the observation time points for each individual
corresponding to y. Each vector should be sorted in ascending order.

Resp

A vector of response values, keep void for trajectory recovery, only necessary
for scalar response prediction task.

p

A fixed positive integer indicating the number of optimal design points requested,
with default: 3.

optns

A list of options control parameters specified by list(name=value) for FPCA,
with default: list().

isRegression

A logical argument, indicating the purpose of the optimal designs: TRUE for
scalar response prediction, FALSE for trajectory recovery, with default value
!missing(Resp).

isSequential

A logical argument, indicating whether to use the sequential optimization procedure for faster computation, recommended for relatively large p (default: FALSE).

RidgeCand

A vector of positive numbers as ridge penalty candidates for regularization. The
final value is selected via cross validation. If only 1 ridge parameter is specified,
CV procedure is skipped.

Details
To select a proper RidgeCand, check with the returned optimal ridge parameter. If the selected
parameter is the maximum/minimum values in the candidates, it is possible that the selected one is
too small/big.
Value
A list containing the following fields:
OptDes

The vector of optimal design points of the regular time grid of the observed data.

R2

Coefficient of determination. (Check the paper for details.)

R2adj

Adjusted coefficient of determination.

OptRidge

The selected ridge parameter.

FPCA

29

References
Ji, H., Mueller, H.G. (2016) "Optimal Designs for Longitudinal and Functional Data" Journal of the
Royal Statistical Society: Series B (Statistical Methodology)
Examples
set.seed(1)
n <- 50
pts <- seq(0, 1, by=0.05)
sampWiener <- Wiener(n, pts)
sampWiener <- MakeFPCAInputs(IDs = rep(1:n, each=length(pts)),
tVec = rep(pts, times = n),
yVec = t(sampWiener))
res <- FOptDes(Ly=sampWiener$Ly, Lt=sampWiener$Lt, p=2,
isSequential=FALSE, RidgeCand = seq(0.02,0.2,0.02))

FPCA

Functional Principal Component Analysis

Description
FPCA for dense or sparse functional data.
Usage
FPCA(Ly, Lt, optns = list())
Arguments
Ly

A list of n vectors containing the observed values for each individual. Missing
values specified by NAs are supported for dense case (dataType='dense').

Lt

A list of n vectors containing the observation time points for each individual
corresponding to y. Each vector should be sorted in ascending order.

optns

A list of options control parameters specified by list(name=value). See ‘Details’.

Details
If the input is sparse data, make sure you check the design plot is dense and the 2D domain is well
covered, using plot or CreateDesignPlot. Some study design such as snippet data (each subject
is observed only on a sub-interval of the period of study) will have an ill-covered design plot, for
which the covariance estimate will be unreliable.
Available control options are
userBwCov The bandwidth value for the smoothed covariance function; positive numeric - default:
determine automatically based on ’methodBwCov’
methodBwCov The bandwidth choice method for the smoothed covariance function; ’GMeanAndGCV’
(the geometric mean of the GCV bandwidth and the minimum bandwidth),’CV’,’GCV’ - default: 10% of the support
userBwMu The bandwidth value for the smoothed mean function (using ’CV’ or ’GCV’); positive
numeric - default: determine automatically based on ’methodBwMu’

30

FPCA
methodBwMu The bandwidth choice method for the mean function; ’GMeanAndGCV’ (the geometric mean of the GCV bandwidth and the minimum bandwidth),’CV’,’GCV’ - default: 5%
of the support
dataType The type of design we have (usually distinguishing between sparse or dense functional
data); ’Sparse’, ’Dense’, ’DenseWithMV’, ’p»n’ - default: determine automatically based on
’IsRegular’
diagnosticsPlot Deprecated. Same as the option ’plot’
plot Plot FPCA results (design plot, mean, scree plot and first K (<=3) eigenfunctions); logical default: FALSE
error Assume measurement error in the dataset; logical - default: TRUE
fitEigenValues Whether also to obtain a regression fit of the eigenvalues - default: FALSE
FVEthreshold Fraction-of-Variance-Explained threshold used during the SVD of the fitted covar.
function; numeric (0,1] - default: 0.9999
kernel Smoothing kernel choice, common for mu and covariance; "rect", "gauss", "epan", "gausvar", "quar" - default: "gauss"; dense data are assumed noise-less so no smoothing is performed.
kFoldMuCov The number of folds to be used for mean and covariance smoothing. Default: 10
lean If TRUE the ’inputData’ field in the output list is empty. Default: FALSE
maxK The maximum number of principal components to consider - default: min(20, N-1), N:# of
curves
methodXi The method to estimate the PC scores; ’CE’ (Condit. Expectation), ’IN’ (Numerical
Integration) - default: ’CE’ for sparse data and dense data with missing values, ’IN’ for dense
data.
methodMuCovEst The method to estimate the mean and covariance in the case of dense functional
data; ’cross-sectional’, ’smooth’ - default: ’cross-sectional’
nRegGrid The number of support points in each direction of covariance surface; numeric - default:
51
numBins The number of bins to bin the data into; positive integer > 10, default: NULL
methodSelectK The method of choosing the number of principal components K; ’FVE’,’AIC’,’BIC’,
or a positive integer as specified number of components: default ’FVE’)
shrink Whether to use shrinkage method to estimate the scores in the dense case (see Yao et al
2003) - default FALSE
outPercent A 2-element vector in [0,1] indicating the outPercent data in the boundary - default
(0,1)
rho The truncation threshold for the iterative residual. ’cv’: choose rho by leave-one-observation
out cross-validation; ’no’: no regularization - default "cv" if error == TRUE, and "no" if error
== FALSE.
rotationCut The 2-element vector in [0,1] indicating the percent of data truncated during sigma^2
estimation; default (0.25, 0.75))
useBinnedData Should the data be binned? ’FORCE’ (Enforce the # of bins), ’AUTO’ (Select the
# of bins automatically), ’OFF’ (Do not bin) - default: ’AUTO’
useBinnedCov Whether to use the binned raw covariance for smoothing; logical - default:TRUE
userCov The user-defined smoothed covariance function; list of two elements: numerical vector
’t’ and matrix ’cov’, ’t’ must cover the support defined by ’Ly’ - default: NULL
userMu The user-defined smoothed mean function; list of two numerical vector ’t’ and ’mu’ of
equal size, ’t’ must cover the support defined ’Ly’ - default: NULL

FPCA

31

userSigma2 The user-defined measurement error variance. A positive scalar. If specified then no
regularization is used (rho is set to ’no’, unless specified otherwise). Default to ‘NULL‘
userRho The user-defined measurement truncation threshold used for the calculation of functional
principal components scores. A positive scalar. Default to ‘NULL‘
useBW1SE Pick the largest bandwidth such that CV-error is within one Standard Error from the
minimum CV-error, relevant only if methodBwMu =’CV’ and/or methodBwCov =’CV’; logical - default: FALSE
verbose Display diagnostic messages; logical - default: FALSE

Value
A list containing the following fields:
sigma2

Variance for measure error.

lambda

A vector of length K containing eigenvalues.

phi

An nWorkGrid by K matrix containing eigenfunctions, supported on workGrid.

xiEst

A n by K matrix containing the FPC estimates.

xiVar

A list of length n, each entry containing the variance estimates for the FPC
estimates.

obsGrid

The (sorted) grid points where all observation points are pooled.

mu

A vector of length nWorkGrid containing the mean function estimate.

workGrid

A vector of length nWorkGrid. The internal regular grid on which the eigen
analysis is carried on.

smoothedCov

A nWorkGrid by nWorkGrid matrix of the smoothed covariance surface.

fittedCov

A nWorkGrid by nWorkGrid matrix of the fitted covariance surface, which is
guaranteed to be non-negative definite.

optns

A list of actually used options.

timings

A vector with execution times for the basic parts of the FPCA call.

bwMu

The selected (or user specified) bandwidth for smoothing the mean function.

bwCov

The selected (or user specified) bandwidth for smoothing the covariance function.

rho

A regularizing scalar for the measurement error variance estimate.

cumFVE

A vector with the percentages of the total variance explained by each FPC. Increase to almost 1.

FVE

A percentage indicating the total variance explained by chosen FPCs with corresponding ’FVEthreshold’.

criterionValue A scalar specifying the criterion value obtained by the selected number of components with specific methodSelectK: FVE,AIC,BIC values or NULL for fixedK.
inputData

A list containting the original ’Ly’ and ’Lt’ lists used as inputs to FPCA. NULL
if ’lean’ was specified to be TRUE.

32

FPCAder

References
Yao, F., Mueller, H.G., Clifford, A.J., Dueker, S.R., Follett, J., Lin, Y., Buchholz, B., Vogel, J.S.
(2003). "Shrinkage estimation for functional principal component scores, with application to the
population kinetics of plasma folate." Biometrics 59, 676-685. (Shrinkage estimates for dense data)
Yao, Fang, Hans-Georg Mueller, and Jane-Ling Wang. "Functional data analysis for sparse longitudinal data." Journal of the American Statistical Association 100, no. 470 (2005): 577-590. (Sparse
data FPCA)
Liu, Bitao, and Hans-Georg Mueller. "Estimating derivatives for samples of sparsely observed functions, with application to online auction dynamics." Journal of the American Statistical Association
104, no. 486 (2009): 704-717. (Sparse data FPCA)
Castro, P. E., W. H. Lawton, and E. A. Sylvestre. "Principal modes of variation for processes with
continuous sample curves." Technometrics 28, no. 4 (1986): 329-337. (Dense data FPCA)
Examples
set.seed(1)
n <- 20
pts <- seq(0, 1, by=0.05)
sampWiener <- Wiener(n, pts)
sampWiener <- Sparsify(sampWiener, pts, 10)
res <- FPCA(sampWiener$Ly, sampWiener$Lt,
list(dataType='Sparse', error=FALSE, kernel='epan', verbose=TRUE))
plot(res) # The design plot covers [0, 1] * [0, 1] well.
CreateCovPlot(res, 'Fitted')

FPCAder

Take derivative of an FPCA object

Description
Take derivative of an FPCA object
Usage
FPCAder(fpcaObj, derOptns = list(p = 1))
Arguments
fpcaObj

A object of class FPCA returned by the function FPCA().

derOptns

A list of options to control the derivation parameters specified by list(name=value).
See ‘Details’. (default = NULL)

Details
Available derivation control options are
method The method used for obtaining the derivatives (default: ’FPC’). ’DPC’: derivatives principal component, with G^(1,1) estimated by first kernel local smoothing G^(1,0), and then
apply a 1D smoother on the second direction; ’FPC’: functional principal component, based
on smoothing the eigenfunctions; ’FPC1’: functional principal component, based on smoothing G^(1,0). May produce better estimate than ’FPC’ but is slower.

FPCAder

33

p The order of the derivatives returned (default: 1, max: 2).
bw Bandwidth for the 1D and the 2D smoothers (default: p * 0.1 * S).
kernelType Smoothing kernel choice; same available types are FPCA(). default(’epan’)

References
Dai, Xiongtao, Hans-Georg Mueller, and Wenwen Tao. "Derivative Principal Component Analysis
for Representing the Time Dynamics of Longitudinal and Functional Data." Submitted (DPC) Liu,
Bitao, and Hans-Georg Mueller. "Estimating derivatives for samples of sparsely observed functions,
with application to online auction dynamics." Journal of the American Statistical Association 104,
no. 486 (2009): 704-717. (FPC)

Examples
bw <- 0.2
kern <- 'epan'
set.seed(1)
n <- 100
M <- 40
pts <- seq(0, 1, length.out=M)
lambdaTrue <- c(1, 0.8, 0.1)^2
sigma2 <- 0.1
samp2 <- MakeGPFunctionalData(n, M, pts, K=length(lambdaTrue),
lambda=lambdaTrue, sigma=sqrt(sigma2), basisType='legendre01')
samp2 <- c(samp2, MakeFPCAInputs(tVec=pts, yVec=samp2$Yn))
fpcaObj <- FPCA(samp2$Ly, samp2$Lt, list(methodMuCovEst='smooth',
userBwCov=bw, userBwMu=bw, kernel=kern, error=TRUE))
CreatePathPlot(fpcaObj, showObs=FALSE)
FPCoptn <- list(bw=bw, kernelType=kern, method='FPC')
DPCoptn <- list(bw=bw, kernelType=kern, method='DPC')
FPC <- FPCAder(fpcaObj, FPCoptn)
DPC <- FPCAder(fpcaObj, DPCoptn)
CreatePathPlot(FPC, ylim=c(-5, 10))
CreatePathPlot(DPC, ylim=c(-5, 10))
# Get the true derivatives
phi <- CreateBasis(K=3, type='legendre01', pts=pts)
basisDerMat <- apply(phi, 2, function(x)
ConvertSupport(seq(0, 1, length.out=M - 1), pts, diff(x) * (M - 1)))
trueDer <- matrix(1, n, M, byrow=TRUE) + tcrossprod(samp2$xi, basisDerMat)
matplot(t(trueDer), type='l', ylim=c(-5, 10))
# DPC is slightly better in terms of RMSE
mean((fitted(FPC) - trueDer)^2)
mean((fitted(DPC) - trueDer)^2)

34

FPCquantile

FPCquantile

Conditional Quantile estimation with functional covariates

Description
Main function to implement conditional Quantile estimation with functional covariates and scalar
response. The method includes 3 steps: 1) FPCA using the PACE method for X(t_x) 2) Computation
of the conditional distribution function through a functional generalized linear model. 3) Prediction
of quantiles for given predictor values
Usage
FPCquantile(Lx, Lt_x, y, outQ = c(0.1, 0.25, 0.5, 0.75, 0.9),
optns_x = NULL, isNewsub = NULL)
Arguments
Lx

A length n list of predictor function where x[[i]] is the row vector of measurements for ith subject, i=1,...,n

Lt_x

A length n list where the observations of x are taken, t_x[[i]] is a row vector of
time points where x[[i]] are observed, i=1,...,n

y

A 1*n vector for scalar response y. y[i] is the response value for the ith subject,
i = 1,...,n.

outQ

A vector of desired quantile levels with default value outQ = c(0.1, 0.25, 0.5,
0.75, 0.9).

optns_x

A list of options for predictor x with control parameters specified by list(name=value)
with default NA. See function FPCA for details.

isNewSub

A 1*n vector of 0s or 1s, where n is the total count of subjects. 0 denotes the
corresponding subject is only used for training and 1 denotes the corresponding
subject is only used for prediction. (default: 0’s)

Value
A list of the following

pred_quantile: a matrix of n*length(outQ) where the the first nn (number of 0s in isNewSub) rows

pred_CDF: a matrix of n*100. The ith row contains the fitted or predicted conditional distributi

P
b: a matrix of 50*(K+1) contains the coefficient functions, defined as F (y|X) = g( ( k = 0)K bk (y

References
Chen, K., M\"uller, H.G. (2011). Conditional quantile analysis when covariates are functions, with
application to growth data. J. Royal Statistical Society B.

FPCReg

35

Examples
set.seed(10)
n = 200
npred = 50
m = 50
xi <- Wiener(n, 0:m/m)
x=list()
t_x=list()
y=numeric(n)
for(i in 1:n){
t_x = c(t_x,list(0:m/m))
x = c(x,list(xi[i,]))
y[i] = 5*rnorm(1)+2*sum(xi[i,])
}
outQ = c(0.1,0.25,0.5,0.75,0.9,0.95)
isNewsub = c(rep(0,150),rep(1,50))
qtreg = FPCquantile(x, t_x, y, outQ,optns_x = NULL,isNewsub)

FPCReg

Function for performing functonal linear regression where the covariates are functions X1(t1),X2(t2),.. and the response is a function
Y(t_y).

Description
Function for performing functonal linear regression where the covariates are functions X1(t1),X2(t2),..
and the response is a function Y(t_y).
Usage
FPCReg(vars, varsOptns = NULL, isNewSub = NULL, method = "AIC",
FVEthreshold = 0.99, alpha = 0.05, Kx = NULL)
Arguments
vars

A list of input functional covariates with name of "X1", "X2",.. and a functional
response with name "Y". Each field should have two fields: ’Lt’, a list (sparse)
or a matrix (Dense) specifying the time of observations, and ’Ly’, a list (Sparse)
or a matrix (Dense) of the observations.

varsOptns

A list of options named by "X1", "X2",..."Y". Each filed specify the paramaters
that control the corresponding variables. (default: see details of FPCA())

isNewSub

A 1*n vector of 0s or 1s, where n is the total count of subjects. 0 denotes the corresponding subject is only used for estimation and 1 denotes the corresponding
subject is only used for prediction. (default: 0’s)

method

The method used for selecting the number of principal components of functional
predictors X’s used in functional regression , including ’AIC’, ’BIC’ and ’FVE’.
(default: "AIC")

36

FPCReg
FVEthreshold

A scalar specifying the proportion used for ’FVE’. (default: 0.99)

alpha

A scalar specifying the level of the confidence bands. (default: 0.05)

Kx

The number of principal components of functional predictors X’s used in functional regression.

Value
A list containing the following fields:
estiBeta

A list with fields of estimated beta_XiY(s,t) defiend on [range(Xi),range(Y)]

predictY

A list containing fitted or predicted (when is NewSub is true) functions for
E(Y|X).

cbandY

A list with confidence bands of E(Y|X).

Q

Quasi R-square

r2

Functional R-square.

varsMean

A list with mean function of covariates and response.

Kx

The number of principal components of functional predictors X’s used in functional regression.

References
Yao, F., Mueller, H.G., Wang, J.L. "Functional Linear Regression Analysis for Longitudinal Data."
Annals of Statistics 33, (2005): 2873-2903.
Examples
set.seed(1000)
#Model: E(Y(t)|X) = int(beta(s,t)*X(s))
n <- 200 #number of subjects
ngrids <- 51 #number of grids in [0,1] for X(s)
ngridt <- 101 #number of grids in [0,1] for Y(t)
grids <- seq(0, 1, length.out=ngrids) #regular grids in [0,1] for X(s)
gridt <- seq(0, 1, length.out=ngridt) #regular grids in [0,1] for Y(t)
#generate X
#{1, sqrt(2)*sin(2*pi*s), sqrt(2)*cos(2*pi*t)} are used to generate X.
eigenFun <- list( function(s){1 + 0 * s},
function(s){sqrt(2) * sin(2*pi*s)},
function(s){sqrt(2) * cos(2*pi*s)})
sig <- matrix(c(1.5, 0.0, 0.0,
0.0, 1.2, 0.0,
0.0, 0.0, 1.0,
0.9, -.3, 0.4,
-.5, 0.8, -.3,
0.1, 0.4, 0.7,
nrow=6,ncol=6)

0.9,
-.3,
0.4,
2.0,
0.0,
0.0,

-.5,
0.8,
-.3,
0.0,
1.5,
0.0,

0.1,
0.4,
0.7,
0.0,
0.0,
1.0),

scoreX <- MASS::mvrnorm(n,mu=rep(0,6),Sigma=sig)
scoreX1 <- scoreX[,1:3]
scoreX2 <- scoreX[,4:6]
basisX1 <- sapply(eigenFun,function(x){x(grids)})

FPCReg

37

latentX1 <- scoreX1 %*% t(basisX1)
measErrX1 <- sqrt(0.03) * matrix(rnorm(n * ngrids), n, ngrids) #0.01 is sigma^2.
denseX1 <- latentX1 + measErrX1
basisX2 <- sapply(eigenFun,function(x){x(grids)})
latentX2 <- scoreX2 %*% t(basisX2)
measErrX2 <- sqrt(0.03) * matrix(rnorm(n * ngrids), n, ngrids) #0.01 is sigma^2.
denseX2 <- latentX2 + measErrX2
#generate Y
#beta(s, t) <- sin(2 * pi * s)*cos(2 * pi * t)
betaEigen1 <- function(t){f <- function(s){
sin(2*pi*s) * cos(2*pi*t) * (1+0*s)}; return(f)}
betaEigen2 <- function(t){f <- function(s){
sin(2*pi*s) * cos(2*pi*t) * (sqrt(2)*sin(2*pi*s))}; return(f)}
betaEigen3 <- function(t){f <- function(s){
sin(2*pi*s) * cos(2*pi*t) * (sqrt(2)*cos(2*pi*s))}; return(f)}
betaEigen <- list(betaEigen1, betaEigen2, betaEigen3)
basisY <- array(0,c(ngridt, 3))
for(i in 1:3){
intbetaEigen <- function (t) {integrate(betaEigen[[i]](t), lower = 0, upper = 1)$value}
basisY[, i] <- sapply(1:ngridt, function(x){intbetaEigen(gridt[x])})
}
latentY <- scoreX1 %*% t(basisY) - scoreX2 %*% t(basisY)
measErrY <- sqrt(0.01) * matrix(rnorm(n*ngridt), n, ngridt) #0.01 is sigma^2
denseY <- latentY + measErrY
#======Dense data===============================================
timeX <- t(matrix(rep(grids, n),length(grids), n))
timeY <- t(matrix(rep(gridt, n),length(gridt), n))
denseVars <- list(X1 = list(Ly = denseX1, Lt = timeX),
X2 = list(Ly = denseX2, Lt = timeX),
Y=list(Ly = denseY,Lt = timeY))
resuDense <- FPCReg(denseVars, method="FVE")
par(mfrow=c(1,2))
estiBetaX1Y_Dense <- resuDense$estiBeta$betaX1Y
args1 <- list(xlab = 's', ylab = 't', zlab = 'estiBetaX1Y_Dense(s, t)',
lighting = FALSE, phi = 45, theta = 45)
args2 <- list(x = 1:ngrids, y = 1:ngridt, z = estiBetaX1Y_Dense[1:ngrids, 1:ngridt])
do.call(plot3D::persp3D,c(args2, args1))
estiBetaX2Y_Dense <- resuDense$estiBeta$betaX2Y
args1 <- list(xlab = 's', ylab = 't', zlab = 'estiBetaX2Y_Dense(s, t)',
lighting = FALSE, phi = 45, theta = 45)
args2 <- list(x = 1:ngrids, y = 1:ngridt, z = estiBetaX2Y_Dense[1:ngrids, 1:ngridt])
# do.call(plot3D::persp3D,c(args2, args1))
#======Sparse data===============================================
## Not run:
sparsity = 5:8
sparseX1 <- Sparsify(denseX1, grids, sparsity)
sparseX2 <- Sparsify(denseX2, grids, sparsity)
sparseY <- Sparsify(denseY, gridt, sparsity)
sparseVars <- list(X1 = sparseX1, X2 = sparseX2, Y = sparseY)

38

FSVD
resuSparse <- FPCReg(sparseVars, method="FVE", FVEthreshold=0.98)
#or resuSparse <- FPCReg(vars = sparseVars,
#
varsOptns = list(X1=list(userBwCov = 0.03)))
par(mfrow=c(1,2))
estiBetaX1Y_Sparse = resuSparse$estiBeta$betaX1Y
args1 = list(xlab = 's', ylab = 't', zlab = 'estiBetaX1Y_Sparse(s,t)',
lighting = FALSE, phi = 45,theta = 45)
args2 = list(x = 1:51, y = 1:51, z = estiBetaX1Y_Sparse[1:51, 1:51])
do.call(plot3D::persp3D, c(args2, args1))
estiBetaX2Y_Sparse = resuSparse$estiBeta$betaX2Y
args1 = list(xlab = 's', ylab = 't', zlab = 'estiBetaX2Y_Sparse(s,t)',
lighting = FALSE, phi = 45,theta = 45)
args2 = list(x = 1:51, y = 1:51, z = estiBetaX2Y_Sparse[1:51, 1:51])
do.call(plot3D::persp3D, c(args2, args1))
par(mfrow=c(2,3))
for(i in 1:6){
plot(sparseVars[['Y']]$Lt[[i]], sparseVars[['Y']]$Ly[[i]],
xlab = 'time', ylab = 'observations', ylim = c(-1.5, 1.5))
lines(seq(0, 1, length.out = 51), resuSparse$predictY[[i]])
lines(seq(0, 1, length.out = 51), resuSparse$cbandY[[i]][,2], lty = 2)
lines(seq(0, 1, length.out = 51), resuSparse$cbandY[[i]][,1], lty = 2)
}
## End(Not run)

FSVD

Functional Singular Value Decomposition

Description
FSVD for a pair of dense or sparse functional data.
Usage
FSVD(Ly1, Lt1, Ly2, Lt2, FPCAoptns1 = NULL, FPCAoptns2 = NULL,
SVDoptns = list())
Arguments
Ly1

A list of n vectors containing the observed values for each individual. Missing
values specified by NAs are supported for dense case (dataType='dense').

Lt1

A list of n vectors containing the observation time points for each individual
corresponding to y. Each vector should be sorted in ascending order.

Ly2

A list of n vectors containing the observed values for each individual. Missing
values specified by NAs are supported for dense case (dataType='dense').

Lt2

A list of n vectors containing the observation time points for each individual
corresponding to y. Each vector should be sorted in ascending order.

FPCAoptns1

A list of options control parameters specified by list(name=value) for the FPC
analysis of sample 1. See ‘?FPCA’.

FSVD

39

FPCAoptns2

A list of options control parameters specified by list(name=value) for the FPC
analysis of sample 2. See ‘?FPCA’.

SVDoptns

A list of options control parameters specified by list(name=value) for the
FSVD analysis of samples 1 & 2. See ‘Details‘.

Details
Available control options for SVDoptns are:
bw1 The bandwidth value for the smoothed cross-covariance function across the direction of sample 1; positive numeric - default: determine automatically based on ’methodBwCov’
bw2 The bandwidth value for the smoothed cross-covariance function across the direction of sample 2; positive numeric - default: determine automatically based on ’methodBwCov’
methodBwCov The bandwidth choice method for the smoothed covariance function; ’GMeanAndGCV’
(the geometric mean of the GCV bandwidth and the minimum bandwidth),’CV’,’GCV’ - default: 10% of the support
userMu1 The user defined mean of sample 1 used to centre it prior to the cross-covariance estimation. - default: determine automatically based by the FPCA of sample 1
userMu2 The user defined mean of sample 2 used to centre it prior to the cross-covariance estimation. - default: determine automatically based by the FPCA of sample 2
maxK The maximum number of singular components to consider; default: min(20, N-1), N:# of
curves.
kernel Smoothing kernel choice, common for mu and covariance; "rect", "gauss", "epan", "gausvar", "quar" - default: "gauss"; dense data are assumed noise-less so no smoothing is performed.
rmDiag Logical describing if the routine should remove diagonal raw cov for cross cov estimation
(default: FALSE)
noScores Logical describing if the routine should return functional singular scores or not (default:
TRUE)
regulRS String describing if the regularisation of the compositie cross-covariance matrix should
be done using ’sigma1’ or ’rho’ (see ?FPCA for details) (default: ’sigma2’)
bwRoutine String specifying the routine used to find the optimal bandwidth ’grid-search’, ’bobyqa’,
’l-bfgs-b’ (default: ’l-bfgs-b’)
flip Logical describing if the routine should flip the sign of the singular components functions or
not after the SVD of the cross-covariance matrix. (default: FALSE)
useGAM Indicator to use gam smoothing instead of local-linear smoothing (semi-parametric option) (default: FALSE)
dataType1 The type of design we have for sample 1 (usually distinguishing between sparse or
dense functional data); ’Sparse’, ’Dense’, ’DenseWithMV’ - default: determine automatically
based on ’IsRegular’
dataType2 The type of design we have for sample 2 (usually distinguishing between sparse or
dense functional data); ’Sparse’, ’Dense’, ’DenseWithMV’ - default: determine automatically
based on ’IsRegular’
Value
A list containing the following fields:
bw1

The selected (or user specified) bandwidth for smoothing the cross-covariance
function across the support of sample 1.

40

FVPA
bw2

The selected (or user specified) bandwidth for smoothing the cross-covariance
function across the support of sample 2.

CrCov

The smoothed cross-covariance between samples 1 & 2.

sValues

A list of length nsvd, each entry containing the singuar value estimates for the
FSC estimates.

nsvd

The number of singular componentes used.

canCorr

The canonical correlations for each dimension.

FVE

A percentage indicating the total variance explained by chosen FSCs with corresponding ’FVEthreshold’.

sFun1

An nWorkGrid by K matrix containing the estimated singular functions for sample 1.

sFun2

An nWorkGrid by K matrix containing the estimated singular functions for sample 2.

grid1

A vector of length nWorkGrid1. The internal regular grid on which the singular
analysis is carried on the support of sample 1.

grid2

A vector of length nWorkGrid2. The internal regular grid on which the singular
analysis is carried on the support of sample 2.

sScores1

A n by K matrix containing the singular scores for sample 1.

sScores2

A n by K matrix containing the singular scores for sample 2.

optns

A list of options used by the SVD and the FPCA’s procedures.

FVPA

Functional Variance Process Analysis for dense functional data

Description
Functional Variance Process Analysis for dense functional data
Usage
FVPA(y, t, q = 0.1, optns = list(error = TRUE, FVEthreshold = 0.9))
Arguments
y

A list of n vectors containing the observed values for each individual. Missing
values specified by NAs are supported for dense case (dataType='dense').

t

A list of n vectors containing the observation time points for each individual
corresponding to y.

q

A scalar defining the percentile of the pooled sample residual sample used for
adjustment before taking log (default: 0.1).

optns

A list of options control parameters specified by list(name=value); by default:
’error’ has to be TRUE, ’FVEthreshold’ is set to 0.90. See ‘Details in ?FPCA’.

GetCrCorYX

41

Value
A list containing the following fields:
sigma2

Variance estimator of the functional variance process.

fpcaObjY

FPCA object for the original data.

fpcaObjR

FPCA object for the functional variance process associated with the original
data.

References
Hans-Georg Mueller, Ulrich Stadtmuller and Fang Yao, "Functional variance processes." Journal of
the American Statistical Association 101 (2006): 1007-1018
Examples
set.seed(1)
n <- 25
pts <- seq(0, 1, by=0.01)
sampWiener <- Wiener(n, pts)
# Data have to dense for FVPA to be relevant!
sampWiener <- Sparsify(sampWiener, pts, 101)
fvpaObj <- FVPA(sampWiener$Ly, sampWiener$Lt)

GetCrCorYX

Make cross-correlation matrix from auto- and cross-covariance matrix

Description
Make cross-correlation matrix from auto- andcross-covariance matrix
Usage
GetCrCorYX(ccXY, ccXX, ccYY)
Arguments
ccXY

The cross-covariance matrix between variables X and Y.

ccXX

The auto-covariance matrix of variable X or the diagonal of that matrix.

ccYY

The auto-covariance matrix of variable Y or the diagonal of that matrix.

Value
A cross-correlation matrix between variables X and Y.

42

GetCrCovYX

GetCrCorYZ

Make cross-correlation matrix from auto- and cross-covariance matrix

Description
Make cross-correlation matrix from auto- andcross-covariance matrix
Usage
GetCrCorYZ(ccYZ, acYY, covZ)
Arguments
ccYZ

The cross-covariance vector between variables Y and Z (n-by-1).

acYY

The auto-covariance n-by-n matrix of variable Y or the (n-by-1) diagonal of that
matrix.

covZ

The (scalar) covariance of variable Z.

Value
A cross-correlation matrix between variables Y (functional) and Z (scalar).

GetCrCovYX

Functional Cross Covariance between longitudinal variable Y and
longitudinal variable X

Description
Calculate the raw and the smoothed cross-covariance between functional predictors using bandwidth bw or estimate that bw using GCV.
Usage
GetCrCovYX(bw1 = NULL, bw2 = NULL, Ly1, Lt1 = NULL, Ymu1 = NULL,
Ly2, Lt2 = NULL, Ymu2 = NULL, useGAM = FALSE, rmDiag = FALSE,
kern = "gauss", bwRoutine = "l-bfgs-b")
Arguments
bw1

Scalar bandwidth for smoothing the cross-covariance function (if NULL it will
be automatically estimated) (Y)

bw2

Scalar bandwidth for smoothing the cross-covariance function (if NULL it will
be automatically estimated) (X)

Ly1

List of N vectors with amplitude information (Y)

Lt1

List of N vectors with timing information (Y)

Ymu1

Vector Q-1 Vector of length nObsGrid containing the mean function estimate
(Y)

Ly2

List of N vectors with amplitude information (X)

GetCrCovYZ

43

Lt2

List of N vectors with timing information (X)

Ymu2

Vector Q-1 Vector of length nObsGrid containing the mean function estimate
(X)

useGAM

Indicator to use gam smoothing instead of local-linear smoothing (semi-parametric
option) (default: FALSE)

rmDiag

Indicator to remove the diagonal element when smoothing (default: FALSE)

kern

String specifying the kernel type (default: FALSE; see ?FPCA for details)

bwRoutine

String specifying the routine used to find the optimal bandwidth ’grid-search’,
’bobyqa’, ’l-bfgs-b’ (default: ’l-bfgs-b’) If the variables Ly1 and Ly2 are in
matrix form the data are assumed dense and only the raw cross-covariance is
returned. One can obtain Ymu1 and Ymu2 from FPCA and ConvertSupport.

Value
A list containing:
smoothedCC

The smoothed cross-covariance as a matrix (currently only 51 by 51)

rawCC

The raw cross-covariance as a list

bw

The bandwidth used for smoohting as a vector of lengh 2

score

The GCV score associated with the scalar used

smoothGrid

The grid over which the smoothed cross-covariance is evaluated

References
Yang, Wenjing, Hans-Georg Mueller, and Ulrich Stadtmueller. "Functional singular component
analysis." Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73.3 (2011):
303-324
Examples
Ly1= list( rep(2.1,7), rep(2.1,3),2.1 );
Lt1 = list(1:7,1:3, 1);
Ly2 = list( rep(1.1,7), rep(1.1,3),1.1);
Lt2 = list(1:7,1:3, 1);
Ymu1 = rep(55,7);
Ymu2 = rep(1.1,7);
AA<-GetCrCovYX(Ly1 = Ly1, Ly2= Ly2, Lt1=Lt1, Lt2=Lt2, Ymu1=Ymu1, Ymu2=Ymu2)

GetCrCovYZ

Functional Cross Covariance between longitudinal variable Y and
scalar variable Z

Description
Calculate the raw and the smoothed cross-covariance between functional and scalar predictors using
bandwidth bw or estimate that bw using GCV

44

GetCrCovYZ

Usage
GetCrCovYZ(bw = NULL, Z, Zmu = NULL, Ly, Lt = NULL, Ymu = NULL,
support = NULL, kern = "gauss")

Arguments
bw

Scalar bandwidth for smoothing the cross-covariance function (if NULL it will
be automatically estimated)

Z

Vector N-1 Vector of length N with the scalar function values

Zmu

Scalar with the mean of Z (if NULL it will be automaticall estimated)

Ly

List of N vectors with amplitude information

Lt

List of N vectors with timing information

Ymu

Vector Q-1 Vector of length nObsGrid containing the mean function estimate

support

Vector of unique and sorted values for the support of the smoothed cross-covariance
function (if NULL it will be automatically estimated)

kern

Kernel type to be used. See ?FPCA for more details. (defult: ’gauss’) If the variables Ly1 is in matrix form the data are assumed dense and only the raw crosscovariance is returned. One can obtain Ymu1 from FPCA and ConvertSupport.

Value
A list containing:
smoothedCC

The smoothed cross-covariance as a vector

rawCC

The raw cross-covariance as a vector

bw

The bandwidth used for smoohting as a scalar

score

The GCV score associated with the scalar used

References
Yang, Wenjing, Hans-Georg Mueller, and Ulrich Stadtmueller. "Functional singular component
analysis." Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73.3 (2011):
303-324

Examples
Ly <- list( runif(5), c(1:3), c(2:4), c(4))
Lt <- list( c(1:5), c(1:3), c(1:3), 4)
Z = rep(4,4) # Constant vector so the covariance has to be zero.
sccObj = GetCrCovYZ(bw=1, Z= Z, Ly=Ly, Lt=Lt, Ymu=rep(4,5))

GetNormalisedSample

45

GetNormalisedSample

Normalise sparse functional sample

Description
Normalise sparse functional sample given in an FPCA object

Usage
GetNormalisedSample(fpcaObj, errorSigma = FALSE)
GetNormalizedSample(...)
Arguments
fpcaObj

An FPCA object.

errorSigma

Indicator to use sigma^2 error variance when normalising the data (default:
FALSE)

...

Passed into GetNormalisedSample

Value
A list containing the normalised sample ’y’ at times ’t’

References
Chiou, Jeng-Min and Chen, Yu-Ting and Yang, Ya-Fang. "Multivariate Functional Principal Component Analysis: A Normalization Approach" Statistica Sinica 24 (2014): 1571-1596

Examples
set.seed(1)
n <- 100
M <- 51
pts <- seq(0, 1, length.out=M)
mu <- rep(0, length(pts))
sampDense <- MakeGPFunctionalData(n, M, mu, K=1, basisType='sin', sigma=0.01)
samp4 <- MakeFPCAInputs(tVec=sampDense$pts, yVec=sampDense$Yn)
res4E <- FPCA(samp4$Ly, samp4$Lt, list(error=TRUE))
sampN <- GetNormalisedSample(res4E, errorSigma=TRUE)
CreatePathPlot(subset=1:20, inputData=samp4, obsOnly=TRUE, showObs=FALSE)
CreatePathPlot(subset=1:20, inputData=sampN, obsOnly=TRUE, showObs=FALSE)

46

kCFC

kCFC

Functional clustering and identifying substructures of longitudinal
data using kCFC.

Description
Functional clustering and identifying substructures of longitudinal data using kCFC.
Usage
kCFC(y, t, k = 3, kSeed = 123, maxIter = 125,
optnsSW = list(methodMuCovEst = "smooth", FVEthreshold = 0.9,
methodBwCov = "GCV", methodBwMu = "GCV"), optnsCS = list(methodMuCovEst
= "smooth", FVEthreshold = 0.7, methodBwCov = "GCV", methodBwMu = "GCV"))
Arguments
y

A list of n vectors containing the observed values for each individual. Missing
values specified by NAs are supported for dense case (dataType='dense').

t

A list of n vectors containing the observation time points for each individual
corresponding to y.

k

A scalar defining the number of clusters to define; default 3. Values that define
very small clusters (eg.cluster size <=3) will potentiall err.

kSeed

A scalar valid seed number to ensure replication; default: 123

maxIter

A scalar defining the maximum number of iterations allowed; default 20, common for both the simple kmeans initially and the subsequent k-centres

optnsSW

A list of options control parameters specified by list(name=value) to be used
for sample-wide FPCA; by default: "list( methodMuCovEst =’smooth’, FVEthreshold= 0.90, methodBwCov = ’GCV’, methodBwMu = ’GCV’ )". See ‘Details in
?FPCA’.

optnsCS

A list of options control parameters specified by list(name=value) to be used
for cluster-specific FPCA; by default: "list( methodMuCovEst =’smooth’, FVEthreshold= 0.70, methodBwCov = ’GCV’, methodBwMu = ’GCV’ )". See ‘Details in
?FPCA’.

Value
A list containing the following fields:
cluster

A vector of levels 1:k, indicating the cluster to which each curve is allocated.

fpcaList

A list with the fpcaObj for each separate cluster.

iterToConv

A number indicating how many iterations where required until convergence.

References
Jeng-Min Chiou, Pai-Ling Li, "Functional clustering and identifying substructures of longitudinal
data." Journal of the Royal Statistical Society 69 (2007): 679-699

Lwls1D

47

Examples
data(medfly25)
Flies <- MakeFPCAInputs(medfly25$ID, medfly25$Days, medfly25$nEggs)
kcfcObj <- kCFC(Flies$Ly[1:250], Flies$Lt[1:250], # using only 250 for speed consideration
optnsSW = list(methodMuCovEst = 'smooth', userBwCov = 2, FVEthreshold = 0.90),
optnsCS = list(methodMuCovEst = 'smooth', userBwCov = 2, FVEthreshold = 0.70))

Lwls1D

One dimensional local linear kernel smoother

Description
One dimensional local linear kernel smoother for longitudinal data.
Usage
Lwls1D(bw, kernel_type, win = rep(1L, length(xin)), xin, yin, xout,
npoly = 1L, nder = 0L)
Arguments
bw

Scalar holding the bandwidth

kernel_type

Character holding the kernel type (see ?FPCA for supported kernels)

win

Vector of length N with weights

xin

Vector of length N with measurement points

yin

Vector of length N with measurement values

xout

Vector of length M with output measurement points

npoly

Scalar (integer) degree of polynomial fitted (default 1)

nder

Scalar (integer) degree of derivative fitted (default 0)

Value
Vector of length M with measurement values at the the point speficied by ’xout’

Lwls2D

Two dimensional local linear kernel smoother.

Description
Two dimensional local weighted least squares smoother. Only local linear smoother for estimating
the original curve is available (no higher order, no derivative).
Usage
Lwls2D(bw, kern = "epan", xin, yin, win = NULL, xout1 = NULL,
xout2 = NULL, xout = NULL, subset = NULL, crosscov = FALSE,
method = ifelse(kern == "gauss", "plain", "sort2"))

48

Lwls2DDeriv

Arguments
bw

A scalar or a vector of length 2 specifying the bandwidth.

kern

Kernel used: ’gauss’, ’rect’, ’gausvar’, ’epan’ (default), ’quar’.

xin

An n by 2 data frame or matrix of x-coordinate.

yin

A vector of y-coordinate.

win

A vector of weights on the observations.

xout1

a p1-vector of first output coordinate grid. Defaults to the input gridpoints if left
unspecified.

xout2

a p2-vector of second output coordinate grid. Defaults to the input gridpoints if
left unspecified.

xout

alternative to xout1 and xout2. A matrix of p by 2 specifying the output points
(may be inefficient if the size of xout is small).

subset

a vector with the indices of x-/y-/w-in to be used (Default: NULL)

crosscov

using function for cross-covariance estimation (Default: FALSE)

method

should one try to sort the values xin and yin before using the lwls smoother? if
yes (’sort2’ - default for non-Gaussian kernels), if no (’plain’ - fully stable; de)

Value
a p1 by p2 matrix of fitted values if xout is not specified. Otherwise a vector of length p corresponding to the rows of xout.

Lwls2DDeriv

Two dimensional local linear kernel smoother with derivatives.

Description
Two dimensional local weighted least squares smoother. Only local linear smoother for estimating
the original curve is available (no higher order, no derivative).
Usage
Lwls2DDeriv(bw, kern = "epan", xin, yin, win = NULL, xout1 = NULL,
xout2 = NULL, xout = NULL, npoly = 1L, nder1 = 0L, nder2 = 0L,
subset = NULL, crosscov = TRUE, method = "sort2")
Arguments
bw

A scalar or a vector of length 2 specifying the bandwidth.

kern

Kernel used: ’gauss’, ’rect’, ’gausvar’, ’epan’ (default), ’quar’.

xin

An n by 2 data frame or matrix of x-coordinate.

yin

A vector of y-coordinate.

win

A vector of weights on the observations.

xout1

a p1-vector of first output coordinate grid. Defaults to the input gridpoints if left
unspecified.

MakeBWtoZscore02y

49

xout2

a p2-vector of second output coordinate grid. Defaults to the input gridpoints if
left unspecified.

xout

alternative to xout1 and xout2. A matrix of p by 2 specifying the output points
(may be inefficient if the size of xout is small).

npoly

The degree of polynomials (include all xa y b terms where a + b <= npoly)

nder1

Order of derivative in the first direction

nder2

Order of derivative in the second direction

subset

a vector with the indices of x-/y-/w-in to be used (Default: NULL)

crosscov

using function for cross-covariance estimation (Default: TRUE)

method

should one try to sort the values xin and yin before using the lwls smoother? if
yes (’sort2’ - default for non-Gaussian kernels), if no (’plain’ - fully stable; de)

Value
a p1 by p2 matrix of fitted values if xout is not specified. Otherwise a vector of length p corresponding to the rows of xout.

MakeBWtoZscore02y

Z-score body-weight for age 0 to 24 months based on WHO standards

Description
Make vector of age and body-weight to z-scores based on WHO standards using LMS

Usage
MakeBWtoZscore02y(sex, age, bw)
Arguments
sex

A character ’M’ or ’F’ indicating the sex of the child.

age

A vector of time points of size Q.

bw

A vector of body-weight readings of size Q.

Value
A vector of Z-scores of size Q.

50

MakeGPFunctionalData

MakeFPCAInputs

Format FPCA input

Description
Turn vector inputs to the list so they can be used in FPCA
Usage
MakeFPCAInputs(IDs = NULL, tVec, yVec, na.rm = FALSE, sort = FALSE)
Arguments
IDs

n-by-1 vector of subject IDs (Default: NULL)

tVec

Either an n-by-1 vector of measurement times, or a p-by-1 vector corresponding
to the common time support

yVec

n-by-1 vector of measurements from the variable of interest, or a n-by-p matrix
with each row corresponding to the dense observations.

na.rm

logical indicating if NA should be omitted (Default: FALSE)

sort

logical indicating if the returned lists Lt and Ly should be ensured to be sorted
(Default: FALSE)

Value
L list containing 3 lists each of length ’m’, ’m’ being the number of unique subject IDs

MakeGPFunctionalData

Make Gaussian Process Dense Functional Data sample

Description
Make a Gaussian process dense functional data sample of size n over a [0,1] support.
Usage
MakeGPFunctionalData(n, M = 100, mu = rep(0, M), K = 2,
lambda = rep(1, K), sigma = 0, basisType = "cos")
Arguments
n

number of samples to generate

M

number of equidistant readings per sample (default: 100)

mu

vector of size M specifying the mean (default: rep(0,M))

K

scalar specifying the number of basis to be used (default: 2)

lambda

vector of size K specifying the variance of each components (default: rep(1,K))

sigma

The standard deviation of the Gaussian noise added to each observation points.

basisType

string specifiying the basis type used; possible options are: ’sin’, ’cos’ and
’fourier’ (default: ’cos’) (See code of ’CreateBasis’ for implementation details.)

MakeHCtoZscore02y

51

Value
Y: X(t_j), Yn: noisy observations

MakeHCtoZscore02y

Z-score head-circumference for age 0 to 24 months based on WHO
standards

Description
Make vector of age and height measurement to z-scores based on WHO standards using mu and
sigma (not LMS)
Usage
MakeHCtoZscore02y(sex, age, hc)
Arguments
sex

A character ’M’ or ’F’ indicating the sex of the child.

age

A vector of time points of size Q.

hc

A vector of head circumference readings of size Q (in cm).

Value
A vector of Z-scores of size Q.

MakeLNtoZscore02y

Z-score height for age 0 to 24 months based on WHO standards

Description
Make vector of age and height measurement to z-scores based on WHO standards using mu and
sigma (not LMS)
Usage
MakeLNtoZscore02y(sex, age, ln)
Arguments
sex

A character ’M’ or ’F’ indicating the sex of the child.

age

A vector of time points of size Q.

ln

A vector of body-length readings of size Q (in cm).

Value
A vector of Z-scores of size Q.

52

medfly25

MakeSparseGP

Make Gaussian Process Sparse Functional Data sample

Description
Make a Gaussian process sparse functional data sample of size n
Usage
MakeSparseGP(n, rdist = runif, sparsity = 2:9, muFun = function(x)
rep(0, length(x)), K = 2, lambda = rep(1, K), sigma = 0,
basisType = "cos", CovFun = NULL)
Arguments
n
rdist
sparsity
muFun
K
lambda
sigma
basisType
CovFun

number of samples to generate.
a sampler for generating the random design time points within [0, 1].
A vector of integers. The number of observation per sample is chosen to be one
of the elements in sparsity with equal chance.
a function that takes a vector input and output a vector of the corresponding
mean (default: zero function).
scalar specifying the number of basis to be used (default: 2).
vector of size K specifying the variance of each components (default: rep(1,K)).
The standard deviation of the Gaussian noise added to each observation points.
string specifiying the basis type used; possible options are: ’sin’, ’cos’ and
’fourier’ (default: ’cos’) (See code of ’CreateBasis’ for implementation details.)
an alternative specification of the covariance structure.

Value
TODO
medfly25

Number of eggs laid daily from medflies

Description
A dataset containing the eggs laid from 789 medflies (Mediterranean fruit flies, Ceratitis capitata)
during the first 25 days of their lives. This is a subset of the dataset used by Carey at al. (1998);
only flies having lived at least 25 days are shown. At the end of the recording period all flies were
still alive.
Format
A data frame with 19725 rows and 3 variables:
ID : Medfly ID according to the orignal dataset
Days : Day of measurement
nEggs : Number of eggs laid at that particular day
nEggsRemain : Remaining total number of eggs laid

MultiFAM

53

Source
http://anson.ucdavis.edu/~mueller/data/medfly1000.html
References
Carey, J.R., Liedo, P., Mueller, H.G., Wang, J.L., Chiou, J.M. (1998). Relationship of age patterns
of fecundity to mortality, longevity, and lifetime reproduction in a large cohort of Mediterranean
fruit fly females. J. of Gerontology –Biological Sciences 53, 245-251.

MultiFAM

Functional Additive Models with Multiple Predictor Processes

Description
Smooth backfitting procedure for functional additive models with multiple predictor processes
Usage
MultiFAM(Y, X, ker = "epan", nEval = 51, XTest = NULL,
bwMethod = 0, alpha = 0.7, supp = c(-2, 2), optnsList = NULL)
Arguments
Y

An n-dimensional vector whose elements consist of scalar responses.

X

A d-dimensional list whose components consist of two lists of Ly and Lt containing obervation times and functional covariate values for each predictor component, respectively. For details of Ly and Lt, see FPCA for detail.

ker

A function object representing the base kernel to be used in the smooth backfitting algorithm (default is ’epan’ which is the only option supported currently).

nEval

The number of evaluation grid points for kernel smoothing (default is 51. If
it is specified as 0, then estimated FPC scores in the training set are used for
evaluation grid instead of equal grid).

XTest

A d-dimensional list for test set of functional predictors (default is NULL). If
XTest is specified, then estimated FPC scores in the test set are used for evalution grid.

bwMethod

The method of initial bandwidth selection for kernel smoothing, a positive value
for designating K-fold cross-validtaion and zero for GCV (default is 0)

alpha

The shrinkage factor (positive number) for bandwidth selection. See Han et al.
(2016) (default is 0.7).

supp

The lower and upper limits of kernel smoothing domain for studentized FPC
scores, which FPC scores are divided by the square roots of eigenvalues (default
is [-2,2]).

optnsList

A d-dimensional list whose components consist of optns for each predictor
component, respectively. (default is NULL which assigns the same default
optns for all components as in FPCA).

54

MultiFAM

Details
MultiFAM fits functional additive models for a scalar response and multiple predictor processes
based on the smooth backfitting algorithm proposed by Han et al. (2016) that
E(Y |X) =

Kj
d X
X

gjk (ξjk ),

j=1 k=1

where ξjk stand for the k-th FPC score of the j-th predictor process. MultiFAM only focuses on
mutiple predictor processes case. In fact, the case of univariate predictor is the same with functional
additive model proposed by Mueller and Yao (2008). Especially in this development, one can
designate an estimation support of additive models when the additive modeling is only allowed
over restricted intervals or one is interested in the modeling over the support (see Han et al., 2016).
Value
A list containing the following fields:
mu

A scalar of centered part of the regression model.

SBFit

An N by (K1 + ... + Kd ) matrix whose column vectors consist of the smooth
backfitting component function estimators at the given N estimation points.

xi

An N by (K1 + ... + Kd ) matrix whose column vectors consist of FPC score
grid vectors that each additive component functional is evluated.

bw

A (K1 + ... + Kd )-dimensional bandwidth vector.

lambda

A (K1 + ... + Kd )-dimensional vector containing eigenvalues.

phi

A d-dimensional list whose components consist of an nWorkGrid by K_j matrix
containing eigenfunctions, supported by WorkGrid. See FPCA.

workGrid

A d-dimensional list whose components consist of an nWorkGrid by K_j working grid, the internal regular grid on which the eigen analysis is carried on. See
FPCA.

References
Mammen, E., Linton, O. and Nielsen, J. (1999), "The existence and asymptotic properties of a
backfitting projection algorithm under weak conditions", Annals of Statistics, Vol.27, No.5, p.14431490.
Mammen, E. and Park, B. U. (2006), "A simple smooth backfitting method for additive models",
Annals of Statistics, Vol.34, No.5, p.2252-2271.
Mueller, H.-G. and Yao, F. (2008), "Functional additive models", Journal of the Americal Statistical
Association, Vol.103, No.484, p.1534-1544.
Han, K., Mueller, H.-G. and Park, B. U. (2016), "Smooth backfitting for additive modeling with
small errors-in-variables, with an application to additive functional regression for multiple predictor
functions", Bernoulli (accepted).
Examples
set.seed(1000)
library(MASS)
f11 <- function(t) t

MultiFAM

55

f12 <- function(t) 2*cos(2*pi*t/4)
f21 <- function(t) 1.5*sin(2*pi*t/4)
f22 <- function(t) 1.5*atan(2*pi*t/4)
n<-100
N<-200
sig <- matrix(c(2.0, 0.0, 0.5,
0.0, 1.2, -.2,
0.5, -.2, 1.7,
-.2, 0.3, 0.0,
nrow=4,ncol=4)

-.2,
0.3,
0.0,
1.0),

scoreX <- mvrnorm(n,mu=rep(0,4),Sigma=sig)
scoreXTest <- mvrnorm(N,mu=rep(0,4),Sigma=sig)
Y <- f11(scoreX[,1]) + f12(scoreX[,2]) + f21(scoreX[,3]) + f22(scoreX[,4]) + rnorm(n,0,0.5)
YTest <- f11(scoreXTest[,1]) + f12(scoreXTest[,2]) +
f21(scoreXTest[,3]) + f22(scoreXTest[,4]) + rnorm(N,0,0.5)
phi11
phi12
phi21
phi22

<<<<-

function(t)
function(t)
function(t)
function(t)

sqrt(2)*sin(2*pi*t)
sqrt(2)*sin(4*pi*t)
sqrt(2)*cos(2*pi*t)
sqrt(2)*cos(4*pi*t)

grid <- seq(0,1,length.out=21)
Lt <- Lx1 <- Lx2 <- list()
for (i in 1:n) {
Lt[[i]] <- grid
Lx1[[i]] <- scoreX[i,1]*phi11(grid) + scoreX[i,2]*phi12(grid) + rnorm(1,0,0.01)
Lx2[[i]] <- scoreX[i,3]*phi21(grid) + scoreX[i,4]*phi22(grid) + rnorm(1,0,0.01)
}
LtTest <- Lx1Test <- Lx2Test <- list()
for (i in 1:N) {
LtTest[[i]] <- grid
Lx1Test[[i]] <- scoreXTest[i,1]*phi11(grid) + scoreXTest[i,2]*phi12(grid) + rnorm(1,0,0.01)
Lx2Test[[i]] <- scoreXTest[i,3]*phi21(grid) + scoreXTest[i,4]*phi22(grid) + rnorm(1,0,0.01)
}
X1 <- list(Ly=Lx1, Lt=Lt)
X2 <- list(Ly=Lx2, Lt=Lt)
X1Test <- list(Ly=Lx1Test, Lt=LtTest)
X2Test <- list(Ly=Lx2Test, Lt=LtTest)
X <- list(X1, X2)
XTest <- list(X1Test, X2Test)
# estimation
sbf <- MultiFAM(Y=Y,X=X)
xi <- sbf$xi
par(mfrow=c(2,2))
j <- 1
p0 <- trapzRcpp(sort(xi[,j]),dnorm(sort(xi[,j]),0,sqrt(sig[j,j])))

56

MultiFAM
g11 <- f11(sort(xi[,j])) trapzRcpp(sort(xi[,j]),f11(sort(xi[,j]))*dnorm(sort(xi[,j]),0,sqrt(sig[j,j])))/p0
tmpSgn <- sign(sum(g11*sbf$SBFit[,j]))
plot(sort(xi[,j]),g11,type='l',col=2,ylim=c(-2.5,2.5),xlab='xi11')
points(sort(xi[,j]),tmpSgn*sbf$SBFit[order(xi[,j]),j],type='l')
legend('top',c('true','SBF'),col=c(2,1),lwd=2,bty='n',horiz=TRUE)
j <- 2
p0 <- trapzRcpp(sort(xi[,j]),dnorm(sort(xi[,j]),0,sqrt(sig[j,j])))
g12 <- f12(sort(xi[,j])) trapzRcpp(sort(xi[,j]),f12(sort(xi[,j]))*dnorm(sort(xi[,j]),0,sqrt(sig[j,j])))/p0
tmpSgn <- sign(sum(g12*sbf$SBFit[,j]))
plot(sort(xi[,j]),g12,type='l',col=2,ylim=c(-2.5,2.5),xlab='xi12')
points(sort(xi[,j]),tmpSgn*sbf$SBFit[order(xi[,j]),j],type='l')
legend('top',c('true','SBF'),col=c(2,1),lwd=2,bty='n',horiz=TRUE)
j <- 3
p0 <- trapzRcpp(sort(xi[,j]),dnorm(sort(xi[,j]),0,sqrt(sig[j,j])))
g21 <- f21(sort(xi[,j])) trapzRcpp(sort(xi[,j]),f21(sort(xi[,j]))*dnorm(sort(xi[,j]),0,sqrt(sig[j,j])))/p0
tmpSgn <- sign(sum(g21*sbf$SBFit[,j]))
plot(sort(xi[,j]),g21,type='l',col=2,ylim=c(-2.5,2.5),xlab='xi21')
points(sort(xi[,j]),tmpSgn*sbf$SBFit[order(xi[,j]),j],type='l')
legend('top',c('true','SBF'),col=c(2,1),lwd=2,bty='n',horiz=TRUE)
j <- 4
p0 <- trapzRcpp(sort(xi[,j]),dnorm(sort(xi[,j]),0,sqrt(sig[j,j])))
g22 <- f22(sort(xi[,j])) trapzRcpp(sort(xi[,j]),f22(sort(xi[,j]))*dnorm(sort(xi[,j]),0,sqrt(sig[j,j])))/p0
tmpSgn <- sign(sum(g22*sbf$SBFit[,j]))
plot(sort(xi[,j]),g22,type='l',col=2,ylim=c(-2.5,2.5),xlab='xi22')
points(sort(xi[,j]),tmpSgn*sbf$SBFit[order(xi[,j]),j],type='l')
legend('top',c('true','SBF'),col=c(2,1),lwd=2,bty='n',horiz=TRUE)
## Not run:
# fitting
sbf <- MultiFAM(Y=Y,X=X,nEval=0)
yHat <- sbf$mu+apply(sbf$SBFit,1,'sum')
par(mfrow=c(1,1))
plot(yHat,Y)
abline(coef=c(0,1),col=2)
# R^2
R2 <- 1-sum((Y-yHat)^2)/sum((Y-mean(Y))^2)
R2
# prediction
sbf <- MultiFAM(Y=Y,X=X,XTest=XTest)
yHat <- sbf$mu+apply(sbf$SBFit,1,'sum')
par(mfrow=c(1,1))
plot(yHat,YTest)
abline(coef=c(0,1),col=2)
## End(Not run)

NormCurvToArea

NormCurvToArea

57

Normalise a curve to a particular area.

Description
Normalise a curve such that \intyNewdx = area (according to trapezoid integration)
Usage
NormCurvToArea(y, x = seq(0, 1, length.out = length(y)), area = 1)
Arguments
y

values of curve at time-points x

x

design time-points (default: seq(0,1, length.out=length(y)))

area

value to normalise the curve onto (default: 1)

Value
yNew values of curve at times x such that [\intyNewdx = area]

predict.FPCA

Predict FPC scores for a new sample given an FPCA object

Description
Return a matrix with the first K FPC scores according to conditional expectation or numerical
integration.
Usage
## S3 method for class 'FPCA'
predict(object, newLy, newLt, sigma2 = NULL, K = 1,
xiMethod = "CE", ...)
Arguments
object

An FPCA object.

newLy

A list of n vectors containing the observed values for each individual.

newLt

A list of n vectors containing the observation time points for each individual
corresponding to y.

sigma2

The user-defined measurement error variance. A positive scalar. (default: rho if
applicable, otherwise sigma2 if applicable, otherwise 0 if no error. )

K

The scalar defining the number of clusters to define; (default: 1).

xiMethod

The integration method used to calculate the functional principal component
scores ( standard numerical integration ’IN’ or conditional expectation ’CE’);
default: ’CE’.

...

Not used.

58

print.FPCA

Value
A matrix of size n-by-K
Examples
set.seed(1)
n <- 50
pts <- seq(0, 1, by=0.05)
# The first n samples are for training and the rest testing
sampWiener <- Wiener(2 * n, pts)
sparsity <- 2:5
train <- Sparsify(sampWiener[seq_len(n), , drop=FALSE], pts, sparsity)
test <- Sparsify(sampWiener[seq(n + 1, 2 * n), , drop=FALSE], pts, sparsity)
res <- FPCA(train$Ly, train$Lt)
pred <- predict(res, test$Ly, test$Lt, K=5)
plot(pred[, 1], pred[, 2])

print.FPCA

Print an FPCA object

Description
Print a simple description of an FPCA object
Usage
## S3 method for class 'FPCA'
print(x, ...)
Arguments
x

An FPCA object.

...

Not used.

Examples
set.seed(1)
n <- 20
pts <- seq(0, 1, by=0.05)
sampWiener <- Wiener(n, pts)
sampWiener <- Sparsify(sampWiener, pts, 10)
res <- FPCA(sampWiener$Ly, sampWiener$Lt)
res

print.FSVD

print.FSVD

59

Print an FSVD object

Description
Print a simple description of an FSVD object
Usage
## S3 method for class 'FSVD'
print(x, ...)
Arguments
x

An FSVD object.

...

Not used.

print.WFDA

Print a WFDA object

Description
Print a simple description of a WFDA object
Usage
## S3 method for class 'WFDA'
print(x, ...)
Arguments
x

A WFDA object.

...

Not used.

SBFitting

Iterative Smooth Backfitting Algorithm

Description
Smooth backfitting procedure for nonparametric additive models
Usage
SBFitting(Y, x, X, h = NULL, K = "epan", supp = NULL)

60

SBFitting

Arguments
Y

An n-dimensional vector whose elements consist of scalar responses.

x

An N by d matrix whose column vectors consist of N vectors of estimation
points for each component function.

X

An n by d matrix whose row vectors consist of multivariate predictors.

h

A d vector of bandwidths for kernel smoothing to estimate each component
function.

K

A function object representing the kernel to be used in the smooth backfitting
(default is ’epan’, the the Epanechnikov kernel.).

supp

A d by 2 matrix whose row vectors consist of the lower and upper limits of
estimation intervals for each component function (default is the d-dimensional
unit rectangle of [0,1]).

Details
SBFitting fits component functions of additive models for a scalar response and a multivariate
predictor based on the smooth backfitting algorithm proposed by Mammen et al. (1999) and intensively studied by Mammen and Park (2006), Yu et al. (2008), Lee et al. (2010, 2012) and so
on. SBFitting only focuses on the local constant smooth backfitting estimator with multivariate
predictor case. In fact, the case of univariate predictor is the same with the local constant kernel regression estimator (Nadaraya-Watson estimator) and the local polynomial version can be extended
similarly, so that those are omitted in the development. Support of the multivariate predictor is
assumed to be a product of closed intervals. Especially in this development, one can designate an
estimation support of additive models when the additive modeling is only allowed over restricted
intervals or one is interested in the modeling over the support (see Han et al., 2016). If one puts
X on the argument of estimation points x, SBFitting returns estimated values of conditional mean
responses given observed predictors.
Value
A list containing the following fields:
SBFit

An N by d matrix whose column vectors consist of the smooth backfitting component function estimators at the given estimation points.

mY

A scalar of centered part of the regression model.

NW

An N by d matrix whose column vectors consist of the Nadaraya-Watson marginal
regression function estimators for each predictor component at the given estimation points.

mgnDens

An N by d matrix whose column vectors consist of the marginal kernel density
estimators for each predictor component at the given estimation points.

jntDens

An N by N by d by d array representing the 2-dimensional joint kernel density
estimators for all pairs of predictor components at the given estimation grid. For
example, [,,j,k] of the object provides the 2-dimensional joint kernel density
estimator of the (j,k)-component of predictor components at the corresponding
N by N matrix of estimation grid.

itemNum

The iteration number that the smooth backfitting algorithm has stopped.

itemErr

The iteration error of the smooth backfitting algorithm that represents the maximum L2 distance among component functions in the last successive updates.

SBFitting

61

References
Mammen, E., Linton, O. and Nielsen, J. (1999), "The existence and asymptotic properties of a
backfitting projection algorithm under weak conditions", Annals of Statistics, Vol.27, No.5, p.14431490.
Mammen, E. and Park, B. U. (2006), "A simple smooth backfitting method for additive models",
Annals of Statistics, Vol.34, No.5, p.2252-2271.
Yu, K., Park, B. U. and Mammen, E. (2008), "Smooth backfitting in generalized additive models",
Annals of Statistics, Vol.36, No.1, p.228-260.
Lee, Y. K., Mammen, E. and Park., B. U. (2010), "backfitting and smooth backfitting for additive
quantile models", Vol.38, No.5, p.2857-2883.
Lee, Y. K., Mammen, E. and Park., B. U. (2012), "Flexible generalized varying coefficient regression models", Annals of Statistics, Vol.40, No.3, p.1906-1933.
Han, K., Mueller, H.-G. and Park, B. U. (2016), "Smooth backfitting for additive modeling with
small errors-in-variables, with an application to additive functional regression for multiple predictor
functions", Bernoulli (accepted).
Examples
set.seed(100)
n <- 100
d <- 2
X <- pnorm(matrix(rnorm(n*d),nrow=n,ncol=d)%*%matrix(c(1,0.6,0.6,1),nrow=2,ncol=2))
f1 <- function(t) 2*(t-0.5)
f2 <- function(t) sin(2*pi*t)
Y <- f1(X[,1])+f2(X[,2])+rnorm(n,0,0.1)
#
N
x
h

component function estimation
<- 101
<- matrix(rep(seq(0,1,length.out=N),d),nrow=N,ncol=d)
<- c(0.12,0.08)

sbfEst <- SBFitting(Y,x,X,h)
fFit <- sbfEst$SBFit
par(mfrow=c(1,2))
plot(x[,1],f1(x[,1]),type='l',lwd=2,col=2,lty=4,xlab='X1',ylab='Y')
points(x[,1],fFit[,1],type='l',lwd=2,col=1)
points(X[,1],Y,cex=0.3,col=8)
legend('topleft',legend=c('SBF','true'),col=c(1,2),lwd=2,lty=c(1,4),horiz=FALSE,bty='n')
abline(h=0,col=8)
plot(x[,2],f2(x[,2]),type='l',lwd=2,col=2,lty=4,xlab='X2',ylab='Y')
points(x[,2],fFit[,2],type='l',lwd=2,col=1)
points(X[,2],Y,cex=0.3,col=8)
legend('topright',legend=c('SBF','true'),col=c(1,2),lwd=2,lty=c(1,4),horiz=FALSE,bty='n')
abline(h=0,col=8)
# prediction
x <- X
h <- c(0.12,0.08)

62

SelectK
sbfPred <- SBFitting(Y,X,X,h)
fPred <- sbfPred$mY+apply(sbfPred$SBFit,1,'sum')
par(mfrow=c(1,1))
plot(fPred,Y,cex=0.5,xlab='SBFitted values',ylab='Y')
abline(coef=c(0,1),col=8)

SelectK

Selects number of functional principal components for given FPCA
output and selection criteria

Description
Selects number of functional principal components for given FPCA output and selection criteria
Usage
SelectK(fpcaObj, criterion = "FVE", FVEthreshold = 0.95, Ly = NULL,
Lt = NULL)
Arguments
fpcaObj

A list containing FPCA related subjects returned by MakeFPCAResults().

criterion

A string or positive integer specifying selection criterion for number of functional principal components, available options: ’FVE’, ’AIC’, ’BIC’, or the
specified number of components - default: ’FVE’

FVEthreshold

A threshold percentage specified by user when using "FVE" as selection criterion: (0,1] - default: NULL

Ly

A list of n vectors containing the observed values for each individual - default:
NULL

Lt

A list of n vectors containing the observation time points for each individual
corresponding to Ly - default: NULL

Value
A list including the following two fields:
K

An integer indicating the selected number of components based on given criterion.

criterion

The calculated criterion value for the selected number of components, i.e. FVE,
AIC or BIC value, NULL for fixedK criterion.

SetOptions

63

SetOptions

Set the PCA option list

Description
Set the PCA option list
Usage
SetOptions(y, t, optns)
Arguments
y
t
optns

Sparsify

A list of n vectors containing the observed values for each individual.
A list of n vectors containing the observation time points for each individual
corresponding to y.
A list of options control parameters specified by list(name=value). See ‘Details’.
See ’?FPCAfor more details. Usually users are not supposed to use this function
directly.

Sparsify densely observed functional data

Description
Given a matrix of densely observed functional data, make a sparsified sample.
Usage
Sparsify(samp, pts, sparsity, aggressive = FALSE, fragment = FALSE)
Arguments
samp
pts
sparsity
aggressive
fragment

A matrix of densely observed functional data, with each row containing one
sample.
A vector of grid points corresponding to the columns of samp.
A vector of integers. The number of observation per sample is chosen to be one
of the elements in sparsity with equal chance.
Sparsify in an "aggressive" manner making sure that near-by readings are excluded.
Sparsify the observations into fragments, which are (almost) uniformly distributed in the time domain. Default to FALSE as not fragmenting. Otherwise a
positive number specifying the approximate length of each fragment.

Value
A list of length 2, containing the following fields:
Lt
Ly

A list of observation time points for each sample.
A list of values for each sample, corresponding to the time points.

64

Stringing

str.FPCA

Compactly display the structure of an FPCA object

Description
Compactly display the structure of an FPCA object
Usage
## S3 method for class 'FPCA'
str(object, ...)
Arguments
object

An FPCA object

...

Not used

Stringing

Stringing for High-Dimensional data

Description
Stringing for High-Dimensional data
Usage
Stringing(X, Y = NULL, standardize = FALSE, disOptns = "euclidean",
disMat = NA)
Arguments
X

A matrix (n by p) of data, where X[i,] is the row vector of measurements for the
ith subject.

Y

A vector (n by 1), where Y[i] is the reponse associated with X[i,]

standardize

A logical variable indicating whether standardization of the input data matrix is
required, with default: FALSE.

disOptns

A distance metric to be used, one of the following: "euclidean" (default), "maximum", "manhattan", "canberra", "binary", "minkowski", "correlation", "spearman", "hamming", "xycor", or "user". If specified as "xycor", the absolute difference of correlation between predictor and response is used. If specified as
"user", a dissimilarity matrix for the argument "disMat" must be provided.

disMat

A user-specified dissimilarity matrix, only necessary when "disOptns" is "user".

trapzRcpp

65

Value
A list containing the following fields:
Ly

A list of n vectors, which are the random trajectories for all subjects identified
by the Stringing method.

Lt

A list of n time points vectors, at which corresponding measurements Ly are
taken.

StringingOrder A vector representing the order of the stringing, s.t. using as column index on X
yields recovery of the underlying process.
Xin

A matrix, corresponding to the input data matrix.

Xstd

A matrix, corresponding to the standardized input data matrix. It is NULL if
standardize is FALSE.

References
Chen, K., Chen, K., Mueller, H. G., and Wang, J. L. (2011). "Stringing high-dimensional data for
functional analysis." Journal of the American Statistical Association, 106(493), 275-284.
Examples
set.seed(1)
n <- 50
wiener = Wiener(n = n)[,-1]
p = ncol(wiener)
rdmorder = sample(size = p, x=1:p, replace = FALSE)
stringingfit = Stringing(X = wiener[,rdmorder], disOptns = "correlation")
diff_norev = sum(abs(rdmorder[stringingfit$StringingOrder] - 1:p))
diff_rev = sum(abs(rdmorder[stringingfit$StringedOrder] - p:1))
if(diff_rev <= diff_norev){
stringingfit$StringingOrder = rev(stringingfit$StringingOrder)
stringingfit$Ly = lapply(stringingfit$Ly, rev)
}
plot(1:p, rdmorder[stringingfit$StringingOrder], pch=18); abline(a=0,b=1)

trapzRcpp

Trapezoid Rule Numerical Integration

Description
Trapezoid Rule Numerical Integration using Rcpp
Usage
trapzRcpp(X, Y)
Arguments
X

Sorted vector of X values

Y

Vector of Y values.

66

TVAM

TVAM

Iterative Smooth Backfitting Algorithm

Description
Smooth backfitting procedure for time-varying additive models
Usage
TVAM(Lt, Ly, LLx, gridT = NULL, x = NULL, ht = NULL, hx = NULL,
K = "epan", suppT = NULL, suppX = NULL)
Arguments
Lt

An n-dimensional list of N_i-dimensional vector whose elements consist of longitudial time points for each i-th subject.

Ly

An n-dimensional list of N_i-dimensional vector whose elements consist of longitudial response observations of each i-subject corresponding to Lt.

LLx

A tuple of d-lists, where each list respresents longitudinal covariate observations
of the j-th component corresponding to Lt and Ly.

gridT

An M-dimensional sequence of evaluation time points for additive surface estimators. (Must be sorted in increasing orders.)

x

An N by d matrix whose column vectors consist of N vectors of evaluation
points for additive surface component estimators at each covariate value.

ht

A bandwidth for kernel smoothing in time component.

hx

A d vector of bandwidths for kernel smoothing covariate components, respectively.

K

A function object representing the kernel to be used in the smooth backfitting
(default is ’epan’, the the Epanechnikov kernel.).

suppT

A 2-dimensional vector consists of the lower and upper limits of estimation
intervals for time component (default is [0,1]).

suppX

A d by 2 matrix whose row vectors consist of the lower and upper limits of
estimation intervals for each component function (default is the d-dimensional
unit rectangle of [0,1]).

Details
TVAM estimates component surfaces of time-varying additive models for londitudinal observations
based on the smooth backfitting algorithm proposed by Zhang et al. (2013). TVAM only focuses
on the local constant smooth backfitting in contrast to the original development as in Zhang et al.
(2013). However, the local polynomial version can be extended similarly, so that those are omitted
in the development. Especially in this development, one can designate an estimation support of
additive surfaces when the additive modeling is only allowed over restricted intervals or one is
interested in the modeling over the support (see Han et al., 2016).

TVAM

67

Value
A list containing the following fields:
tvamComp

A tuple of d-lists, where each list is given by M by N matrix whose elements
represents the smooth backfitting surface estimator of the j-component evaluated
at gridT and the j-th column of x.

tvamMean

An M-dimensional vector whose elelments consist of the marginal time regression function estimated at gridT.

References
Zhang, X., Park, B. U. and Wang, J.-L. (2013), "Time-varying additive models for longitudinal
data", Journal of the American Statistical Association, Vol.108, No.503, p.983-998.
Han, K., Mueller, H.-G. and Park, B. U. (2018), "Smooth backfitting for additive modeling with
small errors-in-variables, with an application to additive functional regression for multiple predictor
functions", Bernoulli, Vol.24, No.2, p.1233-1265.
Examples
set.seed(1000)
n <- 100
Lt <- list()
Ly <- list()
Lx1 <- list()
Lx2 <- list()
for (i in 1:n) {
Ni <- sample(10:15,1)
Lt[[i]] <- sort(runif(Ni,0,1))
Lx1[[i]] <- runif(Ni,0,1)
Lx2[[i]] <- runif(Ni,0,1)
Ly[[i]] <- Lt[[i]]*(cos(2*pi*Lx1[[i]]) + sin(2*pi*Lx2[[i]])) + rnorm(Ni,0,0.1)
}
LLx <- list(Lx1,Lx2)
gridT <- seq(0,1,length.out=41)
x0 <- seq(0,1,length.out=51)
x <- cbind(x0,x0)
ht <- 0.1
hx <- c(0.1,0.1)
tvam <- TVAM(Lt,Ly,LLx,gridT=gridT,x=x,ht=ht,hx=hx,K='epan')
g0Sbf <- tvam$tvamMean
gjSbf <- tvam$tvamComp
par(mfrow=c(1,2))
par(mar=c(1,1,1,1)+0.1)
persp(gridT,x0,gjSbf[[1]],theta=60,phi=30,

68

VCAM
xlab='time',ylab='x1',zlab='g1(t, x1)')
persp(gridT,x0,gjSbf[[2]],theta=60,phi=30,
xlab='time',ylab='x2',zlab='g1(t, x2)')

VCAM

Sieve estimation

Description
B-spline based estimation procedure for time-varying additive models.
Usage
VCAM(Lt, Ly, X, optnAdd = list(), optnVc = list())
Arguments
Lt

An n-dimensional list of N_i-dimensional vectors whose elements consist of
longitudial time points for each i-th subject.

Ly

An n-dimensional list of N_i-dimensional vectors whose elements consist of
longitudial response observations of each i-subject corresponding to Lt.

X

An n by d matrix whose row vectors consist of covariate vector of additive components for each subject.

optnAdd

A list of options controls B-spline parameters for additive components, specified
by list(name=value). See ’Details’.

optnVc

A list of options controls B-spline parameters for varying-coefficient components, specified by list(name=value). See ’Details’.

Details
VCAM provides a simple algorithm based on B-spline basis to estimate its nonparametric additive
and varying-coefficient components.
Available control options for optnAdd are
nKnot A d-dimensional vector which designates the number of knots for each additive component
function estimation (default=10).
order A d-dimensional vector which designates the order of B-spline basis for each additive component function estimation (default=3).
grid A N by d matrix whose column vector consist of evaluation grid points for each component
function estimation.
and control options for optnVc are
nKnot A (d+1)-dimensional vector which designates the number of knots for overall mean function and each varying-coefficient component function estimation (default=10).
order A (d+1)-dimensional vector which designates the order of B-spline basis for overall mean
function and each varying-coefficient component function estimation (default=3).
grid A M by (d+1) matrix whose column vectors consist of evaluation grid points for overall mean
function and each varying-coefficient component function estimation.

VCAM

69

Value
A list containing the following fields:
Lt

The same with input given by Lt

LyHat

Fitted values corresponding to Ly

phiEst

An N by d matrix whose column vectors consist of esimates for each additive
component function evaluated at gridX.

beta0Est

An M-dimensional vector for overall mean function estimates evalueated at
gridT.

betaEst

An M by d matrix whose column vectors consist of esimates for each varyingcoefficient components evalueated at gridT.

gridX

The same with input given by optnAdd$grid

gridT

The same with input given by optnVc$grid

References
Zhang, X. and Wang, J.-L. (2015), "Varying-coefficient additive models for functional data", Biometrika,
Vol.102, No.1, p.15-32.
Examples
library(MASS)
set.seed(100)
n <- 100
d <- 2
Lt <- list()
Ly <- list()
m <- rep(0,2)
S <- matrix(c(1,0.5,0.5,1),nrow=2,ncol=2)
X <- pnorm(mvrnorm(n,m,S))
beta0 <- function(t) 1.5*sin(3*pi*(t+0.5))
beta1 <- function(t) 3*(1-t)^2
beta2 <- function(t) 4*t^3
phi1 <- function(x) sin(2*pi*x)
phi2 <- function(x) 4*x^3-1
for (i in 1:n) {
Ni <- sample(10:20,1)
Lt[[i]] <- sort(runif(Ni,0,1))
Ly[[i]] <- beta0(Lt[[i]]) + beta1(Lt[[i]])*phi1(X[i,1]) + beta2(Lt[[i]])*phi2(X[i,2]) + rnorm(Ni,0,0.1)
}
vcam <- VCAM(Lt,Ly,X)

70

WFDA
par(mfrow=c(1,1))
plot(unlist(vcam$LyHat),unlist(Ly),xlab='observed Y',ylab='fitted Y')
abline(coef=c(0,1),col=8)
par(mfrow=c(1,2))
plot(vcam$gridX[,1],vcam$phiEst[,1],type='l',ylim=c(-1,1),xlab='x1',ylab='phi1')
points(vcam$gridX[,1],phi1(vcam$gridX[,1]),col=2,type='l',lty=2,lwd=2)
legend('topright',c('true','est'),lwd=2,lty=c(1,2),col=c(1,2))
plot(vcam$gridX[,2],vcam$phiEst[,2],type='l',ylim=c(-1,3),xlab='x2',ylab='phi2')
points(vcam$gridX[,2],phi2(vcam$gridX[,2]),col=2,type='l',lty=2,lwd=2)
legend('topleft',c('true','est'),lwd=2,lty=c(1,2),col=c(1,2))
par(mfrow=c(1,3))
plot(vcam$gridT,vcam$beta0Est,type='l',xlab='t',ylab='beta0')
points(vcam$gridT,beta0(vcam$gridT),col=2,type='l',lty=2,lwd=2)
legend('topright',c('true','est'),lwd=2,lty=c(1,2),col=c(1,2))
plot(vcam$gridT,vcam$betaEst[,1],type='l',xlab='t',ylab='beta1')
points(vcam$gridT,beta1(vcam$gridT),col=2,type='l',lty=2,lwd=2)
legend('topright',c('true','est'),lwd=2,lty=c(1,2),col=c(1,2))
plot(vcam$gridT,vcam$betaEst[,2],type='l',xlab='t',ylab='beta2')
points(vcam$gridT,beta2(vcam$gridT),col=2,type='l',lty=2,lwd=2)
legend('topright',c('true','est'),lwd=2,lty=c(1,2),col=c(1,2))

WFDA

Warped Functional Data Analysis

Description
Pairwise curve synchronization for functional data
Usage
WFDA(Ly, Lt, optns = list())
Arguments
Ly

A list of n vectors containing the observed values for each individual.

Lt

A list of n vectors containing the observation time points for each individual
corresponding to y. Each vector should be sorted in ascending order.

optns

A list of options control parameters specified by list(name=value). See ’Details’.

Details
WFDA uses a pairwise warping method to obtain the desired alignment (registration) of the random
trajectories. The data has to be regular. The routine returns the aligned curves and the associated
warping function.
Available control options are

WFDA

71

choice Choice of estimating the warping functions (’weighted’ or ’truncated’). If ’weighted’ then
weighted averages of pairwise warping functions are computed; the weighting is based on
the inverse pairwise distances. If ’truncated’ the pairs with the top 10% largest distances
are truncated and the simple average of the remaining pairwise distances are used - default:
’truncated’
subsetProp Pairwise warping functions are determined by using a subset of the whole sample;
numeric (0,1] - default: 0.50
lambda Penalty parameter used for estimating pairwise warping functions; numeric - default :
V*10^-4, where V is the average L2 norm of y-mean(y).
nknots Number of knots used for estimating the piece-wise linear pairwise warping functions;
numeric - default: 2
isPWL Indicator if the resulting warping functions should piece-wise linear, if FALSE ’nknots’ is
ignored and the resulting warping functions are simply monotonic; logical - default: TRUE
(significantly larger computation time.)
seed Random seed for the selection of the subset of warping functions; numeric - default: 666
verbose Indicator if the progress of the pairwise warping procedure should be displayed; logical default: FALSE
Value
A list containing the following fields:
optns

Control options used.

lambda

Penalty parameter used.

aligned

Aligned curves evaluated at time ’t’

h

Warping functions for ’t’

hInv

Inverse warping functions evaluated at ’t’

costs

The mean cost associated with each curve

timing

The time required by time-warping.

References
Tang, R. and Mueller, H.G. (2008). "Pairwise curve synchronization for functional data." Biometrika
95, 875-889
Tang, R. and Mueller, H.G. (2009) "Time-synchronized clustering of gene expression trajectories."
Biostatistics 10, 32-45
Examples
N = 44;
eps = 0.123;
M = 41;
set.seed(123)
Tfinal = 3
me <- function(t) exp(-Tfinal*(((t/Tfinal^2)-0.5))^2);
T = seq(0,Tfinal,length.out = M)
recondingTimesMat = matrix(nrow = N, ncol = M)
yMat = matrix(nrow = N, ncol = M)
for (i in 1:N){

72

Wiener

}

peak = runif(min = 0.2,max = 0.8,1) * Tfinal
recondingTimesMat[i,] = sort( unique(c( seq(0.0 , peak, length.out = round((M+1)/2)),
seq( peak, Tfinal, length.out = round((M+1)/2))))) * Tfinal
yMat[i,] = me(recondingTimesMat[i,]) * rnorm(1, mean=4.0, sd= eps)
+ rnorm(M, mean=0.0, sd= eps)

Y = as.list(as.data.frame(t(yMat)))
X = rep(list(T),N)
sss = WFDA(Ly = Y, Lt = X, list( choice = 'weighted' ))
par(mfrow=c(1,2))
matplot(x= T, t(yMat), t='l', main = 'Raw', ylab = 'Y'); grid()
matplot(x= T, t(sss$aligned), t='l', main = 'Aligned', ylab='Y'); grid()

Wiener

Simulate standard Wiener processes (Brownian motions)

Description
Simulate n standard Wiener processes on [0, 1], possibly sparsifying the results.
Usage
Wiener(n = 1, pts = seq(0, 1, length = 50), sparsify = NULL,
K = 50)
Arguments
n

Sample size.

pts

A vector of points in [0, 1] specifying the support of the processes.

sparsify

A vector of integers. The number of observations per curve will be uniform
distribution on sparsify.

K

The number of components.

Details
The algorithm is based on Karhunen-Loeve expansion.
Value
If sparsify is not specified, a matrix with n rows corresponding to the samples are returned. Otherwise the sparsified sample is returned.
See Also
Sparsify

Index
BwNN, 3

kCFC, 46

CheckData, 4
CheckOptions, 4
ConvertSupport, 5
CreateBasis, 5
CreateBWPlot, 6
CreateCovPlot, 6
CreateDesignPlot, 7
CreateDiagnosticsPlot, 8
CreateFuncBoxPlot, 9
CreateModeOfVarPlot, 10
CreateOutliersPlot, 11
CreatePathPlot, 12
CreateScreePlot, 13
CreateStringingPlot, 14
cumtrapzRcpp, 14

Lwls1D, 47
Lwls2D, 47
Lwls2DDeriv, 48
MakeBWtoZscore02y, 49
MakeFPCAInputs, 50
MakeGPFunctionalData, 50
MakeHCtoZscore02y, 51
MakeLNtoZscore02y, 51
MakeSparseGP, 52
medfly25, 52
MultiFAM, 53
NormCurvToArea, 57
plot.FPCA (CreateDiagnosticsPlot), 8
predict.FPCA, 57
print.FPCA, 58
print.FSVD, 59
print.WFDA, 59

Dyn_test, 16
DynCorr, 15
FAM, 17
FCCor, 20
FClust, 21
FCReg, 23
fdapace, 24
fdapace-package (fdapace), 24
fitted.FPCA, 25
fitted.FPCAder, 27
FOptDes, 28
FPCA, 29
FPCAder, 32
FPCquantile, 34
FPCReg, 35
FSVD, 38
FVPA, 40

SBFitting, 59
SelectK, 62
SetOptions, 63
Sparsify, 63
str.FPCA, 64
Stringing, 64
trapzRcpp, 65
TVAM, 66
VCAM, 68
WFDA, 70
Wiener, 72

GetCrCorYX, 41
GetCrCorYZ, 42
GetCrCovYX, 42
GetCrCovYZ, 43
GetNormalisedSample, 45
GetNormalizedSample
(GetNormalisedSample), 45
73



Source Exif Data:
File Type                       : PDF
File Type Extension             : pdf
MIME Type                       : application/pdf
PDF Version                     : 1.5
Linearized                      : No
Page Count                      : 73
Page Mode                       : UseOutlines
Author                          : 
Title                           : 
Subject                         : 
Creator                         : LaTeX with hyperref package
Producer                        : pdfTeX-1.40.19
Create Date                     : 2019:04:24 22:13:18-07:00
Modify Date                     : 2019:04:24 22:13:18-07:00
Trapped                         : False
PTEX Fullbanner                 : This is pdfTeX, Version 3.14159265-2.6-1.40.19 (TeX Live 2018) kpathsea version 6.3.0
EXIF Metadata provided by EXIF.tools

Navigation menu