Package 'fdapace'

Description

Input a list of time points Lt, and the number of unique neighbors k. Obtain the minimum bandwidth guaranteeing k unique neighbours.

Usage

BwNN(Lt, k = 3, onlyMean = FALSE, onlyCov = FALSE)

Arguments

Examples

tinyGrid = list(c(1,7), c(2,3),  6,  c(2,4), c(4,5))
BwNN(tinyGrid, k = 2) # c(3,2)

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

Description

Check if the options structure is valid and set the NULL options

Usage

CheckOptions(t, optns, n)

Arguments

Description

Convert the support of a given function 1-D or 2-D function from fromGrid to toGrid. Both grids need to be sorted. This is an interpolation/convenience function.

Usage

ConvertSupport(
  fromGrid,
  toGrid,
  mu = NULL,
  Cov = NULL,
  phi = NULL,
  isCrossCov = FALSE
)

Arguments

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

Value

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

Examples

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

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 variation for 50

Usage

CreateBWPlot(fpcaObj, derOptns = NULL, bwMultipliers = NULL)

Arguments

Examples

set.seed(1)
n <- 40
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)

Description

This function will open a new device if not instructed otherwise.

Usage

CreateCovPlot(
  fpcaObj,
  covPlotType = "Fitted",
  corr = FALSE,
  isInteractive = FALSE,
  colSpectrum = NULL,
  ...
)

Arguments

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) ##plotting the covariance surface
CreateCovPlot(res, corr = TRUE) ##plotting the correlation surface

Description

Create design plots for functional data. See Yao, F., Müller, H.G., Wang, J.L. (2005). Functional data analysis for sparse longitudinal data. J. American Statistical Association 100, 577-590 for interpretation and usage of these plots. This function will open a new device as default.

Usage

CreateDesignPlot(
  Lt,
  obsGrid = NULL,
  isColorPlot = TRUE,
  noDiagonal = TRUE,
  addLegend = TRUE,
  ...
)

Arguments

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))))

Description

Deprecated. Use plot.FPCA instead.

Plotting the results of an FPCA, including printing the design plot, mean function, scree-plot and the first three eigenfunctions for a functional sample. If provided with a derivative options object (?FPCAder), it will return the differentiated mean function and first two principal modes of variation 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

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.

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)

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 gray area to the area spanned by the curves within the 25th and 75-th percentile and the light gray to the area spanned by the curves within the 2.5th and 97.5-th percentile.

Usage

CreateFuncBoxPlot(fpcaObj, optns = list(), ...)

Arguments

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, relevant only for variant='KDE')

addIndx

vector of indices 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

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)) )

Description

Creates 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 variation around the mean: $\pm Q \sqrt{\lambda_k} \phi_k$ for the dark grey area Q = 1, and for the light grey are Q = 2. In the case of 'rainbowPlot' the blue edge corresponds to Q = -3, the green edge to Q = +3 and the red-line to Q = 0 (the mean).

Usage

CreateModeOfVarPlot(
  fpcaObj,
  k = 1,
  rainbowPlot = FALSE,
  colSpectrum = NULL,
  ...
)

Arguments

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)

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

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, relevant 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, relevant 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'), relevant only for groupingType='slice')

groupingType

string specifying if a slice grouping ('slice') or a standard percentile/bagplot grouping ('standard') should be returned (default: 'standard')

fIndices

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.

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

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)

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

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:

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)

Description

This function will open a new device if not instructed otherwise.

Usage

CreateScreePlot(fpcaObj, ...)

Arguments

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)

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

Details

This approach is based on Chen, K., Chen, K., Müller, H.G., Wang, J.L. (2011). Stringing high-dimensional data for functional analysis. J. American Statistical Association 106, 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$StringingOrder] - p:1))
if(diff_rev <= diff_norev){
  stringingfit$StringingOrder = rev(stringingfit$StringingOrder)
  stringingfit$Ly = lapply(stringingfit$Ly, rev)
}
CreateStringingPlot(stringingfit, 1:20)

Description

Cumulative Trapezoid Rule Numerical Integration using Rcpp

Usage

