# A partial least squares approach for function-on-function interaction
regression

Ufuk Beyaztas111Postal address: Department of Economics and Finance, Pirireis
University, Istanbul, Turkey; Email<EMAIL_ADDRESS>ORCID:
0000-0002-5208-4950
Department of Economics and Finance
Piri Reis University

Han Lin Shang
Department of Actuarial Studies and Business Analytics
Macquarie University

###### Abstract

A partial least squares regression is proposed for estimating the function-on-
function regression model where a functional response and multiple functional
predictors consist of random curves with quadratic and interaction effects.
The direct estimation of a function-on-function regression model is usually an
ill-posed problem. To overcome this difficulty, in practice, the functional
data that belong to the infinite-dimensional space are generally projected
into a finite-dimensional space of basis functions. The function-on-function
regression model is converted to a multivariate regression model of the basis
expansion coefficients. In the estimation phase of the proposed method, the
functional variables are approximated by a finite-dimensional basis function
expansion method. We show that the partial least squares regression
constructed via a functional response, multiple functional predictors, and
quadratic/interaction terms of the functional predictors is equivalent to the
partial least squares regression constructed using basis expansions of
functional variables. From the partial least squares regression of the basis
expansions of functional variables, we provide an explicit formula for the
partial least squares estimate of the coefficient function of the function-on-
function regression model. Because the true forms of the models are generally
unspecified, we propose a forward procedure for model selection. The finite
sample performance of the proposed method is examined using several Monte
Carlo experiments and two empirical data analyses, and the results were found
to compare favorably with an existing method.

Keywords: Basis function; Function-on-function regression; Functional partial
least squares; Interaction effects; Quadratic term.

##  1 Introduction

Recent advances in computer storage and data collection have enabled
researchers in diverse branches of science such as, for instance,
chemometrics, meteorology, medicine and finance, recording data of
characteristics varying over a continuum (time, space, depth, wavelength,
etc). This is the case of data collected by spectrometer, rain gauges,
electroencephalographs, or even just by means of high-performance computer. In
all these cases, a number of subjects are observed densely over time, space or
both. Through the application of interpolation or smoothing techniques, these
data become functions. Thus, these data can be represented in a form of curve,
image or shape. Functional data analysis has arisen as a field of statistics
which provides statistical tools for analyzing this type of information
(Ramsay and Silverman, 2002, 2006; Ferraty and Vieu, 2006; Horvath and
Kokoszka, 2012; Cuevas, 2014). Hsing and Eubank (2015) and Kokoszka and
Reimherr (2017) provide a thorough overview of the research on theoretical
results and case studies of functional data analysis methods. The focus of
this study will be on constructing function-on-function regression models,
where a functional response and multiple functional predictors both consist of
random curves.

Function-on-function regression models have been well studied in the
literature (see e.g., Ramsay and Dalzell, 1991; Yao et al., 2005; Ramsay and
Silverman, 2006; Müller and Yao, 2008; Matsui et al., 2009; He et al., 2010;
Wang, 2014; Ivanescu et al., 2015; Chiou et al., 2016, and references
therein). The main objective of such models is to explain the association
between the functional response and one or more functional predictors. Most of
the function-on-function regression models that have been proposed concentrate
on main effects only. Let $\mathcal{Y}(t)$ and $\mathcal{X}_{m}(s)$
($m=1,\cdots,M$) respectively denote a functional response and $M$ sets of
functional predictors, where $s\in S$ and $t\in T$ are closed and bounded
intervals on the real line. For ease of notation, $s$ and $t$ can be
standardized to the interval $[0,1]$. The function-on-function regression
model with main effects is defined as follows:

(1.1)
$\mathcal{Y}(t)=\beta_{0}(t)+\sum_{m=1}^{M}\int_{0}^{1}\mathcal{X}_{m}(s)\beta_{m}(s,t)ds+\epsilon(t),$

where $\beta_{0}(t)$ is the intercept function, $\beta_{m}(s,t)$ is the
bivariate coefficient function, and $\epsilon(t)$ is the random error function
with mean zero.

A regression model with only main effects is less flexible compared with the
one that incorporates the quadratic term and interaction effects. A scalar-on-
function regression model with the quadratic term of the functional predictor
was proposed by Yao and Müller (2010). This model was later extended to the
case where there is more than one functional predictor by Fuchs et al. (2015)
and Usset et al. (2016), who consider interactions between multiple functional
predictors. Matsui (2020) proposed a simple function-on-function regression
model with a quadratic term of the functional predictor as follows:

(1.2)
$\mathcal{Y}(t)=\beta_{0}(t)+\int_{0}^{1}\mathcal{X}(s)\beta(s,t)ds+\int_{0}^{1}\int_{0}^{1}\mathcal{X}(s)\mathcal{X}(r)\gamma(s,r,t)dsdr+\epsilon(t),$

where $\gamma(s,r,t)$ is the three-dimensional coefficient function of the
quadratic term. Matsui (2020) estimated model parameters by the penalized
maximum likelihood (PML) method. Sun and Wang (2020) proposed a functional
principal component (PC) regression to estimate the parameters in model (1.2)
and provided a procedure to test the significance of the quadratic term. Luo
and Qi (2019) proposed a multiple function-on-function regression model with
quadratic and interaction effects as follows:

(1.3)
$\mathcal{Y}(t)=\beta_{0}(t)+\sum_{m\in\mathcal{M}}\int_{0}^{1}\mathcal{X}_{m}(s)\beta_{m}(s,t)ds+\sum_{m,n\in\mathcal{I}}\int_{0}^{1}\int_{0}^{1}\mathcal{X}_{m}(s)\mathcal{X}_{n}(r)\gamma_{mn}(s,r,t)dsdr+\epsilon(t),$

where the error function $\epsilon(t)$ has mean zero and is independent of
$\mathcal{X}_{m}(s)$’s, and $\mathcal{M}$ and $\mathcal{I}$, respectively,
denote the index sets of main and quadratic/interaction effects. The index
sets $\mathcal{M}$ and $\mathcal{I}$ are assumed to be any subset of
$\left\\{1,\cdots,M\right\\}$ and $\left\\{\left(m,n\right):1\leq m\leq n\leq
M\right\\}$, respectively. In addition, Luo and Qi (2019) proposed a stepwise
procedure for model selection. The numerical results produced from the
previous studies have shown that functional regression models including
quadratic term and interaction effects perform better than standard functional
regression models in the presence of interaction.

We propose a functional partial least squares (PLS) regression to estimate the
multiple function-on-function regression model with quadratic and interaction
effects given in (1.3). Several functional PLS methods have been proposed to
estimate parameters of the functional regression models. For example, Preda
and Saporta (2005), Reiss and Ogden (2007), Krämer et al. (2008), Aguilera et
al. (2010), and Aguilera et al. (2016) proposed functional PLS methods for
scalar-on-function regression models. Delaigle and Hall (2012b) examined the
theoretical properties of functional PLS methods in scalar-on-function
regression models. In addition, Febrero-Bande et al. (2017) compared the
finite sample performances of these methods and discussed their advantages and
disadvantages. Hyndman and Shang (2009) proposed a functional PLS regression
based on a lagged functional predictor and a functional response for
forecasting functional time series. Escabias et al. (2007) proposed a
functional PLS method for logit regression. Functional PLS has also been used
for other regression problems. For example, Ferraty and Vieu (2006) defined a
functional PLS-based semi-metric for nonparametric functional predictor
variables and Preda et al. (2007) and Delaigle and Hall (2012a) proposed
functional PLS-based functional data classification methods. Our method
differs from these existing functional PLS methods in terms of the number of
functional predictors in the model: while the existing functional PLS methods
allow only one functional predictor in the model, the proposed method allows
for more than one functional predictor and their quadratic/interaction
effects. The main contributions of the proposed method include:  i)
vectorization of function values,  ii) an application of PLS (for vector-
valued predictors and responses), and  iii) a forward stepwise variable
selection.  The interaction terms, including the quadratic terms, are treated
as linear terms once they are appropriately transformed.

We use the idea of functional PLS proposed by Preda and Saporta (2005) and
Aguilera et al. (2010) to estimate model (1.3) in Section 2. In our proposed
method, all the functional predictors and their quadratic/interaction terms
are first stacked into a two-dimensional function
($\boldsymbol{\mathcal{X}}(s,r)$; see (2.4)). Then, a functional PLS
regression of $\mathcal{Y}(t)$ on $\boldsymbol{\mathcal{X}}(s,r)$ is
considered. Functional random variables, by their nature, belong to infinite-
dimensional space. Thus, the direct estimation of function-on-function
regression models is generally an ill-posed problem. One remedy to estimate
function-on-function regression models is the approximation of functional data
via finite-dimensional basis function expansion methods. In this way, a
function-on-function regression model is converted into a multivariate
regression model of basis expansion coefficients. Then, the functional PLS
regression of $\mathcal{Y}(t)$ on $\boldsymbol{\mathcal{X}}(s,r)$ is estimated
using the basis expansion approximations of $\mathcal{Y}(t)$ and
$\boldsymbol{\mathcal{X}}(s,r)$. We show the equality of the functional PLS
and the PLS regression constructed via basis expansions of $\mathcal{Y}(t)$
and $\boldsymbol{\mathcal{X}}(s,r)$ using the metrics associated with basis
functions. In addition, we provide an explicit formula for the coefficient
function of the functional PLS regression model from the PLS estimate of the
regression model constructed via basis coefficient expansions. When
considering model (1.3), the form of the true model is usually unspecified and
the index sets $\mathcal{M}$ and $\mathcal{I}$ are unknown in practice (Luo
and Qi, 2019). To this end, we further propose a forward procedure to
determine the significant functional predictors and their
quadratic/interaction terms. Our numerical results, which will be discussed in
detail in Sections 3 and 4, reveal that the proposed method based on the
forward procedure produces better prediction performance compared with the
full model, which includes all the predictors and their quadratic and
interaction terms. Section 5 concludes the paper.

##  2 Methodology

Let us consider a set of random functional response
$\\{Y_{1}(t),\dots,Y_{N}(t)\\}$ containing $N$ independent realization
$\mathcal{Y}_{i}$ and $M$ sets of random functional predictors
$\\{\mathcal{X}_{mi}(s);m=1,\cdots,M,i=1,\cdots,N\\}$. With these
realizations, we aim to construct the multiple function-on-function regression
model with quadratic and interaction effects given in (1.3). We assume that
the functional random variables are the elements of square-integrable and
real-valued functions $\mathcal{L}_{2}[0,1]$, and are of second-order
stochastic processes. Without loss of generality, we assume that
$\mathbb{E}\left[\mathcal{Y}(t)\right]=0$ and
$\mathbb{E}\left[\mathcal{X}_{m}(s)\right]=0$ for $m=1,\cdots,M$ and $\forall
t,s\in[0,1]$, and the functional data pairs
$\left\\{\mathcal{X}_{mi}(s),Y_{i}(t)\right\\}$ are fully observed. For the
sake of clarity, we consider the full model where all the functional
predictors and their quadratic/interaction terms are included in the model;
that is $m=1,\cdots,M$ and
$(m,n)=(1,1),(1,2),\cdots,(1,M),\allowbreak(2,3),\cdots,(M-1,M),(M,M)$. Using
the continuity of functional objects, we express model (1.3) as follows:

