cv.sigcom: Cross-validation for linear function-on-function regression
In FRegSigCom: Functional Regression using Signal Compression Approach

Description Usage Arguments Details Value Author(s) References Examples

This function is used to perform cross-validation and build the final model using the signal compression approach for the following linear function-on-function regression model with or without scalar predictors:

Y(t)= μ(t)+Z'α(t)+\int X_1(s)β_1(s,t)ds+...+\int X_p(s)β_p(s,t)ds+ε(t),

where μ(t) is the intercept function, Z is a multivariate predictor and α(t) is the vector of corresponding coefficient functions. When Z=NULL, there is no scalar predictor. The {X_i(s),1≤ i≤ p} are p functional predictors and {β_i(s,t),1≤ i≤ p} are their corresponding coefficient functions, where p is a positive integer. When p is large (e.g., greater than 6), in addition to smoothness penalty, one may want to consider to impose sparsity penalty (see cv.hd()). The ε(t) is the noise function.

We require that all the sample curves of each functional predictor are observed in a common dense grid of time points, but the grid can be different for different predictors. All the sample curves of the functional response are observed in a common dense grid.

1
2
3

cv.sigcom(X, Y, t.x, t.y, Z = NULL, s.n.basis = 50, t.n.basis = 50,
           K.cv = 5, upper.comp = 20, thresh = 0.001,
           basis.type.x="Bspline", basis.type.y="Bspline")

`X`	a list of length p, the number of functional predictors. Its i-th element is the nm_i* data matrix for the i-th functional predictor X_i(s), where n is the sample size and m_i is the number of observation time points for X_i(s).
`Y`	the nm* data matrix for the functional response Y(t), where n is the sample size and m is the number of the observation time points for Y(t).
`t.x`	a list of length p. Its i-th element is the vector of obesrvation time points of the i-th functional predictor X_i(s), 1≤ i≤ p.
`t.y`	the vector of obesrvation time points of the functional response Y(t).
`Z`	the nq* data matrix for multivariate scalar predictors, where n is the sample size and q is the number of the scalar predictors. Default is NULL, indicating no scalar predictors.
`s.n.basis`	the number of B-spline basis functions used for estimating the functions ψ_{ik}(s)'s. Default is 50.
`t.n.basis`	the number of B-spline basis functions used for estimating the functions w_{k}(t)'s. Default is 50.
`K.cv`	the number of CV folds. Default is 5.
`upper.comp`	the upper bound for the maximum number of components to be calculated. Default is 20.
`thresh`	a number between 0 and 1 used to determine the maximum number of components we need to calculate. The maximum number is between one and the "upp.comp" above. The optimal number of components will be chosen between 1 and this maximum number, together with other tuning parameters by cross-validation. A smaller thresh value leads to a larger maximum number of components and a longer running time. A larger thresh value needs less running time, but may miss some important components and lead to a larger prediction error. Default is 0.001.
`basis.type.x`	the type of basis functions ψ_{ik}(s)'s . Only "BSpline" (default) and "Fourier" are supported.
`basis.type.y`	the type of basis functions w_{k}(t)'s. Only "BSpline" (default) and "Fourier" are supported.

We consider the model with p=1 and Z=NULL. Details for the model with scalar and multiple functional predictors can be found in the reference and its supplementary material. We consider a special decomposition of the coefficient function:

β(s,t)=ψ_{1}(s)w_1(t)+...+ψ_{K}(s)w_K(t).

Let \{(x_1(s), y_1(t)),...,(x_n(s),y_n(t))\} denote the samples. We first estimate ψ_{k}(s) for each k>0 by solving a penalized generalized functional eigenvalue problem

max_{ψ} \frac{\int\int ψ(s)\hat{B}(s,s')ψ(s')dsds'}{\int\int ψ(s)\hat{Σ}(s,s')ψ(s')dsds'+P(ψ)}

{\rm{s.t.}}\quad \int\int ψ(s) \hat{Σ}(s,s')ψ(s')dsds'=1