cumtrapzRcpp(X, Y)

Arguments

Description

Perform one sample (H0: Dynamic correlation = 0) or two sample (H0:Dynamic_correlation_1 = Dynamic_correlation_2) bootstrap test of H_0: Dynamical Correlation=0.

Usage

Dyn_test(x1, y1, t1, x2, y2, t2, B = 1000)

Arguments

Value

a list of the following

References

Dubin J A, Müller H G. (2005) Dynamical correlation for multivariate longitudinal data. Journal of the American Statistical Association 100(471): 872-881.

Liu S, Zhou Y, Palumbo R, Wang, J.L. (2016). Dynamical correlation: A new method for quantifying synchrony with multivariate intensive longitudinal data. Psychological Methods 21(3): 291.

Examples

n=20             # 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)   # 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=500) # using B=500 for speed consideration

Description

Calculates the Dynamical Correlation for 2 paired dense regular functional data observed on the same grid.

Usage

DynCorr(x, y, t)

Arguments

Value

A length n vector of individual dynamic correlations. The dynamic correlation can be obtained by taking average of this vector.

References

Dubin J A, Müller H G. Dynamical correlation for multivariate longitudinal data (2005). Journal of the American Statistical Association 100(471): 872-881. Liu S, Zhou Y, Palumbo R, Wang, J.L. (2016). Dynamical correlation: A new method for quantifying synchrony with multivariate intensive longitudinal data. Psychological methods 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)

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

Details

FAM fits functional additive models for a scalar response and single predictor process proposed by Müller and Yao (2007) that

$E(Y | \mathbf{X}) = \sum_{k=1}^K g_{k}(\xi_{k}),$

where $\xi_{k}$ stand for the k-th FPC score of the the predictor process.

Value

A list containing the following fields:

References

Müller, H.-G. and Yao, F. (2005), "Functional additive models", JASA, Vol.103, No.484, p.1534-1544.

Examples

set.seed(1000)

library(MASS)

f1 <- function(t) 0.5*t
f2 <- function(t) 2*cos(2*pi*t/4)
f3 <- function(t) 1.5*sin(2*pi*t/4)
f4 <- function(t) 2*atan(2*pi*t/4)

n <- 100
N <- 100

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 <- function(t) sqrt(2)*sin(2*pi*t)
phi2 <- function(t) sqrt(2)*sin(4*pi*t)
phi3 <- function(t) sqrt(2)*cos(2*pi*t)
phi4 <- function(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)
}

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

op <- 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')
par(op)

# fitting
fit <- FAM(Y=Y,Lx=Lx,Lt=Lt,nEval=0)
yHat <- fit$mu+apply(fit$fam,1,'sum')
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')
plot(yHat,YTest,xlim=c(-10,10))
abline(coef=c(0,1),col=2)

Description

Calculation of functional correlation between two simultaneously observed processes.

Usage

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

Arguments

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:

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)

Description

Default: Cluster functional data using the functional principal component (FPC) scores obtained from the data 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

Details

Within EMCluster, uses the model initiated EMCluster::em.EM and returns the optimal model based on EMCluster::emcluster. See ?EMCluster::emcluster for details.

Value

A list containing the following fields:

References

Wei-Chen Chen and Ranjan Maitra, "EMCluster: EM Algorithm for Model-Based Clustering 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

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%

Description

Functional concurrent regression with dense or sparse functional data for scalar or functional dependent variables. Note: function-to-scalar regression can also be handled using the VCAM function in fdapace.

Usage

FCReg(
  vars,
  userBwMu,
  userBwCov,
  outGrid,
  kern = "gauss",
  measurementError = TRUE,
  diag1D = "none",
  useGAM = FALSE,
  returnCov = TRUE
)

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. WARNING! For very sparse functional data, setting measurementError = TRUE is not recommended.

Value

A list containing the following fields:

References

Yao, F., Müller, H.G., Wang, J.L. "Functional Linear Regression Analysis for Longitudinal Data." Annals of Statistics 33, (2005): 2873-2903.(Dense data) Sentürk, D., Müller, H.G. "Functional varying coefficient models for longitudinal data." J. American Statistical Association, 10, (2010): 1256–1264. Sentürk, 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)