$\displaystyle\mathcal{Y}(t)$
$\displaystyle=\sum_{m\in\mathcal{M}}\int_{0}^{1}\mathcal{X}_{m}(s)\beta_{m}(s,t)ds+\sum_{m\in\mathcal{M}}\sum_{n\in\mathcal{M}}\int_{0}^{1}\int_{0}^{1}\mathcal{X}_{m}(s)\mathcal{X}_{n}(r)\gamma_{mn}(s,r,t)dsdr+\epsilon(t),$
$\displaystyle\mathcal{Y}(t)$
$\displaystyle=\int_{0}^{1}\sum_{m\in\mathcal{M}}\mathcal{X}_{m}(s)\beta_{m}(s,t)ds+\int_{0}^{1}\sum_{m\in\mathcal{M}}\mathcal{X}_{m}(s)\int_{0}^{1}\sum_{n\in\mathcal{M}}\mathcal{X}_{n}(r)\gamma_{mn}(s,r,t)dsdr+\epsilon(t),$
$\displaystyle\mathcal{Y}(t)$
$\displaystyle=\int_{0}^{1}\sum_{m\in\mathcal{M}}\mathcal{X}_{m}(s)\left[\beta_{m}(s,t)+\int_{0}^{1}\sum_{n\in\mathcal{M}}\mathcal{X}_{n}(r)\gamma_{mn}(s,r,t)dr\right]ds+\epsilon(t),$
$\displaystyle\mathcal{Y}(t)$
$\displaystyle=\sum_{m\in\mathcal{M}}\int_{0}^{1}\mathcal{X}_{m}(s)\left[\beta_{m}(s,t)+\sum_{n\in\mathcal{M}}\int_{0}^{1}\mathcal{X}_{n}(r)\gamma_{mn}(s,r,t)dr\right]ds+\epsilon(t).$

Let
$\boldsymbol{\mathcal{X}}(s,r):\mathbb{R}^{2}\rightarrow\mathbb{R}^{M(M+3)/2}$
be a two-dimensional vector-valued function as follows:

$\displaystyle\boldsymbol{\mathcal{X}}(s,r)$
$\displaystyle=(\mathcal{X}_{1}(s),\cdots,\mathcal{X}_{M}(s)$
$\displaystyle\mathcal{X}_{1}(s)\mathcal{X}_{1}(r),\mathcal{X}_{1}(s)\mathcal{X}_{2}(r),\cdots,\mathcal{X}_{1}(s)\mathcal{X}_{M}(r)$
$\displaystyle
X_{2}(s)\mathcal{X}_{2}(r),\cdots,X_{M-1}(s)X_{M}(r),\mathcal{X}_{M}(s)\mathcal{X}_{M}(r))^{\top}$

Herein, the $\boldsymbol{\mathcal{X}}(s,r)$ is considered as a single two-
dimensional vector-valued function but including $M(M+3)/2$ ($M$ sets of
functional predictors and $M\times(M+1)/2$ sets of functional
quadratic/interaction terms) sets of functional predictors and their
quadratic/interaction terms. Hereafter, we will consider the following
regression model:

(2.4)
$\mathcal{Y}(t)=\int_{0}^{1}\int_{0}^{1}\boldsymbol{\mathcal{X}}(s,r)\Xi(s,r,t)dsdr+\epsilon(t),$

where $\Xi(s,r,t)$ is a single three-dimensional coefficient function
comprising $M(M+3)/2$ sets of coefficient functions of both main and
quadratic/interaction effects, as follows:

$\displaystyle\Xi(s,r,t)$
$\displaystyle=(\beta_{1}(s,t),\cdots,\beta_{M}(s,t)$
$\displaystyle\gamma_{11}(s,r,t),\gamma_{12}(s,r,t),\cdots,\gamma_{1M}(s,r,t)$
$\displaystyle\gamma_{22}(s,r,t),\cdots,\gamma_{M-1,M}(s,r,t),\gamma_{MM}(s,r,t))^{\top}$

The coefficient identifiability of (2.4) is discussed in Appendix.

In practice, the random functions of a functional dataset
$\left\\{x_{i}(t):i=1,\cdots,N\right\\}$ are observed in the set of finite
discrete time points $\left\\{t_{\ell}:\ell=1,\cdots,L\right\\}$; that is
$\left\\{x_{i}(t_{\ell})\right\\}$. Thus, the discretely observed data points
are first converted into their functional forms using a suitable method such
as basis function expansion (Ramsay and Silverman, 2006) and nonparametric
smoothing of functions (Ferraty and Vieu, 2006) before fitting the functional
regression model. In this study, we consider the basis function expansion
approach. Briefly, a random function $x(t)$ can be approximated by the linear
combinations of basis functions, $\phi_{k}(t)$, and associated coefficients,
$c_{k}$, for a sufficiently large number of basis functions, $K$:

$x(t)\approx\sum_{k=1}^{K}c_{k}\phi_{k}(t).$

Several basis function expansion methods such as $B$-spline, Fourier, radial,
and wavelet bases have been proposed (see, e.g., Ramsay and Silverman, 2006).
In this study, we consider the $B$-spline basis function expansion method.

Let $K_{\mathcal{Y}}$ and $K_{\boldsymbol{\mathcal{X}}}$ denote the number of
basis functions used for approximating $\mathcal{Y}(t)$ and
$\boldsymbol{\mathcal{X}}(s,r)$, respectively. They are expressed by basis
function expansion as follows:

(2.5) $\displaystyle\mathcal{Y}(t)$
$\displaystyle=\sum_{k=1}^{K_{\mathcal{Y}}}c_{k}\phi_{k}(t)=\boldsymbol{c}^{\top}\boldsymbol{\Phi}(t),$
(2.6) $\displaystyle\boldsymbol{\mathcal{X}}(s,r)$
$\displaystyle=\sum_{j=1}^{K_{\boldsymbol{\mathcal{X}}}}\sum_{l=0}^{K_{\boldsymbol{\mathcal{X}}}}d_{jl}\psi_{jl}(s,r)=\boldsymbol{d}^{\top}\boldsymbol{\Psi}(s,r),$

where $\boldsymbol{\Phi}(t)$ and $\boldsymbol{\Psi}(s,r)$ respectively are the
$K_{\mathcal{Y}}$ and $K_{\boldsymbol{\mathcal{X}}}\times M(M+3)/2$
dimensional basis functions vector/matrix, as follows:

$\displaystyle\boldsymbol{\Phi}(t)$
$\displaystyle=(\phi_{1}(t),\cdots,\phi_{K_{\mathcal{Y}}}(t))^{\top}$
$\displaystyle\boldsymbol{\Psi}(s,r)$
$\displaystyle=(\psi_{1}(s),\cdots,\psi_{M}(s)$
$\displaystyle\psi_{1}(s)\psi_{1}(r),\psi_{1}(s)\psi_{2}(r),\cdots,\psi_{1}(s),\psi_{M}(r)$
$\displaystyle\psi_{2}(s)\psi_{2}(r),\cdots,\psi_{M-1}(s)\psi_{M}(r),\psi_{M}(s)\psi_{M}(r))^{\top}$

with
$\psi_{i}(s)=\left(\psi_{i1},\cdots,\psi_{iK_{\boldsymbol{\mathcal{X}}}}\right)^{\top}$
and $\boldsymbol{c}$ and $\boldsymbol{d}$ are the $K_{\mathcal{Y}}$ and
$K_{\boldsymbol{\mathcal{X}}}$ dimensional coefficient vectors, respectively.
In addition, let us assume that the three-dimensional coefficient function
$\Xi(s,r,t)$ in (2.4) is expressed by basis expansion as follows:

(2.7)
$\Xi(s,r,t)=\sum_{j=1}^{K_{\boldsymbol{\mathcal{X}}}}\sum_{l=1}^{K_{\boldsymbol{\mathcal{X}}}}\sum_{k=1}^{K_{\mathcal{Y}}}\xi_{jlk}\psi_{jl}(s,r)\phi_{k}(t)=\boldsymbol{\Psi}^{\top}(s,r)\Gamma_{(3)}^{\top}\boldsymbol{\Phi}(t),$

where $\boldsymbol{\Gamma}_{(3)}$ is a $K_{\mathcal{Y}}\times
K_{\boldsymbol{\mathcal{X}}}^{2}$ dimensional matrix obtained by matricizing a
three-dimensional $K_{\boldsymbol{\mathcal{X}}}\times
K_{\boldsymbol{\mathcal{X}}}\times K_{\mathcal{Y}}$ tensor
$\underline{\boldsymbol{\Gamma}}=\left(\xi_{jlk}\right)_{jlk}$ with respect to
the 3rd array (de Lathuauwer et al., 2000; Matsui, 2020). Without loss of
generality, we further assume that the error function in (2.4) is represented
as a basis function expansion as follows:

$\epsilon(t)=\boldsymbol{e}\boldsymbol{\Phi}(t),$

where
$\boldsymbol{e}=\left(\boldsymbol{e}_{1},\cdots,\boldsymbol{e}_{N}\right)^{\top}$
is an $N$-dimensional vector comprising random variables
$\boldsymbol{e}_{i}=\left(e_{i1},\cdots,e_{iK_{\mathcal{Y}}}\right)^{\top}$.
Now, replacing (2.5) to (2.7) with (2.4) yields:

$\displaystyle\boldsymbol{c}^{\top}\boldsymbol{\Phi}(t)$
$\displaystyle=\boldsymbol{d}^{\top}\boldsymbol{\Psi}(s,r)\boldsymbol{\Psi}^{\top}(s,r)\Gamma_{(3)}^{\top}\boldsymbol{\Phi}(t)+\boldsymbol{e}\boldsymbol{\Phi}(t),$
(2.8) $\displaystyle\boldsymbol{c}$
$\displaystyle=\boldsymbol{d}^{\top}\boldsymbol{\Psi}\Gamma_{(3)}^{\top}+\boldsymbol{e},$

where
$\boldsymbol{\Psi}=\int_{0}^{1}\int_{0}^{1}\boldsymbol{\Psi}(s,r)\boldsymbol{\Psi}^{\top}(s,r)dsdr$.
Thus, the problem of estimating the function-on-function regression model in
(2.4) is reduced to estimating the multivariate regression model (2.8) with
coefficient matrix $\Gamma_{(3)}^{\top}$.

When considering (2.8), the dimensions of the matrices increase exponentially
when a large number of predictors are used in the model or when a large number
of basis functions are used to approximate the functional objects. In such
cases, an ill-posed problem (multicollinearity) arises and estimating
$\Gamma_{(3)}^{\top}$ using traditional methods, such as least squares and
maximum likelihood results in infinitely many solutions. Regularisation
methods such as the PML method proposed by Matsui (2020) may be used to
overcome this problem and to produce an accurate estimate for
$\Gamma_{(3)}^{\top}$. However, it is computationally time-consuming, and it
may not be possible to obtain the PML estimate of $\Gamma_{(3)}^{\top}$ using
a computer with standard memory. One remedy to obtain a stable estimate for
$\Gamma_{(3)}^{\top}$ is to use dimension-reduction techniques, such as PC or
PLS regression models. PC regression has been proposed by Aguilera et al.
(1999) and Yao et al. (2005) to estimate standard function-on-function
regression models and by Sun and Wang (2020) to estimate univariate function-
on-function quadratic regression models. However, this method does not take
into account the relationship between response and predictor variables when
deciding the PCs. Further Aguilera et al. (2010) showed that in the functional
case, the parameter function estimated with PLS is much more accurate than
with PC. Therefore, we consider functional PLS regression to estimate model
(2.4).

###  2.1 PLS for the functional interaction regression model

In the multivariate case, the main goal of PLS regression is to construct a
linear model between an $n\times m$ dimensional fixed full-rank design matrix
$\boldsymbol{X}=\left(X_{1},\cdots,X_{m}\right)^{\top}$ and an $n\times r$
dimensional matrix for response variables
$\boldsymbol{Y}=\left(Y_{1},\cdots,Y_{r}\right)^{\top}$ as follows:

$\boldsymbol{Y}=\boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{E},$

where $\boldsymbol{\beta}$ and $\boldsymbol{E}$ are the $m\times r$
dimensional coefficient matrix and the error matrix with dimensions $n\times
r$, respectively. The PLS regression produces orthogonal latent components
$\boldsymbol{\mathbb{T}}$ formed as linear combinations of predictor matrix
$\boldsymbol{X}$ and an appropriate weight matrix $\boldsymbol{W}$ such that
$\boldsymbol{\mathbb{T}}=\boldsymbol{X}\boldsymbol{W}$. Herein, the weight
matrix is computed so that it maximizes the covariance between
$\boldsymbol{Y}$ and $\boldsymbol{\mathbb{T}}$. Then, an ordinary least
squares regression is used to capture the relationship between
$\boldsymbol{Y}$ and $\boldsymbol{\mathbb{T}}$ such that
$\boldsymbol{Y}=\boldsymbol{\mathbb{T}}\boldsymbol{Q}+\boldsymbol{E}$. In
turn, it produces the loading matrix $\boldsymbol{Q}$ for $\boldsymbol{Y}$.
Then, the PLS regression is constructed as follows:

$\boldsymbol{Y}=\boldsymbol{X}\boldsymbol{\Omega}+\boldsymbol{E},$

where $\boldsymbol{\Omega}=\boldsymbol{W}\boldsymbol{Q}$.

Several functional versions of PLS regression have been proposed to explore
the relationship between a functional predictor and a scalar and/or functional
response (see e.g., Preda and Saporta, 2005; Hyndman and Shang, 2009; Aguilera
et al., 2010, 2016; Beyaztas and Shang, 2020). Similar to the discrete case,
the functional PLS method produces orthogonal latent components as linear
combinations of the functional predictors by maximizing the squared covariance
between the response and orthogonal latent components (Tucker, 1938). The PLS
components are determined by an iterative process. In this study, we extend
the functional PLS methods proposed by Preda and Saporta (2005) and Aguilera
et al. (2010) to estimate the model given in (2.4). Our extension differs from
the methods of Preda and Saporta (2005) and Aguilera et al. (2010) in the
context of the number of functional predictors in the model. The original
functional PLS methods allow for only one functional predictor in the model.
However, our proposed method is possible to construct a regression model with
more than one functional predictor and their quadratic/interaction effects.
The proposed method first stacks all the main and quadratic/interaction
effects into a two-dimensional function, $\boldsymbol{\mathcal{X}}(s,r)$, it
then applies a PLS method on the vectorized functional variables to estimate
the model parameters.

The PLS components of regression model (2.4) can be obtained as solutions of
Tucker’s criterion extended to functional-type data as follows:

$\underset{\begin{subarray}{c}\kappa\in\mathcal{L}_{2}([0,1]\times[0,1]),~{}\|\kappa\|_{\mathcal{L}_{2}[0,1]\times[0,1]}=1\\\
\zeta\in\mathcal{L}_{2}[0,1],~{}\|\zeta\|_{\mathcal{L}_{2}[0,1]}=1\end{subarray}}{\max}\text{Cov}^{2}\left(\int_{0}^{1}\int_{0}^{1}\boldsymbol{\mathcal{X}}(s,r)\kappa(s,r)dsdr,~{}\int_{0}^{1}\mathcal{Y}(t)\zeta(t)dt\right),$

where $\left\|\cdot\right\|_{\mathcal{L}_{2}}$ denotes $\mathcal{L}_{2}$ norm,
which is approximated by the Riemann sum (Luo and Qi, 2019). The functional
PLS components are also the eigenvectors of the Escoufier operator (Escoufier,
1970). Let $W^{\boldsymbol{\mathcal{X}}}$ and $W^{\mathcal{Y}}$, respectively,
denote the Escoufier’s operators associated with
$\boldsymbol{\mathcal{X}}(s,r)$, with respect to $\mathcal{Y}(t)$ as follows:

$W^{\boldsymbol{\mathcal{X}}}=\int_{0}^{1}\int_{0}^{1}\mathbb{E}\left[\boldsymbol{\mathcal{X}}(s,r)Z\right]\boldsymbol{\mathcal{X}}(s,r)dsdr,\qquad
W^{\mathcal{Y}}=\int_{0}^{1}\mathbb{E}\left[\mathcal{Y}(t)Z\right]\mathcal{Y}(t)dt,\qquad\forall
Z\in\mathcal{L}_{2}[0,1].$

Then, the first PLS component, say $\eta_{1}$, of model (2.4) is equal to the
eigenvector corresponding to the largest eigenvalue ($\lambda_{\max}$) of
$W^{\boldsymbol{\mathcal{X}}}W^{\mathcal{Y}}$:

$W^{\boldsymbol{\mathcal{X}}}W^{\mathcal{Y}}\eta_{1}=\lambda_{\max}\eta_{1}.$

Accordingly, the first PLS component is defined as follows:

$\eta_{1}=\int_{0}^{1}\int_{0}^{1}\kappa_{1}(s,r)\boldsymbol{\mathcal{X}}(s,r)dsdr,$

where the weight function $\kappa_{1}(s,r)$ associated with $\eta_{1}$ is
given by:

$\kappa_{1}(s,r)=\frac{\int_{0}^{1}\mathbb{E}\left[\mathcal{Y}(t)\boldsymbol{\mathcal{X}}(s,r)\right]dt}{\sqrt{\int_{0}^{1}\int_{0}^{1}\left(\int_{0}^{1}\mathbb{E}\left[\mathcal{Y}(t)\boldsymbol{\mathcal{X}}(s,r)\right]dt\right)^{2}dsdr}}.$

PLS is an iterative stepwise procedure. Let $h=1,2,\cdots,$ denote the
iteration number. Define $\boldsymbol{\mathcal{X}}_{h}(s,r)$ and
$\mathcal{Y}_{h}(t)$ by the residuals obtained from the following regressions:

$\displaystyle\boldsymbol{\mathcal{X}}_{h}(s,r)$
$\displaystyle=\boldsymbol{\mathcal{X}}_{h-1}(s,r)-p_{h}(s,r)\eta_{h},$
$\displaystyle\mathcal{Y}_{h}(t)$
$\displaystyle=\mathcal{Y}_{h-1}(t)-\zeta_{h}(t)\eta_{h},$

where $\boldsymbol{\mathcal{X}}_{0}(s,r)=\boldsymbol{\mathcal{X}}(s,r)$,
$\mathcal{Y}_{0}(t)=\mathcal{Y}(t)$,
$p_{h}(s,r)=\frac{\mathbb{E}\left[\boldsymbol{\mathcal{X}}_{h-1}(s,r)\eta_{h}\right]}{\mathbb{E}\left[\eta_{h}^{2}\right]}$,
and
$\zeta_{h}(t)=\frac{\mathbb{E}\left[\mathcal{Y}_{h-1}(t)\eta_{h}\right]}{\mathbb{E}\left[\eta_{h}^{2}\right]}$.
Then, the $h^{\text{th}}$ PLS component, $\eta_{h}$, is equal to the
eigenvector corresponding to largest eigenvalue of
$W_{h-1}^{\boldsymbol{\mathcal{X}}}W_{h-1}^{\mathcal{Y}}$:

$W_{h-1}^{\boldsymbol{\mathcal{X}}}W_{h-1}^{\mathcal{Y}}\eta_{h}=\lambda_{\max}\eta_{h},$

where $W_{h-1}^{\boldsymbol{\mathcal{X}}}$ and $W_{h-1}^{\mathcal{Y}}$ are the
Escoufier’s operators associated with $\boldsymbol{\mathcal{X}}_{h-1}(s,r)$
and $\mathcal{Y}_{h-1}(t)$, respectively. That is, $\eta_{h}$ is given as
follows:

$\eta_{h}=\int_{0}^{1}\int_{0}^{1}\kappa_{h}(s,r)\boldsymbol{\mathcal{X}}_{h-1}(s,r)dsdr,$

where the weight function associated with $\eta_{h}$ is given by:

$\kappa_{h}(s,r)=\frac{\int_{0}^{1}\mathbb{E}\left[\mathcal{Y}_{h-1}(t)\boldsymbol{\mathcal{X}}_{h-1}(s,r)\right]dt}{\sqrt{\int_{0}^{1}\int_{0}^{1}\left(\int_{0}^{1}\mathbb{E}\left[\mathcal{Y}_{h-1}(t)\boldsymbol{\mathcal{X}}_{h-1}(s,r)\right]dt\right)^{2}dsdr}}.$

Following Aguilera et al. (2010), the properties of the PLS components of
model (2.4) can be summarized by the following proposition:

######  Proposition 2.1.

For any $h\geq 1$

  * i)

$\left\\{\eta_{h}\right\\}_{h\geq 1}$ forms an orthogonal system in the linear
space spanned by $\boldsymbol{\mathcal{X}}(s,r)$,

  * ii)

$\mathcal{Y}(t)=\zeta_{1}(t)\eta_{1}+\cdots+\zeta_{h}(t)\eta_{h}+\mathcal{Y}_{h}(t)$,

  * iii)

$\boldsymbol{\mathcal{X}}(s,r)=p_{1}(s,r)\eta_{1}+\cdots+p_{h}(s,r)\eta_{h}+\boldsymbol{\mathcal{X}}_{h}(s,r)$,

  * iv)

$\mathbb{E}\left[\mathcal{Y}_{h}(t)\eta_{j}\right]=0$, $\forall j=1,\cdots,h$,

  * v)

$\mathbb{E}\left[\boldsymbol{\mathcal{X}}_{h}(s,r)\eta_{j}\right]=0$, $\forall
j=1,\cdots,h$.

###  2.2 PLS for basis expansion of the functional interaction regression
model

Functional PLS components are estimated using an observed sample of the
functional dataset. However, the random functions of functional variables are
observed in the set of finite discrete time points. In this case, the main
problem related to the sample estimation of PLS components is that of
estimating the Escoufier’s operators from discrete time observations (Aguilera
et al., 2010). To overcome this problem, we consider the basis expansions of
the observed functional dataset, and construct the PLS regression using the
basis expansion coefficients of $\mathcal{Y}(t)$ and
$\boldsymbol{\mathcal{X}}(s,r)$. Let us revisit model (2.4) and its basis
function expansion in (2.8). Then, we provide the following proposition to
show the relationship between the functional PLS of $\mathcal{Y}(t)$ on
$\boldsymbol{\mathcal{X}}(s,r)$ and the PLS of $\boldsymbol{c}$ on
$\boldsymbol{d}$. The proof of the proposition has been relegated to the
Appendix.

######  Proposition 2.2.

Under the assumption that the functional variables $\mathcal{Y}(t)$ and
$\boldsymbol{\mathcal{X}}(s,r)$ are spanned with finite number of basis
functions, let $\boldsymbol{\Phi}$ and $\boldsymbol{\Psi}$, respectively,
denote $K_{\mathcal{Y}}\times K_{\mathcal{Y}}$ and
$K_{\boldsymbol{\mathcal{X}}}^{2}\times K_{\boldsymbol{\mathcal{X}}}^{2}$
dimensional symmetric matrices with entries the inner products of the basis
functions $\int_{0}^{1}\boldsymbol{\Phi}(t)\boldsymbol{\Phi}^{\top}(t)dt$ and
$\int_{0}^{1}\int_{0}^{1}\boldsymbol{\Psi}(s,r)\boldsymbol{\Psi}^{\top}(s,r)dsdr$,
respectively. Also, let $\boldsymbol{\Phi}^{1/2}$ and
$\boldsymbol{\Psi}^{1/2}$ denote the square roots of $\boldsymbol{\Phi}$ and
$\boldsymbol{\Psi}$, respectively. Then, the functional PLS regression of
$\mathcal{Y}(t)$ on $\boldsymbol{\mathcal{X}}(s,r)$ is equivalent to the PLS
regression of $\boldsymbol{\Phi}^{1/2}\boldsymbol{c}$ on
$\boldsymbol{\Psi}^{1/2}\boldsymbol{d}$ in the sense that at each step $h$ of
the PLS algorithm, the same PLS components are obtained for both PLS
regressions.

The PLS approximation of the coefficient function $\Xi(s,r,t)$ in (2.4) can be
expressed in the bases $\boldsymbol{\Phi}(t)$ and $\boldsymbol{\Psi}(s,r)$.
Denote by
$\boldsymbol{\Theta}^{h}=\left(\Theta_{1}^{h},\cdots,\Theta_{K_{\boldsymbol{\mathcal{X}}}^{2}}^{h}\right)$
the regression coefficients of $\boldsymbol{\Phi}^{1/2}\boldsymbol{c}$ on
$\boldsymbol{\Psi}^{1/2}\boldsymbol{d}$ obtained by the PLS regression at step
$h$. Then, we have

$\displaystyle\boldsymbol{\Phi}^{1/2}\boldsymbol{c}$
$\displaystyle=\boldsymbol{d}\boldsymbol{\Psi}^{1/2}\boldsymbol{\Theta}^{h},$
$\displaystyle\widehat{\mathcal{Y}}(t)$
$\displaystyle=\int_{0}^{1}\int_{0}^{1}\boldsymbol{\mathcal{X}}(s,r)\widehat{\Xi}(s,r,t)dsdr,$

