cv.sigcom: Cross-validation for linear function-on-function regression

Description

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)= \mu(t)+Z^T\alpha(t)+\sum_{i=1}^p\int_{a_i}^{b_i}X_i(s)\beta_i(s,t)ds+\varepsilon(t)$$ where $\mu(t)$ is the intercept function, $Z$ is a multivariate predictor and $\alpha(t)$ is the vector of corresponding coefficient functions. When $Z=$NULL, there is no scalar predictor. The $\{X_i(s),1\le i\le p\}$ are $p$ functional predictors and $\{\beta_i(s,t),1\le i\le 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 $\epsilon(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.

Usage

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

Arguments

a list of length $p$, the number of functional predictors. Its $i$-th element is the $n\times m_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)$.

the $n\times m$ 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\le i\le p$.

t.y

the vector of obesrvation time points of the functional response $Y(t)$.

the $n\times q$ 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 $\psi_{ik}(s)$. Default is 50.

t.n.basis

the number of B-spline basis functions used for estimating the functions $w_{k}(t)$. 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 $\psi_{ik}(s)$ . Only "BSpline" (default) and "Fourier" are supported.

basis.type.y

the type of basis functions $w_{k}(t)$. Only "BSpline" (default) and "Fourier" are supported.

Value

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.

the input Z

opt.index

index of the optimal $lambda$.

opt.lambda

optimal value for $lambda$, the parameter tuning the penalty on ${\psi_{ik}(s), 1\le i \le p}$.

opt.tau

optimal value for $tau$, the smoothness parameter for ${\psi_{ik}(s), 1\le i \le 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.

Details

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: $$\beta(s,t)=\sum_{k=1}^\infty\psi_{k}(s)w_k(t).$$ Let $\{(x_1(s), y_1(t)),...,(x_n(s),y_n(t))\}$ denote the samples. We first estimate $\psi_{k}(s)$ for each $k>0$ by solving a penalized generalized functional eigenvalue problem $$max_{\psi} \frac{\int\int \psi(s)\hat{B}(s,s')\psi(s')dsds'}{\int\int \psi(s)\hat{\Sigma}(s,s')\psi(s')dsds'+P(\psi)}$$ $$ {\rm{s.t.}}\quad \int\int \psi(s) \hat{\Sigma}(s,s')\psi(s')dsds'=1$$ $$ {\rm{and}}\quad \int\int \psi(s) \hat{\Sigma}(s,s')\psi_{k'}(s')dsds'=0 \quad {\rm{for}} \quad k'<k$$ where $\hat{B}(s,s')=\sum_{\ell=1}^n\sum_{\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{\Sigma}(s,s')=\sum_{\ell=1}^n \{x_{\ell}(s)-\bar{x}(s)\}\{x_{\ell}(s')-\bar{x}(s')\}/n$, and penalty $$P(\psi)=\lambda\{|| \psi_{k}||^2+ \tau ||\psi''_{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 $\kappa \sum_{k=0}^K \|w''_k\|^2$ tuned by the smoothness parameter $\kappa$. Here $\hat{z}_{\ell k}= \int (x_{\ell}(s)-\bar{x}(s))\hat{\psi}_{k}(s)ds$.

References

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

Examples

Run this code

# NOT RUN {
#########################################################################
# 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)
# }

Run the code above in your browser using DataLab