Description

fdapace 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. This core algorithm yields covariance and mean functions, eigenfunctions and principal component (scores), for both functional data and derivatives, for both dense (functional) and sparse (longitudinal) sampling designs. For sparse designs, it provides fitted continuous trajectories with confidence bands, even for subjects with very few longitudinal observations. PACE is a viable and flexible alternative to random effects modeling of longitudinal data. There is also a Matlab version (PACE) that contains some methods not available on fdapace and vice versa.

Details

Links for fdapace/PACE: Matlab version of pace at http://anson.ucdavis.edu/~mueller/data/pace.html Papers and background at http://anson.ucdavis.edu/~mueller/ and http://www.stat.ucdavis.edu/~wang/

PACE is based on the idea that observed functional data are generated by a sample of underlying (but usually not fully observed) random trajectories that are realizations of a stochastic process. It does not rely on pre-smoothing of trajectories, which is problematic if functional data are sparsely sampled.

The functional principal components can be used for further statistical analysis depending on the demands of a user, for example if one has densely sampled functional predictors and a generalized response, such as in a GLM, the predictor functions can be replaced by their first couple of principal component scores that will then be used as predictors; one can also easily fit polynomial functional models by using powers (usually squares) and interactions of functional principal components among the predictors for a scalar response.

fdapace is a comprehensive package that directly implements fitting of the following models:

• functional linear regression

• functional additive regression

• functional covariance and correlation (via dynamic correlation)

• functional clustering

• concurrent (varying coefficient) regression models for sparse and dense designs

• varying coefficient additive models

• multivariate functional data analysis (normalization and functional singular component analysis)

• variance processes and volatility processes (the latter of interest in finance)

• optimal designs for longitudinal data analysis (for trajectory prediction and for functional linear regression)

• stringing, a method to convert high-dimensional data into functional data

• quantile regression, with functions as predictors

Author(s)

Maintainer: Yidong Zhou [email protected] (ORCID)

Authors:

• Han Chen

• Su I Iao

• Poorbita Kundu

• Hang Zhou

• Satarupa Bhattacharjee

• Cody Carroll (ORCID)

• Yaqing Chen

• Xiongtao Dai

• Jianing Fan

• Alvaro Gajardo

• Pantelis Z. Hadjipantelis

• Kyunghee Han

• Hao Ji

• Changbo Zhu

• Hans-Georg Müller [email protected] [copyright holder, thesis advisor]

• Jane-Ling Wang [email protected] [copyright holder, thesis advisor]

Other contributors:

• Paromita Dubey [contributor]

• Shu-Chin Lin [contributor]

References

• Wang, J.L., Chiou, J., Müller, H.G. (2016). Functional data analysis. Annual Review of Statistics and Its Application 3, 257–295;

• Chen, K., Zhang, X., Petersen, A., Müller, H.G. (2017). Quantifying infinite-dimensional data: Functional Data Analysis in action. Statistics in Biosciences 9, 582–-604.

See Also

Useful links:

• https://github.com/functionaldata/tPACE

• Report bugs at https://github.com/functionaldata/tPACE/issues

_PACKAGE

Description

Combines 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. For details consult the example.

Usage

## S3 method for class 'FPCA'
fitted(
  object,
  K = NULL,
  derOptns = list(p = 0),
  ciOptns = list(alpha = NULL, cvgMethod = NULL),
  ...
)

Arguments

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 Müller (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 between 'interval' (pointwise coverage) and 'band' (simultaneous coverage). 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:

References

Yao, F., Müller, 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 Müller. "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

op <- par(mfrow=c(2,3))
ind <- sample(1:n,6)
for (i in 1:6) {
 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]])
}
par(op)

Description

Combines the zero-meaned fitted values and the mean derivative to get the fitted values for the derivative 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

Value

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

References

Liu, Bitao, and Hans-Georg Müller. "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)

Description

Functional linear models for scalar or functional responses and scalar and/or functional predictors.

Usage

FLM(Y, X, XTest = NULL, optnsListY = NULL, optnsListX = NULL, nPerm = NULL)