{\rm{and}}\quad \int\int ψ(s) \hat{Σ}(s,s')ψ_{k'}(s')dsds'=0 \quad {\rm{for}} \quad k'<k

where \hat{B}(s,s')=∑_{\ell=1}^n∑_{\ell'=1}^n \{x_{\ell}(s)-\bar{x}(s)\}\int\{y_{\ell}(t)-\bar{y}(t)\}\{y_{\ell'}(t)-\bar{y}(t)\}dt \{x_{\ell'}(s')-\bar{x}(s')\}/n^2, \hat{Σ}(s,s')=∑_{\ell=1}^n \{x_{\ell}(s)-\bar{x}(s)\}\{x_{\ell}(s')-\bar{x}(s')\}/n, and penalty

λ{|| ψ_{k}||^2} + λτ{|| ψ''_{k}||^2}.

Then we estimate {w_{k}(t), k>0} by regressing \{y_{\ell}(t)\}_{\ell=1}^n on \{\hat{z}_{\ell 1},... \hat{z}_{\ell K}\}_{\ell=1}^n with penalty κ ∑_{k=0}^K \|w''_k\|^2 tuned by the smoothness parameter κ. Here \hat{z}_{\ell k}= \int (x_{\ell}(s)-\bar{x}(s))\hat{ψ}_{k}(s)ds.

An object of the “cv.sigcom” class, which is used in the function pred.sigcom for prediction and getcoef.sigcom for extracting the estimated coefficient functions.

`t.x.list`	the input t.x.
`t.y`	the input t.y.
`Z`	the input Z
`opt.index`	index of the optimal lambda.
`opt.lambda`	optimal value for lambda, the parameter tuning the penalty on {ψ_{ik}(s), 1≤ i ≤ p}.
`opt.tau`	optimal value for tau, the smoothness parameter for {ψ_{ik}(s), 1≤ i ≤ p}.
`opt.kappa`	optimal value for kappa, the smoothness parameter for w_k(t).
`opt.K`	optimal number of components to be selected.
`min.error`	minimum CV error.
`errors`	list for CV errors.
`x.smooth.params`	a list for internal use.
`y.smooth.params`	a list for internal use.
`fit.1`	a list for internal use.
`is.null.Z`	a logic value indicating whether Z is null.

Ruiyan Luo and Xin Qi

Ruiyan Luo and Xin Qi (2017) Function-on-Function Linear Regression by Signal Compression, Journal of the American Statistical Association. 112(518), 690-705. https://doi.org/10.1080/01621459.2016.1164053

#########################################################################
# Example 1: linear function-on-function regresion with one predictor
#             curve and two scalar predictors
#########################################################################
ptm <- proc.time()
library(FRegSigCom)
library(refund)
data(DTI)
tmp=1*(DTI$sex=="male")
Z.0=cbind(DTI$case, tmp)

I=which(is.na(apply(DTI$cca,1,mean)))
Y=DTI$cca[-I,] # functional response
X=list(DTI$rcst[-I,-(1:12)]) #functional predictor
Z=Z.0[-I,] # scalar predictor

t.x <- list(seq(0,1,length=dim(X[[1]])[2]))
t.y <- seq(0,1,length=dim(Y)[2])
# randomly split all the observations into a training set with 200 observations
# and a test set.
train.id=sample(1:nrow(Y), 100)
X.train <- list(X[[1]][train.id,])
Y.train <- Y[train.id, ]
Z.train <- Z[train.id, ]
X.test <- list(X[[1]][-(train.id),])
Y.test <- Y[-(train.id), ]
Z.test <- Z[-(train.id), ]
fit.cv=cv.sigcom(X.train, Y.train, t.x, t.y, Z=Z.train)
Y.pred=pred.sigcom(fit.cv, X.test, Z.test=Z.test)
error<- mean((Y.pred-Y.test)^2)
print(c(" prediction error=", error))
coef.est.list=getcoef.sigcom(fit.cv)
mu.est=coef.est.list[[1]] # intercept curve
beta.est=coef.est.list[[2]] # coefficient functions of functional predictors
gamma.est=coef.est.list[[3]] # coefficient functions of scalar predictors
print(proc.time()-ptm)



#########################################################################
# Example 2: linear function-on-function regresion with four predictor curves
#########################################################################


ptm <- proc.time()
library(FRegSigCom)
data(ocean)
Y=ocean[[1]]
Y.train=Y[1:80,]
Y.test=Y[-(1:80),]
t.y=seq(0,1, length.out = ncol(Y))
X.list=list()
X.train.list=list()
X.test.list=list()
t.x.list=list()
for(i in 1:4)
{
  X.list[[i]]=ocean[[i+1]]
  X.train.list[[i]]=X.list[[i]][1:80,]
  X.test.list[[i]]=X.list[[i]][-(1:80),]
  t.x.list[[i]]=seq(0,1, length.out = ncol(X.list[[i]]))
}
fit.cv=cv.sigcom(X.train.list, Y.train, t.x.list, t.y)
Y.pred=pred.sigcom(fit.cv, X.test.list)
error<- mean((Y.pred-Y.test)^2)
print(c(" prediction error=", error))
coef.list=getcoef.sigcom(fit.cv)
mu.est=coef.list[[1]] # intercept curve
beta.est=coef.list[[2]] # coefficient functions of functional predictors
print(proc.time()-ptm)