where

$\widehat{\Xi}(s,r,t)=\sum_{j=1}^{K_{\boldsymbol{\mathcal{X}}}}\sum_{l=1}^{K_{\boldsymbol{\mathcal{X}}}}\sum_{k=1}^{K_{\mathcal{Y}}}\left(\left(\boldsymbol{\Psi}^{1/2}\right)^{-1}\boldsymbol{\Theta}^{h}\left(\boldsymbol{\Phi}^{1/2}\right)^{-1}\right)_{(jl)}\psi_{jl}(s,r)\phi_{k}(t),\qquad
s,r,t\in\mathcal{L}_{2}[0,1].$

Herein, $\widehat{\Xi}(s,r,t)$ denotes the PLS approximation of the three-
dimensional coefficient function $\Xi(s,r,t)$. These results demonstrate that
the infinite-dimensional estimation problem can be reduced to a simple finite-
dimensional PLS regression setting using the particular metrics
$\boldsymbol{\Phi}$ and $\boldsymbol{\Psi}$ in the spaces of expansion
coefficients $\boldsymbol{c}$ and $\boldsymbol{d}$, respectively. Through the
basis expansion, one can implement a finite multivariate PLS regression
setting and a PLS algorithm, such as nonlinear iterative partial least squares
(NIPALS) of Wold (1974), the SIMPLS of de de Jong (1993), and the improved
kernel PLS of Dayal and MacGregor (1997) to approximate the three-dimensional
coefficient function $\Xi(s,r,t)$. In this study, the NIPALS algorithm is used
to obtain the PLS approximation of $\Xi(s,r,t)$.

###  2.3 Model selection

####  2.3.1 Variable selection procedure

The functional regression model in (2.4) includes too many interaction terms
when a large number of functional predictors are used in the model. However,
not all the functional terms may be significant for the model. Thus, a
variable selection procedure is needed to select the best subset of functional
predictors and quadratic/interaction terms. In this context, several criteria,
such as the Akaike information criterion, Bayesian information criterion, and
the predicted residual error sum of squares, have been proposed. However, a
variable selection procedure based on these criteria will be computationally
intensive, since the number of all possible models is very large when there
are a large number of predictors in the model. Various variable selection
procedures, such as backward elimination and forward selection, have been
proposed to overcome this problem. In this study, we propose the following
forward selection procedure to select the significant functional predictors
and quadratic/interaction terms.

  * Step 1.

We consider selecting the main effect terms first, since they are generally
more important for the model compared with the quadratic term and interaction
effects. To start with, first, $M$ function-on-function regression models,
each of which includes the common response and a predictor, are constructed as
follows:

$\mathcal{Y}(t)=\int_{0}^{1}\int_{0}^{1}\boldsymbol{\mathcal{X}}_{m}(s,r)\Xi_{m}(s,r,t)dsdr+\epsilon_{m}(t),$

where $\boldsymbol{\mathcal{X}}_{m}(s,r)=\mathcal{X}_{m}(s)$ and $r=0$ for
$m=1,\cdots,M$. Among these models, the one having the smallest mean squared
error (MSE),

$\text{MSE}=\frac{1}{N}\sum_{j=1}^{N}\left\|\mathcal{Y}_{j}(t)-\widehat{\mathcal{Y}}_{j}(t)\right\|^{2}_{\mathcal{L}_{2}},$

where $\widehat{\mathcal{Y}}_{j}(t)$ is the prediction of $\mathcal{Y}_{j}(t)$
for $j=1,\cdots,N$, is chosen as an initial model. Let
$\boldsymbol{\mathcal{X}}^{(1)}(s,r)$ and $\text{MSE}^{(1)}$ denote the
predictor variable in the initial model and the MSE obtained from this model,
respectively. Then, $M-1$ function-on-function regression models, each of
which includes the response and a predictor, are constructed as follows:

$\mathcal{Y}(t)=\int_{0}^{1}\int_{0}^{1}\boldsymbol{\mathcal{X}}_{m}(s,r)\Xi_{m}(s,r,t)dsdr+\epsilon_{m}(t),$

where
$\boldsymbol{\mathcal{X}}_{m}(s,r)=\left\\{\boldsymbol{\mathcal{X}}^{(1)}(s,r),\mathcal{X}_{m}(s)\right\\}$,
$\mathcal{X}_{m}(s)\neq\boldsymbol{\mathcal{X}}^{(1)}(s,r)$, and $r=0$ for
$m=1,\cdots,M-1$, and the MSE is calculated for each of these models. The
predictor variable that produces the smallest MSE,
$\boldsymbol{\mathcal{X}}^{(2)}(s,r)$, is chosen as the predictor variable for
the current model if $\text{MSE}^{(2)}<\text{MSE}^{(1)}$ where
$\text{MSE}^{(2)}$ is the calculated MSE when
$\boldsymbol{\mathcal{X}}^{(2)}(s,r)$ is used to estimate the current model.
This process is repeated until all the significant variables are included in
the model.

  * Step 2.

We consider an iterative procedure similar to that described in Step 1 to
select the significant quadratic and interaction effects terms. Let
$\mathcal{M}^{*}$ ($\mathcal{M}^{*}\subseteq\mathcal{M}$) denote the index set
of all significant main effect terms found in the previous step. Denote by
$M^{*}$ the number of items in $\mathcal{M}^{*}$, so there are
$M^{*}\times(M^{*}+1)/2$ quadratic/interaction terms to be considered. Let
$\text{MSE}^{(\mathcal{M}^{*})}$ denote the calculated MSE when all the
significant main effect terms are included in the model. To select
quadratic/interaction effect terms, first, $M^{*}\times(M^{*}+1)/2$ function-
on-function regression models including the response and a predictor are
constructed as follows:

$\mathcal{Y}(t)=\int_{0}^{1}\int_{0}^{1}\boldsymbol{\mathcal{X}}_{m}(s,r)\Xi_{m}(s,r,t)dsdr+\epsilon_{m}(t),$

where
$\boldsymbol{\mathcal{X}}_{m}(s,r)=\left[\left\\{\mathcal{X}_{m}(s)\right\\}_{m\in\mathcal{M}^{*}}\left(1+\mathcal{X}_{n}(r)\right)\right]$
for $n=1,\cdots,M^{*}$, and the MSE is calculated for each of these models.
Then, the first quadratic/interaction effect term that leads to the greatest
improvement to the $\text{MSE}^{(\mathcal{M}^{*})}$ is included in the model.
Other quadratic/interaction effect terms are iteratively selected to enter the
model based on their improvements to the MSE of the current model.

####  2.3.2 Selection of the optimum number of PLS components

To select the optimum number of PLS components, we consider a cross-
validation-based mean squared prediction error (MSPE) as follows:

  * 1)

First, the data are randomly divided into training and testing samples with
sample sizes $N_{train}$ and $N_{test}$, respectively. Note that in our
numerical analyses, the training and testing samples are determined based on
the 50% – 50% rule.

  * 2)

Using the training set, $B_{h}$ PLS regression models based on
$h=1,\cdots,B_{h}$ number of PLS components are constructed.

  * 3)

The values of the response variable in the testing sample are predicted using
the values of the predictor variables in the testing sample and the
coefficient matrices obtained from the previous step, and the MSPE is computed
for each $B_{h}$:

$\text{MSPE}(h)=\frac{1}{N_{test}}\sum_{i=1}^{N_{test}}\left\|\mathcal{Y}_{i}(t)-\widehat{\mathcal{Y}}_{i}^{h}(t)\right\|^{2}_{\mathcal{L}_{2}},$

where $\widehat{\mathcal{Y}}_{i}^{h}(t)$ is the prediction of
$\mathcal{Y}_{i}(t)$ based on $h$ PLS components.

Based on the above algorithm, the optimum number of PLS components is
equivalent to one which produces the minimum MSPE. Note that the cross-
validation approach should only be applied to the training set to select the
number of components. Within the training set, we can further split the
samples into a training and a test sample. We note that our procedure for the
selection of the optimum number of PLS components consists of two steps. In
the first step, the significant main, quadratic, and interaction effects using
a fixed number of PLS components, say 8, are determined as discussed in
Section 2.3.1. Then, in the second step, the cross-validation-based MSPE
criterion given above is used to determine the number of optimum PLS
components. Our numerical analyses (not reported in the paper) have shown that
the choice of a fixed number of PLS components in the first step does not have
a significant effect on the determination of main, quadratic, and interaction
effects when $h\geq 4$.

##  3 Numerical results

Several Monte-Carlo simulations were conducted to investigate the finite
sample performances of the proposed method, and the results were compared with
the method proposed by Luo and Qi (2019) (”LQ”, hereafter). Two different
simulations settings, which are modified versions proposed by Luo and Qi
(2019), were considered. Throughout the simulations, $m=5$ functional
predictors $\mathcal{X}_{m}$ were generated from the following process:

$\mathcal{X}_{m}(s)=\sum_{i=0}^{\text{Lag}}V_{m+i}(s)/\sqrt{\text{Lag}+1},$

where $V_{m}(s)\text{s}~{}(m=1,\cdots,9)$ were generated from the Gaussian
process with mean zero and a positive definite covariance function
$\boldsymbol{\Sigma}_{V}(s,s^{\prime})=e^{-100(s-s^{\prime})^{2}}$. Herein,
the parameter Lag ($\text{Lag}>0$) controls the correlation between the
predictor functions, and its large values correspond to high correlations
between the predictor functions. Two different correlation levels, namely
$\text{Lag}=2$ and $\text{Lag}=4$, were considered. The error process
$\epsilon(t)$ was generated from the normal distribution with mean zero and
variance four; $N(0,4)$. Then, the realizations of the response function were
generated as follows:

$\mathcal{Y}(t)=\sum_{m\in\mathcal{M}}\int_{0}^{1}\mathcal{X}_{m}(s)\beta_{m}(s,t)ds+\sum_{m,n\in\mathcal{I}}\int_{0}^{1}\int_{0}^{1}\mathcal{X}_{m}(s)\mathcal{X}_{n}(r)\gamma_{mn}(s,r,t)dsdr+\epsilon(t),$

where $\mathcal{M}=\left\\{2,3,4,5\right\\}$ and
$\mathcal{I}=\left\\{(2,2),(3,4),(4,5)\right\\}$ for the first simulation
setting and $\mathcal{M}=\left\\{1,2,4,5\right\\}$ and
$\mathcal{I}=\left\\{(1,1),(1,2),(1,5),(2,4),(4,5),(5,5)\right\\}$ for the
second simulation setting. For the first simulation setting, the following
coefficient functions were considered to generate the response:

$\displaystyle\beta_{2}(s,t)$ $\displaystyle=e^{-3(s-1)^{2}-5(t-0.5)^{2}},$
$\displaystyle\beta_{3}(s,t)$
$\displaystyle=e^{-5(s-0.5)^{2}-5(t-0.5)^{2}}+8e^{-5(s-1.5)^{2}-5(t-0.5)^{2}},$
$\displaystyle\beta_{4}(s,t)$ $\displaystyle=\sin(1.5\pi s)\sin(\pi t),$
$\displaystyle\beta_{5}(s,t)$ $\displaystyle=\sqrt{st},$
$\displaystyle\gamma_{22}(s,r,t)$ $\displaystyle=5sr\sqrt{t},$
$\displaystyle\gamma_{34}(s,r,t)$ $\displaystyle=5\cos(\pi s)\sin(2\pi
r)\cos(2\pi t),$ $\displaystyle\gamma_{45}$ $\displaystyle=0.5e^{s+2r-t}.$

For the second simulation setting, the response functions were generated using
the following coefficient functions:

$\displaystyle\beta_{1}(s,t)$ $\displaystyle=(s-2t)^{2}/3,$
$\displaystyle\beta_{2}(s,t)$
$\displaystyle=2\left(\ln(1+s)\right)^{2}\sin\left(2\pi(t-0.5)\right),$
$\displaystyle\beta_{4}(s,t)$
$\displaystyle=\left(\cos\left(1-s\right)+\sqrt{t}\right)/3,$
$\displaystyle\beta_{5}(s,t)$
$\displaystyle=(1+s)^{2}/\left(3(1+t^{2})\right),$
$\displaystyle\gamma_{11}(s,r,t)$ $\displaystyle=2(s+r)t^{2},$
$\displaystyle\gamma_{1,2}(s,r,t)$ $\displaystyle=0.01(s^{2}-r^{3}+t),$
$\displaystyle\gamma_{15}(s,r,t)$ $\displaystyle=0.01e^{2s-r+3t},$
$\displaystyle\gamma_{24}(s,r,t)$ $\displaystyle=0.01(2s-r+3t),$
$\displaystyle\gamma_{45}(s,r,t)$ $\displaystyle=0.01\frac{\ln(1+2s)}{1+t},$
$\displaystyle\gamma_{55}(s,r,t)$
$\displaystyle=\cos\left(\pi(s+r)\right)+3\sqrt{t}.$

Before fitting the regression model, the functional predictors were distorted
by the Gaussian noise $u(s)\sim N(0,4)$. All the functions for both the
response and predictors variables were generated at 100 equally spaced points
in the interval $[0,1]$. Examples of the true and noisy functions for the
generated response and predictor variables are presented in Figure 1.

Figure 1: Plots of the generated 20 true (black lines) and noisy (gray lines)
functions when Lag = 4. The functions were generated under the second
simulation setting.

For each simulation setting, MC = 250 Monte Carlo simulations were performed,
and in each simulation, $N=300$ functions were generated for all random
variables. The finite sample performances of the methods were evaluated using
three metrics: MSPE, root mean squared prediction error (RMSPE), and mean
absolute prediction error (MAPE). The performance metrics were calculated
using the following procedure;  1) the first $N_{train}=100$ functions of the
data were used to construct the model, and  2) the last $N_{test}=200$
functions of the data were used to evaluate the prediction performances of the
methods based on the constructed model using the first $N_{train}=100$
functions of the data:

MSPE
$\displaystyle=\frac{1}{\text{MC}}\frac{1}{N_{test}}\sum_{i=1}^{\text{MC}}\sum_{j=1}^{N_{test}}\left\|\mathcal{Y}_{j,i}^{\text{pred}}(t)-\mathcal{Y}_{j,i}(t)\right\|^{2}_{\mathcal{L}_{2}},$
RMSPE
$\displaystyle=\sqrt{\frac{1}{\text{MC}}\frac{1}{N_{test}}\sum_{i=1}^{\text{MC}}\sum_{j=1}^{N_{test}}\left\|\frac{\mathcal{Y}_{j,i}^{\text{pred}}(t)-\mathcal{Y}_{j,i}(t)}{\mathcal{Y}_{j,i}(t)}\right\|^{2}_{\mathcal{L}_{2}}},$
MAPE
$\displaystyle=\frac{1}{\text{MC}}\frac{1}{N_{test}}\sum_{i=1}^{\text{MC}}\sum_{j=1}^{N_{test}}\left\|\frac{\left|\mathcal{Y}_{j,i}^{\text{pred}}(t)-\mathcal{Y}_{j,i}(t)\right|}{\mathcal{Y}_{j,i}(t)}\right\|_{\mathcal{L}_{2}},$

where $\mathcal{Y}_{j,i}(t)$ and $\mathcal{Y}_{j,i}^{\text{pred}}(t)$ are the
generated true response functions for $j^{\text{th}}$ individual in
$i^{\text{th}}$ simulation and its prediction, respectively. The proposed PLS
and LQ methods were compared based on four models;  1) the full model, which
includes all main, quadratic, and interaction terms;  2) the true model, which
includes only the terms used in the data-generating process;  3) the selected
model, which includes the terms selected by the variable selection procedure;
and  4) the main effect model, which includes only all five main effect terms.

Before presenting our results, we note that Luo and Qi (2019) used 40
$B$-spline and $20\times 20=400$ tensor product $B$-spline basis functions for
estimating the model. For the proposed method, the number of basis functions
used to approximate the response and predictor variables was chosen from a
small set number of basis functions according to its prediction performance.
For this purpose, first, $K_{\mathcal{Y}}=\left[4,6,8,10\right]$ and
$K_{\boldsymbol{\mathcal{X}}}=\left[4,6,8,10,15\right]$ number of basis
functions were used to approximate the response and predictor variables,
respectively. Then, the number of basis functions which provide the smalles
MSEs were chosen to be used in the analyses. For the proposed method, the
boxplots of the selected number of basis functions for both the response and
predictor variables are presented in Figure 2.

Figure 2: Boxplots of the selected numbers of basis functions for the
predictor (X) and response (Y) variables from 250 Monte Carlo simulations. The
first and second columns, respectively, correspond to the results obtained
under the first and second simulation settings while the first and second rows
correspond to the case when the correlation parameter equals Lag = 2 and Lag =
4, respectively.

Our findings for the MSPE, RMSPE, and MAPE metrics are presented in Table
LABEL:tab:pm_all. The results demonstrate that:  1) The correlation parameter
Lag has slight effects on the error measures.  2) Compared with the main
effect models, both LQ and the proposed PLS methods generally produce better
prediction performances (for all performance metrics) when
quadratic/interaction effects are used in the model.  3) For the MSPE metric,
both the proposed and LQ methods produce competitive results under the first
simulation setting, while LQ outperforms the proposed PLS method when the data
are generated under the second simulation setting. For the RMSPE and MAPE
metrics, the proposed method shows better performance than that of LQ for both
simulation settings.  4) Table LABEL:tab:pm_all presents 12 different results
(two simulation settings $\times$ two lag parameters $\times$ three error
measures) for each specific model (main effect, full, true, and selected). We
observe that the selected LQ model has better prediction performances than
those of the full model in 58% of all results. In contrast, the selected
proposed method performs better than the full model in 91% of all results.

Table 1: Computed average MSPE, RMSPE, and MAPE values of the LQ and PLS methods. The subscripts ${}_{\text{main}}$, ${}_{\text{true}}$, ${}_{\text{full}}$, and ${}_{\text{selected}}$ correspond to the case where the model is estimated using the main effect, true, full, and selected models. The values given in brackets are the estimated standard errors for the calculated performance metrics. The smallest errors are highlighted in bold. Setting |  Lag |  Method |  MSPE |  RMSPE |  MAPE
---|---|---|---|---|---
Setting 1 |  Lag = 2 |  LQ${}_{\text{main}}$ |  0.428 (0.086) |  180.353 (326.525) |  4.519 (2.437)
|  |  LQ${}_{\text{full}}$ |  0.270 (0.067) |  190.236 (523.532) |  3.458 (3.807)
|  |  LQ${}_{\text{true}}$ |  0.234 (0.045) |  151.666 (274.335) |  3.531 (1.424)
|  |  LQ${}_{\text{selected}}$ |  0.281 (0.051) |  197.873 (634.627) |  3.370 (4.571)
|  |  PLS${}_{\text{main}}$ |  0.433 (0.087) |  80.861 (150.235) |  3.366 (1.051)
|  |  PLS${}_{\text{full}}$ |  0.311 (0.049) |  96.755 (284.596) |  2.376 (1.545)
|  |  PLS${}_{\text{true}}$ |  0.238 (0.036) |  67.417 (189.168) |  1.799 (1.053)
|  |  PLS${}_{\text{selected}}$ |  0.244 (0.048) |  64.563 (130.164) |  1.926 (1.015)
|  Lag = 4 |  LQ${}_{\text{main}}$ |  0.522 (0.128) |  229.012 (354.072) |  5.980 (2.809)
|  |  LQ${}_{\text{full}}$ |  0.284 (0.055) |  142.713 (266.679) |  3.194 (1.449)
|  |  LQ${}_{\text{true}}$ |  0.262 (0.049) |  126.574 (239.987) |  2.833 (1.410)
|  |  LQ${}_{\text{selected}}$ |  0.298 (0.061) |  114.923 (220.162) |  2.912 (1.098)
|  |  PLS${}_{\text{main}}$ |  0.527 (0.124) |  136.552 (217.661) |  3.871 (1.520)
|  |  PLS${}_{\text{full}}$ |  0.288 (0.054) |  45.590 (41.800) |  2.086 (0.992)
|  |  PLS${}_{\text{true}}$ |  0.258 (0.042) |  40.495 (38.179) |  1.889 (0.889)
|  |  PLS${}_{\text{selected}}$ |  0.260 (0.045) |  40.406 (41.207) |  2.142 (1.371)
Setting 2 |  Lag = 2 |  LQ${}_{\text{main}}$ |  0.273 (0.059) |  256.474 (392.914) |  6.065 (3.726)
|  |  LQ${}_{\text{full}}$ |  0.115 (0.035) |  139.166 (224.285) |  2.582 (1.616)
|  |  LQ${}_{\text{true}}$ |  0.101 (0.027) |  146.398 (230.747) |  2.341 (1.104)
|  |  LQ${}_{\text{selected}}$ |  0.103 (0.028) |  158.856 (382.274) |  2.932 (3.491)
|  |  PLS${}_{\text{main}}$ |  0.255 (0.054) |  169.522 (250.788) |  4.877 (2.072)
|  |  PLS${}_{\text{full}}$ |  0.157 (0.029) |  89.768 (192.956) |  2.847 (1.058)
|  |  PLS${}_{\text{true}}$ |  0.140 (0.028) |  89.620 (195.483) |  2.724 (1.027)
|  |  PLS${}_{\text{selected}}$ |  0.118 (0.027) |  72.743 (137.590) |  2.452 (1.011)
|  Lag = 4 |  LQ${}_{\text{main}}$ |  0.284 (0.075) |  282.228 (490.338) |  6.713 (4.571)
|  |  LQ${}_{\text{full}}$ |  0.129 (0.028) |  160.912 (264.747) |  3.385 (1.919)
|  |  LQ${}_{\text{true}}$ |  0.111 (0.027) |  138.492 (222.025) |  3.001 (1.258)
|  |  LQ${}_{\text{selected}}$ |  0.109 (0.024) |  157.069 (268.081) |  3.218 (1.769)
|  |  PLS${}_{\text{main}}$ |  0.267 (0.069) |  157.322 (280.071) |  5.045 (2.326)
|  |  PLS${}_{\text{full}}$ |  0.154 (0.040) |  98.047 (168.614) |  3.174 (1.950)
|  |  PLS${}_{\text{true}}$ |  0.136 (0.029) |  71.935 (106.868) |  2.751 (1.263)
|  |  PLS${}_{\text{selected}}$ |  0.131 (0.035) |  70.139 (109.836) |  2.747 (1.201)

Moreover, we compared both methods in terms of their computing times. The
computing time of the LQ method is mainly based on the number of one-
dimensional basis functions used to construct the tensor product basis
functions to estimate the functions of quadratic/interaction components. Our
records showed that, compared with the LQ method using the default number of
basis functions (40 $B$-spline and $20\times 20=400$ tensor product $B$-spline
basis functions), the proposed method requires a similar computing time when a
small number of basis functions are used to estimate the model, while our
proposed method requires more computing time than the LQ method as the number
of basis functions used in the model increases. This is because, for the
proposed method, the dimension of $\boldsymbol{\Psi}$,
$\text{dim}\left(\boldsymbol{\Psi}\right)=K_{\boldsymbol{\mathcal{X}}}^{2}\times
K_{\boldsymbol{\mathcal{X}}}^{2}$, increases exponentially as the number of
basis functions used in the model increases. For example, the proposed method
requires $[1.21,4.18,7.51]$ times more computing times than the LQ method when
$K_{\mathcal{Y}}=K_{\boldsymbol{\mathcal{X}}}=\left[8,12,20\right]$ number of
basis functions are used to estimate the model. In summary, our findings have
demonstrated that the proposed PLS method generally results in improved
accuracy at the expense of longer computational time in comparison with the LQ
method.