Arguments

Value

A list of the following:

References

Yao, F., Müller, H.G., Wang, J.L. (2005). Functional linear regression analysis for longitudinal data. Annals of Statistics 33, 2873–2903. Hall, P., Horowitz, J.L. (2007). Methodology and convergence rates for functional linear regression. The Annals of Statistics, 35(1), 70–91.

Examples

set.seed(1000)

library(MASS)

### functional covariate
phi1 <- function(t,k) sqrt(2)*sin(2*pi*k*t)
phi2 <- function(t,k) sqrt(2)*cos(2*pi*k*t)

lambdaX <- c(1,0.7)

# training set
n <- 50
Xi <- matrix(rnorm(2*n),nrow=n,ncol=2)

denseLt <- list(); denseLy <- list()
sparseLt <- list(); sparseLy <- list()

t0 <- seq(0,1,length.out=51)
for (i in 1:n) {
 denseLt[[i]] <- t0
  denseLy[[i]] <- lambdaX[1]*Xi[i,1]*phi1(t0,1) + lambdaX[2]*Xi[i,2]*phi1(t0,2)
  
  ind <- sort(sample(1:length(t0),3))
  sparseLt[[i]] <- t0[ind]
  sparseLy[[i]] <- denseLy[[i]][ind]
}

denseX <- list(Ly=denseLy,Lt=denseLt)
sparseX <- list(Ly=sparseLy,Lt=sparseLt)

denseX <- list(X=denseX)
sparseX <- list(X=sparseX)

# test set
N <- 30

XiTest <- matrix(rnorm(2*N),nrow=N,ncol=2)

denseLtTest <- list(); denseLyTest <- list()

sparseLtTest <- list(); sparseLyTest <- list()

t0 <- seq(0,1,length.out=51)
for (i in 1:N) {
  denseLtTest[[i]] <- t0
  denseLyTest[[i]] <- lambdaX[1]*XiTest[i,1]*phi1(t0,1) + lambdaX[2]*XiTest[i,2]*phi1(t0,2)
  
  ind <- sort(sample(1:length(t0),5))
  sparseLtTest[[i]] <- t0[ind]
  sparseLyTest[[i]] <- denseLyTest[[i]][ind]
}

denseXTest <- list(Ly=denseLyTest,Lt=denseLtTest)
sparseXTest <- list(Ly=sparseLyTest,Lt=sparseLtTest)

denseXTest <- list(X=denseXTest)
sparseXTest <- list(X=sparseXTest)


### scalar response
beta <- c(1, -1)
Y <- c(Xi%*%diag(lambdaX)%*%beta) + rnorm(n,0,0.5)
YTest <- c(XiTest%*%diag(lambdaX)%*%beta) + rnorm(N,0,0.5)

## dense
denseFLM <- FLM(Y=Y,X=denseX,XTest=denseXTest,optnsListX=list(FVEthreshold=0.95))

trueBetaList <- list()
trueBetaList[[1]] <- cbind(phi1(denseFLM$workGridX[[1]],1),phi1(denseFLM$workGridX[[1]],2))%*%beta

# coefficient function estimation error (L2-norm)
plot(denseFLM$workGridX[[1]],denseFLM$betaList[[1]],type='l',xlab='t',ylab=paste('beta',1,sep=''))
points(denseFLM$workGridX[[1]],trueBetaList[[1]],type='l',col=2)

denseEstErr <-
  sqrt(trapzRcpp(denseFLM$workGridX[[1]],(denseFLM$betaList[[1]] - trueBetaList[[1]])^2))
denseEstErr

op <- par(mfrow=c(1,2))
plot(denseFLM$yHat,Y,xlab='fitted Y', ylab='observed Y')
abline(coef=c(0,1),col=8)
plot(denseFLM$yPred,YTest,xlab='predicted Y', ylab='observed Y')
abline(coef=c(0,1),col=8)
par(op)

# prediction error
densePredErr <- sqrt(mean((YTest - denseFLM$yPred)^2))
densePredErr

Description

Bootstrap pointwise confidence intervals for the coefficient functions in functional linear models.

Usage

FLMCI(Y, X, level = 0.95, R = 999, optnsListY = NULL, optnsListX = NULL)

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. WARNING! For very sparse functional data, setting measurementError=TRUE is not recommended.

Value

A list containing the following fields:

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

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:

References

Ji, H., Müller, H.G. (2017) "Optimal Designs for Longitudinal and Functional Data" Journal of the Royal Statistical Society: Series B 79, 859-876.

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))

Description

FPCA for dense or sparse functional data.

Usage

FPCA(Ly, Lt, optns = list())

Arguments

Details

If the input is sparse data, make sure you check the design plot is dense and the 2D domain is well covered by support points, using plot or CreateDesignPlot. Some study design such as snippet data (where each subject is observed only on a sub-interval of the period of study) will have an ill-covered design plot, in which case the nonparametric covariance estimate will be unreliable. WARNING! Slow computation times may occur if the dataType argument is incorrect. If FPCA is taking a while, please double check that a dense design is not mistakenly coded as 'Sparse'. Applying FPCA to a mixture of very dense and sparse curves may result in computational issues.

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'

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 covariance function; numeric (0,1] - default: 0.99

FVEfittedCov

Fraction-of-Variance explained by the components that are used to construct fittedCov; numeric (0,1] - default: NULL (all components available will be used)

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-2,nRegGrid-2), N:# of curves, nRegGrid:# of support points in each direction of covariance surface

methodXi

The method to estimate the PC scores; 'CE' (Conditional Expectation), 'IN' (Numerical Integration) - default: 'CE' for sparse data and dense data with missing values, 'IN' for dense data. If time points are irregular but spacing is small enough, 'IN' method is utilized by default.

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 percentages of the time range to be considered as left and right boundary regions of the time window of observation - default (0,1) which corresponds to no boundary

methodRho

The method of regularization (add to diagonal of covariance surface) in estimating principal component scores; 'trunc': rho is truncation of sigma2, 'ridge': rho is a ridge parameter, 'vanilla': vanilla approach - default "vanilla".

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

usergrid

Whether to use observation grid for fitting, if false will use equidistant grid. logical - default:FALSE

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

userSigma2

The user-defined measurement error variance. A positive scalar. If specified then the vanilla approach is used (methodRho is set to 'vanilla', 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

imputeScores

Whether to impute the FPC scores or not; default: 'TRUE'

verbose

Display diagnostic messages; logical - default: FALSE

Value

A list containing the following fields:

References

Yao, F., Müller, 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, Müller, Hans-Georg and Wang, Jane-Ling (2005). "Functional data analysis for sparse longitudinal data." Journal of the American Statistical Association 100, no. 470 577-590. (Sparse data FPCA)

Liu, Bitao and Müller, Hans-Georg (2009). "Estimating derivatives for samples of sparsely observed functions, with application to online auction dynamics." Journal of the American Statistical Association 104, no. 486 704-717. (Sparse data FPCA)

Castro, P. E., Lawton, W.H. and Sylvestre, E.A. (1986). "Principal modes of variation for processes with continuous sample curves." Technometrics 28, no. 4, 329-337. (modes of variation for 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')
CreateCovPlot(res, corr = TRUE)

Description

Obtain the derivatives of eigenfunctions/ eigenfunctions of derivatives (note: these two are not the same)

Usage

FPCAder(fpcaObj, derOptns = list(p = 1))

Arguments

Details

Available derivative options are

method

The method used for obtaining the derivatives – default is 'FPC', which is the derivatives of eigenfunctions; 'DPC': eigenfunctions of derivatives, with G^(1,1) estimated by an initial kernel local smoothing step for G^(1,0), then applying a 1D smoother in the second direction; 'FPC': functional principal component, based on smoothing the eigenfunctions; 'FPC1': functional principal component, based on smoothing G^(1,0). The latter may produce better estimates than 'FPC' but is slower.

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, where S is the length of the domain).

kernelType

Smoothing kernel choice; same available types are FPCA(). default('epan')

References

Dai, X., Tao, W., Müller, H.G. (2018). Derivative principal components for representing the time dynamics of longitudinal and functional data. Statistica Sinica 28, 1583–1609. (DPC) Liu, Bitao, and Hans-Georg Müller. "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 <- 50 
M <- 30
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)

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