We note that our numerical analyses were performed using R 3.6.0 (an example R
code can be found at https://github.com/UfukBeyaztas/FLM_interaction.)

##  4 Empirical data examples

###  4.1 Hawaii ocean data

Available in the R package “FRegSigCom” from Luo and Qi (2018), the Hawaii
ocean dataset consists of five functional variables: salinity, potential
density, temperature, oxygen, and chloropigment. These variables are measured
every 2m between 0 and 200m below the sea surface and are viewed as functions
of depth (see Luo and Qi, 2019, for a more detailed description of the
dataset). The functional variables in this dataset consist of 116 curves,
which are observed at 101 equally spaced points in the interval $[0,200]$. For
this dataset, we considered predicting salinity using the other four
variables. The plots of the functional variables are presented in Figure 3.

Figure 3: Plots of the functional variables: salinity, potential density,
temperature, oxygen, and chloropigment in the Hawaii ocean dataset. The
observations are the functions of depth below the sea surface and $0\leq
s,t\leq 200$. PSS-78 and ITS-90 denote the practical salinity scale of 1978
and the international temperature scale of 1990, respectively.

In the prediction process, the dataset was divided into the following training
and testing samples; the model was constructed based on 90 randomly selected
curves to predict the remaining 26 curves of the salinity variable. This
procedure was repeated 100 times, and for each replication, the performance
metrics (MSPE, RMSPE, and MAPE) were calculated for the main effect, full, and
selected models. In addition, the following two coefficients of determination
were calculated for each replication to measure the goodness of fit ($R^{2}$)
and the relative predictive performances ($R^{2}_{pred}$) of the models (Luo
and Qi, 2019):

$R^{2}=1-\frac{\sum_{i=1}^{90}\left\|\mathcal{Y}^{fit}_{i}(t)-\mathcal{Y}^{train}_{i}(t)\right\|^{2}_{\mathcal{L}_{2}}}{\sum_{i=1}^{90}\left\|\mathcal{Y}^{train}_{i}(t)-\overline{\mathcal{Y}}^{train}(t)\right\|^{2}_{\mathcal{L}_{2}}},\qquad
R^{2}_{pred}=1-\frac{\sum_{i=1}^{26}\left\|\mathcal{Y}^{pred}_{i}(t)-\mathcal{Y}^{test}_{i}(t)\right\|^{2}_{\mathcal{L}_{2}}}{\sum_{i=1}^{26}\left\|\mathcal{Y}^{test}_{i}(t)-\overline{\mathcal{Y}}^{test}(t)\right\|^{2}_{\mathcal{L}_{2}}},$

where $\left\\{\mathcal{Y}^{train}_{i}(t),\mathcal{Y}^{fit}_{i}(t)\right\\}$,
for $i=1,\cdots,90$, are respectively, the $i^{\text{th}}$ observed and fitted
curves of the salinity variable in the training sample, and
$\left\\{\mathcal{Y}^{test}_{i}(t),\mathcal{Y}^{pred}_{i}(t)\right\\}$, for
$i=1,\cdots,26$, are the $i^{\text{th}}$ observed and predicted curves of the
salinity variable in the testing sample, respectively. For this dataset,
$K_{\mathcal{Y}}=K_{\boldsymbol{\mathcal{X}}}=20$ number of basis functions
were used to calculate the performance metrics for the proposed PLS method.

Our findings are reported in Table 2. The results show that the selected
models of both LQ and the proposed method have better performances than their
full and main effect models. The coefficients of determination values produced
by both methods are very close to 1, and these results demonstrate that both
methods have good prediction performances. However, as shown in Table 2, the
selected model of the proposed PLS method has an improved prediction
performance compared with the selected model of LQ.

Table 2: Computed average MSPE, RMSPE, MAPE, $R^{2}$, and $R^{2}_{pred}$ values of the LQ and PLS methods for the Hawaii ocean data. The subscripts ${}_{\text{main}}$, ${}_{\text{full}}$, and ${}_{\text{selected}}$ respectively correspond to the case where the model is estimated using the main effect, full, and selected models. The values given in brackets are the estimated standard errors for the calculated performance metrics. Method |  MSPE |  RMSPE |  MAPE |  $R^{2}$ |  $R^{2}_{pred}$
---|---|---|---|---|---
LQ${}_{\text{main}}$ |  2.6 $\times 10^{-4}$ |  4.6 $\times 10^{-4}$ |  3.3 $\times 10^{-4}$ |  0.994 |  0.999
|  (7.1 $\times 10^{-5}$) |  (5.8 $\times 10^{-5}$) |  (3.7 $\times 10^{-5}$) |  (1.3 $\times 10^{-3}$) |  (4.2 $\times 10^{-5}$)
LQ${}_{\text{full}}$ |  2.1 $\times 10^{-4}$ |  4.1 $\times 10^{-4}$ |  2.8 $\times 10^{-4}$ |  0.994 |  0.999
|  (7.5 $\times 10^{-5}$) |  (6.7 $\times 10^{-5}$) |  (4.9 $\times 10^{-5}$) |  (1.2 $\times 10^{-3}$) |  (3.1 $\times 10^{-5}$)
LQ${}_{\text{selected}}$ |  2.0 $\times 10^{-4}$ |  3.9 $\times 10^{-4}$ |  2.7 $\times 10^{-4}$ |  0.994 |  0.999
|  (7.5 $\times 10^{-5}$) |  (6.8 $\times 10^{-5}$) |  (4.6 $\times 10^{-5}$) |  (1.2 $\times 10^{-3}$) |  (3.1 $\times 10^{-5}$)
PLS${}_{\text{main}}$ |  4.6 $\times 10^{-4}$ |  5.8 $\times 10^{-4}$ |  4.3 $\times 10^{-4}$ |  0.991 |  0.990
|  (2.1 $\times 10^{-4}$) |  (1.2 $\times 10^{-4}$) |  (7.0 $\times 10^{-5}$) |  (2.6 $\times 10^{-5}$) |  (3.1 $\times 10^{-4}$)
PLS${}_{\text{full}}$ |  4.3 $\times 10^{-4}$ |  5.8 $\times 10^{-4}$ |  4.0$\times 10^{-4}$ |  0.993 |  0.993
|  (2.0 $\times 10^{-4}$) |  (1.2 $\times 10^{-4}$) |  (4.7 $\times 10^{-5}$) |  (2.6 $\times 10^{-5}$) |  (5.5 $\times 10^{-4}$)
PLS${}_{\text{selected}}$ |  1.9 $\times 10^{-5}$ |  1.2 $\times 10^{-4}$ |  8.6 $\times 10^{-5}$ |  0.995 |  0.999
|  (5.5 $\times 10^{-6}$) |  (1.7 $\times 10^{-5}$) |  (8.8 $\times 10^{-6}$) |  (2.2 $\times 10^{-5}$) |  (1.8 $\times 10^{-5}$)

To determine the significant main and quadratic/interaction effects terms for
this dataset, a regression model was constructed based on all 116 curves.
According to our results, the potential density $\mathcal{X}_{1}(s)$ and
temperature $\mathcal{X}_{2}(s)$ were selected as the significant main effects
by the forward selection procedures of both the LQ and PLS methods. For the
quadratic and interaction effects terms,
$\mathcal{X}_{1}(s)\mathcal{X}_{1}(r)$,
$\mathcal{X}_{1}(s)\mathcal{X}_{2}(r)$, and
$\mathcal{X}_{2}(s)\mathcal{X}_{2}(r)$ were selected as significant by LQ,
whereas $\mathcal{X}_{1}(s)\mathcal{X}_{1}(r)$ and
$\mathcal{X}_{1}(s)\mathcal{X}_{2}(r)$ were selected as significant by our
proposed method.

We show the estimated coefficient functions of the linear $\beta_{1}(s,t)$ and
$\beta_{2}(s,t)$ and quadratic/interaction ($\gamma_{11}(s,r,t)$ and
$\gamma_{12}(s,r,t)$ when $t=0$, only) terms by our proposed method in Figure
4. The plots of $\beta_{1}(s,t)$ and $\beta_{2}(s,t)$ show that both potential
density and temperature have greater impacts on salinity when $s<25$ or
$s>175$. From the plots of $\gamma_{11}(s,r,t)$ and $\gamma_{12}(s,r,t)$, both
potential density and temperature have negative quadratic/interaction effects
on salinity at $t=0$.

Figure 4: Estimated coefficient functions for the linear ($\beta_{i}(s,t)$ for
$i=1,2$), quadratic ($\gamma_{11}(s,r,t)$, at $t=0$, only), and interaction
($\gamma_{12}(s,r,t)$, at $t=0$, only) terms using the proposed PLS method for
the Hawaii ocean data.

###  4.2 North Dakota weather data

The North Dakota weather data (dataset is available from the North Dakota
Agricultural Weather Network Center: https://ndawn.ndsu.nodak.edu) consists of
five functional variables: solar radiation (MJ/m2), temperature (∘C), bare
soil temperature (∘C), wind speed (m/sec), and wind chill temperature (∘C).
These variables are calculated by averaging daily values for the defined 7-day
periods. The dataset were collected from 91 stations across North Dakota (see
Table 4 in Appendix) from January 1, 2018, to December 31, 2018. For this
dataset, the variables are viewed as functions of weeks, and there are 91
curves observed at 53 equally spaced points in the interval $[1,53]$ for each
variable. Herein, we considered predicting solar radiation using the other
four variables. The plots of the functional variables are presented in Figure
5.

Figure 5: Plots of the functional variables; solar radiation, temperature,
bare soil temperature, wind speed, and wind chill in the North Dakota weather
dataset. The observations are the functions of weeks and $1\leq s,t\leq 53$.

A similar procedure as that for the analysis of the Hawaii ocean data was
followed to calculate the performance metrics and coefficients of
determination values for both the LQ and proposed PLS methods. The dataset was
randomly split into the following two parts: the model was constructed based
on 50 randomly selected curves to predict remaining 41 curves for the solar
radiation variable, and this process was repeated 100 times. For the proposed
method, the performance metrics were calculated based on
$K_{\mathcal{Y}}=K_{\boldsymbol{\mathcal{X}}}=40$ number of basis functions.
The results are presented in Table 3. This table shows that the models
obtained by the forward selection procedures of both the LQ and proposed
methods outperform their full and main effect models. The coefficient of
determination values produced by both methods are very close to 1, which
indicates that both methods have good prediction performances. However, the
results, reported in Table 3, show that the selected model of the proposed PLS
method has better prediction performances compared with the selected model of
LQ.

Table 3: Computed average MSPE, RMSPE, MAPE, $R^{2}$, and $R^{2}_{pred}$ values of the LQ and PLS methods for the North Dakota weather data. The subscripts ${}_{\text{main}}$, ${}_{\text{full}}$, and ${}_{\text{selected}}$ respectively correspond to the case where the model is estimated using the main effect, full, and selected models. The values given in brackets are the estimated standard errors for the calculated performance metrics. Method | MSPE | RMSPE | MAPE | $R^{2}$ | $R^{2}_{pred}$
---|---|---|---|---|---
LQ${}_{\text{main}}$ | 1.445 (0.080) | 0.094 (0.003) | 0.069 (0.002) |  0.988 (9.8 $\times 10^{-4}$) |  0.986 (7.8 $\times 10^{-4}$)
LQ${}_{\text{full}}$ | 1.422 (0.080) | 0.093 (0.003) | 0.068 (0.002) |  0.989 (8.1 $\times 10^{-4}$) |  0.986 (7.5 $\times 10^{-4}$)
LQ${}_{\text{selected}}$ | 1.389 (0.103) | 0.092 (0.004) | 0.067 (0.002) |  0.990 (8.6 $\times 10^{-4}$) |  0.987 (9.7 $\times 10^{-4}$)
PLS${}_{\text{main}}$ | 1.343 (0.078) | 0.090 (0.003) | 0.066 (0.002) |  0.989 (1.1 $\times 10^{-3}$) |  0.987 (7.6 $\times 10^{-4}$)
PLS${}_{\text{full}}$ | 1.352 (0.078) | 0.090 (0.003) | 0.066 (0.002) |  0.989 (1.1 $\times 10^{-3}$) |  0.987 (7.8 $\times 10^{-4}$)
PLS${}_{\text{selected}}$ | 1.300 (0.093) | 0.087 (0.003) | 0.064 (0.002) |  0.990 (9.0 $\times 10^{-4}$) |  0.989 (8.9 $\times 10^{-4}$)

For this dataset, a model using all 91 observations was constructed to
determine the significant main and quadratic/interaction effect terms. Only
the temperature $\mathcal{X}_{1}(s)$ and wind chill temperature
$\mathcal{X}_{4}(s)$ were selected into the final model by both methods. The
quadratic effect $\mathcal{X}_{1}(s)\mathcal{X}_{1}(r)$ was selected as
significant by the LQ while the interaction effect
$\mathcal{X}_{1}(s)\mathcal{X}_{4}(r)$ and the interaction effect
$\mathcal{X}_{4}(s)\mathcal{X}_{4}(r)$ were selected as significant by the
proposed PLS method. The estimated coefficient functions for the main
$\beta(s,t)$ and quadratic/interaction ($\gamma(s,r,t)$ at $t=1$, only)
effects of the temperature and wind chill temperature are presented in Figure
6. It is clear from these plots that temperature and wind chill temperature
have a greater effect on solar radiation when $s>20$, and their quadratic and
interactions have positive and negative quadratic effects along with $s$ and
$r$ at $t=0$.

Figure 6: Estimated coefficient functions for the linear ($\beta_{i}(s,t)$ for
$i=1,2$), quadratic ($\gamma_{11}(s,r,t)$, at $t=0$, only), and interaction
($\gamma_{12}(s,r,t)$, at $t=0$, only) terms using the proposed PLS method for
the North Dakota weather data.

##  5 Conclusion

We propose a PLS method to estimate the function-on-function regression model
where the multiple functional predictors have a quadratic term and interaction
effects. In the proposed method, all the functional predictors and their
quadratic/interaction effects are stacked into a vector, and the functional
PLS regression of the functional response on the vectors of the functional
predictors are considered. Since the direct estimation of a functional
regression model is generally an ill-posed problem, the proposed method is
estimated using the $B$-spline basis expansion of the functional objects.
Basis expansion reduces the infinite-dimensional estimation problem to a
simple finite-dimensional PLS regression setting using particular metrics in
the spaces of the expansion coefficients. Using the metrics associated with
the basis functions, the PLS regression constructed in this way is equivalent
to the PLS regression constructed using a functional response and functional
predictors. The finite sample performances of the proposed method are
evaluated via several Monte Carlo experiments and two empirical data analyses.
The performance of our method are compared with those of an existing method.
Our results have shown that the proposed method produces improved prediction
accuracy than those of the existing method. In addition, we propose a forward
selection procedure to determine the significant main and
quadratic/interaction effects. Our numerical results have proven that the
forward selection procedure further improves the prediction performances of
the PLS method when the true form of the model is even unspecified.

In the proposed method, we consider only $B$-spline basis expansion to
approximate the functional objects. However, the specific forms of bivariate
basis functions (probably tensor products of univariate B-splines) are
generally unknown. In smoothing 2-dimensional surfaces, tensor-product type
splines may ignore interior holes (or gaps) and yield extra bias as well as
oversmoothing. Alternatively, the basis functions of non-tensor-product type
for bivariate smoothing, such as Green’s function (Wood et al., 2008) and
bivariate B-splines (Pan, 2011) may be used.

## Appendix

### Identifiability of the proposed model

Model identifiability, which is an important theme in the estimation phase,
for the function-on-function regression model with only main effects given in
(1.1) was theoretically discussed by He et al. (2000) and Chiou et al. (2004).
Later, Scheipl and Greven (2016) discussed the identifiability of such models
with realistic applications. As is stated in the Proposition 3.1 of Scheipl
and Greven (2016), the coefficient function $\beta_{m}(s,t)$ in (1.1) is
identifiable if and only if the kernel space of the covariance operator of the
functional predictors ($K^{\mathcal{X}}$) is empty except the zero function,
i.e., $ke(K^{\mathcal{X}})=\\{0\\}$. Luo and Qi (2019) extended the
theoretical condition of model identifiability for function-on-function
regression with only main effects to the model with interaction and quadratic
effects. They defined the following $\\{p(p+3)/2\\}\times\\{p(p+3)/2\\}$
dimensional covariance matrix of the predictor functions:

(5.9)
$\boldsymbol{\Sigma}=\begin{bmatrix}\boldsymbol{A}(s,r)&\boldsymbol{B}(s,r,u)\\\
\boldsymbol{B}^{\top}(s,r,u)&\boldsymbol{C}(s,r,s^{\prime},r^{\prime})\end{bmatrix},$

where the $p\times p$ and $p\times\\{p(p+1)/2\\}$ dimensional submatrices
$\boldsymbol{A}(s,r)$ and $\boldsymbol{B}(s,r,u)$ equal to
$\text{Cov}\left(\mathcal{X}(s)\mathcal{X}(r)\right)$ and
$\text{Cov}\left(\mathcal{X}(s),\mathcal{X}(r)\mathcal{X}(u)\right)$,
respectively. The $\\{p(p+1)/2\\}\times\\{p(p+1)/2\\}$ dimensional matrix
$\boldsymbol{C}(s,r,s^{\prime},r^{\prime})$ equals to
$\text{Cov}\left(\mathcal{X}(s)\mathcal{X}(r),\mathcal{X}(s^{\prime})\mathcal{X}(r^{\prime})\right)$.
Then, Luo and Qi (2019) proved that the coefficients in model (1.3) are
identifiable if the kernel space of $\boldsymbol{\Sigma}$ in (5.9) only has
the zero function. Please see Proposition S.1. of Luo and Qi (2019) for the
proof. Based on the above, our proposed method given in (2.4) is also
identifiable since the proposed method is a reformulated version of the model
(1.3).

### Proof of Proposition (2.2)

###### Proof.

The proof of Proposition (2.2) is an extended version of the proof of
Proposition 2 of Aguilera et al. (2010). Denote by
$\Lambda=\boldsymbol{\Phi}^{1/2}\boldsymbol{c}$ and
$\Pi=\boldsymbol{\Psi}^{1/2}\boldsymbol{d}$ the column random vectors.

For any random variable $Z\in\mathcal{L}_{2}[0,1]$, the followings hold when
$h=1$:

$\displaystyle W^{\mathcal{Y}}$
$\displaystyle=\int_{0}^{1}\sum_{k=1}^{K_{\mathcal{Y}}}c_{k}\phi_{k}(t)\mathbb{E}\left[\sum_{i=1}^{K_{\mathcal{Y}}}Zc_{i}\phi_{i}(t)\right]dt=\Lambda^{\top}\mathbb{E}\left[\Lambda
Z\right]=W^{\Lambda}Z,$ $\displaystyle W^{\boldsymbol{\mathcal{X}}}$
$\displaystyle=\int_{0}^{1}\int_{0}^{1}\sum_{j=1}^{K_{\boldsymbol{\mathcal{X}}}}\sum_{l=1}^{K_{\boldsymbol{\mathcal{X}}}}d_{jl}\psi_{jl}(s,r)\mathbb{E}\left[\sum_{\iota=1}^{K_{\boldsymbol{\mathcal{X}}}}\sum_{\tau=1}^{K_{\boldsymbol{\mathcal{X}}}}Zd_{\iota\tau}\psi_{\iota\tau}(s,r)\right]dsdr=\Pi^{\top}\mathbb{E}\left[\Pi
Z\right]=W^{\Pi}Z.$

The residuals obtained from the first iteration are given by:

$\begin{cases}\mathcal{Y}_{1}(t)=\mathcal{Y}(t)-\zeta_{1}(t)\eta_{1},&\zeta_{1}(t)\in\mathcal{L}_{2}[0,1]\\\
\Lambda_{1}=\Lambda-\widetilde{\zeta}_{1}\eta_{1},&\widetilde{\zeta_{1}}\in\mathbb{R}^{K_{\mathcal{Y}}}\end{cases}$

$\begin{cases}\boldsymbol{\mathcal{X}}_{1}(s,r)=\boldsymbol{\mathcal{X}}(s,r)-p_{1}(s,r)\eta_{1},&p_{1}(s,r)\in\mathcal{L}_{2}[0,1]\\\
\Pi_{1}=\Pi-\widetilde{p}_{1}\eta_{1},&\widetilde{p}_{1}\in\mathbb{R}^{K_{\boldsymbol{\mathcal{X}}}^{2}}\end{cases}$

Then, we have $\Lambda_{1}=\boldsymbol{\Phi}^{1/2}\boldsymbol{c}_{1}$ and
$\Pi_{1}=\boldsymbol{\Psi}^{1/2}\boldsymbol{d}_{1}$, where
$\boldsymbol{c}_{1}$ and $\boldsymbol{d}_{1}$ are the random vectors of the
basis coefficients of $\mathcal{Y}_{1}(t)$ and
$\boldsymbol{\mathcal{X}}_{1}(s,r)$, respectively. For functional PLS, the
residuals calculated from the first iteration are as follows:

$\displaystyle\mathcal{Y}_{1}(t)$
$\displaystyle=\sum_{k=1}^{K_{\mathcal{Y}}}c_{1,k}\phi_{k}(t)=\boldsymbol{c}_{1}^{\top}\boldsymbol{\Phi}(t),\qquad\forall
t\in\mathcal{L}_{2}[0,1],$ $\displaystyle\boldsymbol{\mathcal{X}}_{1}(s,r)$
$\displaystyle=\sum_{j=1}^{K_{\boldsymbol{\mathcal{X}}}}\sum_{l=1}^{K_{\boldsymbol{\mathcal{X}}}}d_{1,jl}\psi_{jl}(s,r)=\boldsymbol{d}_{1}^{\top}\boldsymbol{\Psi}(s,r),\qquad\forall
s,r\in\mathcal{L}_{2}[0,1].$

On the other hand, the residuals from the PLS regression are obtained as
follows:

$\displaystyle\mathcal{Y}_{1}(t)$
$\displaystyle=\boldsymbol{\Phi}^{\top}(t)\boldsymbol{c}-\eta_{1}\frac{\mathbb{E}\left[\boldsymbol{\Phi}^{\top}(t)\boldsymbol{c}\eta_{1}\right]}{\mathbb{E}\left[\eta_{1}^{2}\right]},$
$\displaystyle\boldsymbol{\mathcal{X}}_{1}(s,r)$
$\displaystyle=\boldsymbol{\Psi}^{\top}(s,r)\boldsymbol{d}-\eta_{1}\frac{\mathbb{E}\left[\boldsymbol{\Psi}^{\top}(s,r)\boldsymbol{d}\eta_{1}\right]}{\mathbb{E}\left[\eta_{1}^{2}\right]},$
$\displaystyle\Lambda_{1}$
$\displaystyle=\boldsymbol{\Phi}^{1/2}\boldsymbol{c}-\eta_{1}\frac{\mathbb{E}\left[\boldsymbol{\Phi}^{1/2}\boldsymbol{c}\eta_{1}\right]}{\mathbb{E}\left[\eta_{1}^{2}\right]},$
$\displaystyle\Pi_{1}$
$\displaystyle=\boldsymbol{\Psi}^{1/2}\boldsymbol{d}-\eta_{1}\frac{\mathbb{E}\left[\boldsymbol{\Psi}^{1/2}\boldsymbol{d}\eta_{1}\right]}{\mathbb{E}\left[\eta_{1}^{2}\right]}$

and thus we have
$\boldsymbol{c}_{1}=\boldsymbol{c}-\eta_{1}\frac{\mathbb{E}\left(\boldsymbol{c}\eta_{1}\right)}{\mathbb{E}\left[\eta_{1}^{2}\right]}$
and
$\boldsymbol{d}_{1}=\boldsymbol{d}-\eta_{1}\frac{\mathbb{E}\left[\boldsymbol{d}\eta_{1}\right]}{\mathbb{E}\left[\eta_{1}^{2}\right]}$,
that is, $W^{\mathcal{Y}_{1}}=W^{\Lambda_{1}}$ and
$W^{\boldsymbol{\mathcal{X}}_{1}}=W^{\Pi_{1}}$. Now, assuming
$W^{\mathcal{Y}_{\ell}}=W^{\Lambda_{\ell}}$ and
$W^{\boldsymbol{\mathcal{X}}_{\ell}}=W^{\Pi_{\ell}}$ for each $\ell\leq h$.
Also, let us assume that $W^{\mathcal{Y}_{h+1}}=W^{\Lambda_{h+1}}$ and
$W^{\boldsymbol{\mathcal{X}}_{h+1}}=W^{\Pi_{h+1}}$. Then, at step $h+1$ the
followings hold:

$\begin{cases}\mathcal{Y}_{h+1}(t)=\mathcal{Y}(t)-\zeta_{1}(t)\eta_{1}-\cdots-\zeta_{h}(t)\eta_{h},&\zeta_{i}(t)\in\mathcal{L}_{2}[0,1]\\\
\Lambda_{h+1}=\Lambda-\widetilde{\zeta}_{1}\eta_{1}-\cdots-\widetilde{\zeta}_{h}\eta_{h},&\widetilde{\zeta_{i}}\in\mathbb{R}^{K_{\mathcal{Y}}}\end{cases}$

$\begin{cases}\boldsymbol{\mathcal{X}}_{h+1}(s,r)=\boldsymbol{\mathcal{X}}(s,r)-p_{1}(s,r)\eta_{1}-\cdots-
p_{h}(s,r)\eta_{h},&p_{i}(s,r)\in\mathcal{L}_{2}[0,1]\\\
\Pi_{h+1}=\Pi-\widetilde{p}_{1}\eta_{1}-\cdots-\widetilde{p}_{h}\eta_{h},&\widetilde{p}_{i}\in\mathbb{R}^{K_{\boldsymbol{\mathcal{X}}}^{2}}\end{cases}$

As for $h=1$ and using the orthogonality of $\eta_{i}$ for $i=1,\cdots,h$, it
is shown that $\Lambda_{h+1}=\boldsymbol{\Phi}^{1/2}\boldsymbol{c}_{h+1}$ and
$\Pi_{h+1}=\boldsymbol{\Psi}^{1/2}\boldsymbol{d}_{h+1}$, which concludes the
proof. ∎

### Station names for the North Dakota weather data

Table 4: Station names for the North Dakota weather data. Station | Station | Station | Station | Station | Station | Station
---|---|---|---|---|---|---
Ada | Carson | Eldred | Harvey | Linton | Pekin | St. Thomas
Alamo | Cavalier | Fargo | Hazen | Lisbon | Perley | Stephen
Baker | Cooperstown | Fingal | Hettinger | Mandan | Pillsbury | Streeter
Beach | Crary | Finley | Hillsboro | Marion | Plaza | Tappen
Berthold | Crosby | Forest River | Hofflund | Mavie | Prosper | Turtle Lake
Bottineau | Dagmar | Fort Yates | Hope | Mayville | Redstone | Ulen
Bowbells | Dazey | Fox | Humboldt | McHenry | Robinson | Wahpeton
Bowman | Dickinson | Froid | Inkster | Michigan | Rolla | Warren
Brampton | Dooley | Galesburg | Jamestown | Minot | Roseau | Watford City
Brorson | Dunn | Garrison | Karlsruhe | Mohall | Ross | Waukon
Campbell | Edgeley | Grafton | Kennedy | Mooreton | Rugby | Williams
Cando | Edmore | Grand Forks | Langdon | Mott | Sabin | Williston
Carrington | Ekre | Greenbush | Leonard | Oakes | Sidney | Wishek

## References

  * (1)
  * Aguilera et al. (2010) Aguilera, A. M., Escabias, M., Preda, C. and Saporta, G. (2010), ‘Using basis expansions for estimating functional PLS regression: Applications with chemometric data’, Chemometrics and Intelligent Laboratory Systems 104(2), 289–305.
  * Aguilera et al. (2016) Aguilera, A. M., Escabias, M., Preda, C. and Saporta, G. (2016), ‘Penalized versions of functional PLS regression’, Chemometrics and Intelligent Laboratory Systems 154, 80–92.
  * Aguilera et al. (1999) Aguilera, A. M., Ocana, F. A. and Valderrama, M. J. (1999), ‘Forecasting with unequally spaced data by a functional principal component approach’, Test 8(1), 233–254.
  * Beyaztas and Shang (2020) Beyaztas, U. and Shang, H. L. (2020), ‘On function-on-function regression: Partial least squares approach’, Environmental and Ecological Statistics 27(1), 95–114.
  * Chiou et al. (2004) Chiou, J.-M., Müller, H. G. and Wang, J. L. (2004), ‘Functional response models’, Statistica Sinica 14, 675–693.
  * Chiou et al. (2016) Chiou, J.-M., Yang, Y.-F. and Chen, Y.-T. (2016), ‘Multivariate functional linear regression and prediction’, Journal of Multivariate Analysis 146, 301–312.
  * Cuevas (2014) Cuevas, A. (2014), ‘A partial overview of the theory of statistics with functional data’, Journal of Statistical Planning and Inference 147, 1–23.
  * Dayal and MacGregor (1997) Dayal, B. S. and MacGregor, J. F. (1997), ‘Improved PLS algorithms’, Journal of Chemometrics 11(1), 73–85.
  * de Jong (1993) de Jong, S. (1993), ‘SIMPLS: An alternative approach to partial least squares regression’, Chemometrics and Intelligent Laboratory Systems 18(3), 251–263.
  * de Lathuauwer et al. (2000) de Lathuauwer, L., de Moor, B. and Vandewalle, J. (2000), ‘A multilinear singular value decomposition’, SIAM Journal on Matrix Analysis and Applications 21(4), 1253–1278.
  * Delaigle and Hall (2012a) Delaigle, A. and Hall, P. (2012a), ‘Achieving near perfect classification for functional data’, Journal of the Royal Statistical Society Series B 74(2), 267–286.
  * Delaigle and Hall (2012b) Delaigle, A. and Hall, P. (2012b), ‘Methodology and theory for partial least squares applied to functional data’, The Annals of Statistics 40(1), 322–352.
  * Escabias et al. (2007) Escabias, M., Aguilera, A. M. and Valderrama, M. J. (2007), ‘Functional PLS logit regression model’, Computational Statistics & Data Analysis 51(10), 4891–4902.
  * Escoufier (1970) Escoufier, Y. (1970), ‘Echantillonnage dans une population de variables aléatories réelles’, Publications de I’Institut de Statistique de I’Université de Paris 19(4), 1–47.
  * Febrero-Bande et al. (2017) Febrero-Bande, M., Galeano, P. and Gozalez-Manteiga, W. (2017), ‘Functional principal component regression and functional partial least-squares regression: An overview and a comparative study’, International Statistical Review 85(1), 61–83.
  * Ferraty and Vieu (2006) Ferraty, F. and Vieu, P. (2006), Nonparametric Functional Data Analysis, Springer, New York.
  * Fuchs et al. (2015) Fuchs, K., Scheipl, F. and Greven, S. (2015), ‘Penalized scalar-on-functions regression with interaction term’, Computational Statistics & Data Analysis 81, 38–51.
  * He et al. (2000) He, G., Müller, H. G. and Wang, J. L. (2000), Extending correlation and regression from multivariate to functional data, in M. L. Puri, ed., ‘Asymptotics in Statistics and Probability’, VSP International Science Publishers, pp. 197–210.
  * He et al. (2010) He, G., Müller, H. G., Wang, J. L. and Yang, W. (2010), ‘Functional linear regression via canonical analysis’, Bernoulli 16(3), 705–729.
  * Horvath and Kokoszka (2012) Horvath, L. and Kokoszka, P. (2012), Inference for Functional Data with Applications, Springer, New York.
  * Hsing and Eubank (2015) Hsing, T. and Eubank, R. (2015), Theoretical Foundations of Functional Data Analysis, with an Introduction to Linear Operators, John Wiley & Sons, Chennai, India.
  * Hyndman and Shang (2009) Hyndman, R. J. and Shang, H. L. (2009), ‘Forecasting functional time series’, Journal of the Korean Statistical Society 38(3), 199–211.
  * Ivanescu et al. (2015) Ivanescu, A. E., Staicu, A.-M., Scheipl, F. and Greven, S. (2015), ‘Penalized function-on-function regression’, Computational Statistics 30(2), 539–568.
  * Kokoszka and Reimherr (2017) Kokoszka, P. and Reimherr, M. (2017), Introduction to Functional Data Analysis, CRC Press, Boca Raton.
  * Krämer et al. (2008) Krämer, N., Boulesteix, A.-L. and Tutz, G. (2008), ‘Penalized partial least squares with applications to B-spline transformations and functional data’, Chemometrics and Intelligent Laboratory Systems 94(1), 60–69.
  * Luo and Qi (2018) Luo, R. and Qi, X. (2018), FRegSigCom: Functional Regression using Signal Compression Approach.  R package version 0.3.0.
https://CRAN.R-project.org/package=FRegSigCom

  * Luo and Qi (2019) Luo, R. and Qi, X. (2019), ‘Interaction model and model selection for function-on-function regression’, Journal of Computational and Graphical Statistics 28(2), 309–322.
  * Matsui (2020) Matsui, H. (2020), ‘Quadratic regression for functional response models’, Econometrics and Statistics 13, 125–136.
  * Matsui et al. (2009) Matsui, H., Kawano, S. and Konishi, S. (2009), ‘Regularized functional regression modeling for functional response and predictors’, Journal of Math-for-Industry 1(A3), 17–25.
  * Müller and Yao (2008) Müller, H.-G. and Yao, F. (2008), ‘Functional additive models’, Journal of the American Statistical Association 103(484), 1534–1544.
  * Pan (2011) Pan, H. (2011), Bivariate B-splines and its Applications in Spatial Data Analysis, PhD thesis, Texas A & M University.
  * Preda and Saporta (2005) Preda, C. and Saporta, G. (2005), ‘PLS regression on a stochastic process’, Computational Statistics & Data Analysis 48(1), 149–158.
  * Preda et al. (2007) Preda, C., Saporta, G. and Lévéder, C. (2007), ‘PLS classification of functional data’, Computational Statistics 22(2), 223–235.
  * Ramsay and Dalzell (1991) Ramsay, J. O. and Dalzell, C. J. (1991), ‘Some tools for functional data analysis’, Journal of the Royal Statistical Society, Series B 53(3), 539–572.
  * Ramsay and Silverman (2002) Ramsay, J. O. and Silverman, B. W. (2002), Applied Functional Data Analysis, Springer, New York.
  * Ramsay and Silverman (2006) Ramsay, J. O. and Silverman, B. W. (2006), Functional Data Analysis, Springer, New York.
  * Reiss and Ogden (2007) Reiss, P. T. and Ogden, T. R. (2007), ‘Functional principal component regression and functional partial least squares’, Journal of the American Statistical Association 102(479), 984–996.
  * Scheipl and Greven (2016) Scheipl, F. and Greven, S. (2016), ‘Identifiability in penalized function-on-function regression models’, Electronic Journal of Statistics 10, 495–526.
  * Sun and Wang (2020) Sun, Y. and Wang, Q. (2020), ‘Function-on-function quadratic regression models’, Computational Statistics & Data Analysis 142, Article number 106814.
  * Tucker (1938) Tucker, R. S. (1938), ‘The reasons for price rigidity’, The American Economic Review 28(1), 41–54.
  * Usset et al. (2016) Usset, J., Staicu, A.-M. and Maity, A. (2016), ‘Interaction models for functional regression’, Computational Statistics & Data Analysis 94, 317–329.
  * Wang (2014) Wang, W. (2014), ‘Linear mixed function‐on‐function regression models’, Biometrics 70(4), 794–801.
  * Wold (1974) Wold, H. (1974), ‘Causal flows with latent variables: Partings of the ways in the light of NIPALS modelling’, European Economic Review 5(1), 67–86.
  * Wood et al. (2008) Wood, S. N., Bravington, M. V. and Hedley, S. L. (2008), ‘Soap film smoothing’, Journal of the Royal Statistical Society, Series B 70(Part 5), 931–955.
  * Yao and Müller (2010) Yao, F. and Müller, H.-G. (2010), ‘Functional quadratic regression’, Biometrika 97(1), 49–64.
  * Yao et al. (2005) Yao, F., Müller, H.-G. and Wang, J.-L. (2005), ‘Functional linear regression analysis for longitudinal data’, The Annals of Statistics 33(6), 2873–2903.