Value

A list of the following

References

Chen, K., Müller, H.G. (2011). Conditional quantile analysis when covariates are functions, with application to growth data. J. Royal Statistical Society B 74, 67-89

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)

Description

FSVD for a pair of dense or sparse functional data.

Usage

FSVD(
  Ly1,
  Lt1,
  Ly2,
  Lt2,
  FPCAoptns1 = NULL,
  FPCAoptns2 = NULL,
  SVDoptns = list()
)

Arguments

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 composite 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:

References

Yang, W., Müller, H.G., Stadtmüller, U. (2011). Functional singular component analysis. J. Royal Statistical Society B73, 303-324.

Description

Functional Variance Process Analysis for dense functional data

Usage

FVPA(y, t, q = 0.1, optns = list(error = TRUE, FVEthreshold = 0.9))

Arguments

Value

A list containing the following fields:

References

Hans-Georg Müller, Ulrich Stadtmüller 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)

Description

Covariance surface estimation for dense or sparse functional data.

Usage

GetCovSurface(Ly, Lt, optns = list())

Arguments

Value

A list containing the following fields:

Examples

set.seed(1)
n <- 20
pts <- seq(0, 1, by=0.025)
sampWiener <- Wiener(n, pts)
mu = sin(2*pi*pts)
sampWiener <- Sparsify(t(t(sampWiener) + mu), pts, 10)
res = GetCovSurface(Ly = sampWiener$Ly, Lt = sampWiener$Lt)

Description

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

Usage

GetCrCorYX(ccXY, ccXX, ccYY)

Arguments

Value

A cross-correlation matrix between variables X and Y.

Description

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

Usage

GetCrCorYZ(ccYZ, acYY, covZ)

Arguments

Value

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

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

Value

A list containing:

References

Yang, Wenjing, Hans-Georg Müller, and Ulrich Stadtmüller. "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)

Description

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

Usage

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

Arguments

Value

A list containing:

References

Yang, Wenjing, Hans-Georg Müller, and Ulrich Stadtmüller. "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))

Description

Note that bootstrap pointwise confidence intervals do not work for sparsely observed data.

Usage

GetMeanCI(Ly, Lt, level = 0.95, R = 999, optns = list())

Arguments

Value

A list of two elements:

.

Examples

n <- 30
tgrid <- seq(0,1,length.out=21)
phi1 <- function(t) sqrt(2)*sin(2*pi*t)
phi2 <- function(t) sqrt(2)*sin(4*pi*t)
Lt <- rep(list(tgrid), n)
Ly <- lapply(1:n, function(i){
 tgrid + rnorm(1,0,2) * phi1(tgrid) + rnorm(1,0,0.5) * phi2(tgrid) + rnorm(1,0,0.01)
 })
res <- GetMeanCI(Lt = Lt, Ly = Ly, level = 0.9)

Description

Mean curve calculation for dense or sparse functional data.

Usage

GetMeanCurve(Ly, Lt, optns = list())

Arguments

Value

A list containing the following fields:

Examples

set.seed(1)
n <- 20
pts <- seq(0, 1, by=0.025)
sampWiener <- Wiener(n, pts)
mu = sin(2*pi*pts)
sampWiener <- Sparsify(t(t(sampWiener) + mu), pts, 10)
res = GetMeanCurve(Ly = sampWiener$Ly, Lt = sampWiener$Lt, optns = list(plot = TRUE))

Description

Normalise sparse functional sample given in an FPCA object

Usage

GetNormalisedSample(fpcaObj, errorSigma = FALSE)

GetNormalizedSample(...)

Arguments

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)

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

Value

A list containing the following fields:

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

Examples

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

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

Value

Vector of length M with measurement values at the the point specified by 'xout'

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")
)

Arguments

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.

Description

Two dimensional local weighted least squares smoother. Only a local linear smoother for estimating the original curve is available (no higher order)

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