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

Ufuk Beyaztas
Department of Statistics
Marmara University

Han Lin Shang
Department of Actuarial Studies and Business Analytics
Macquarie University Postal address: Department of Statistics, Marmara
Üniversitesi Göztepe Yerleşkesi, 34722, Kadiköy-Istanbul, Turkey; Email:
<EMAIL_ADDRESS>0000-0002-5208-4950 0000-0003-1769-6430

###### Abstract

The function-on-function linear regression model in which the response and
predictors consist of random curves has become a general framework to
investigate the relationship between the functional response and functional
predictors. Existing methods to estimate the model parameters may be sensitive
to outlying observations, common in empirical applications. In addition, these
methods may be severely affected by such observations, leading to undesirable
estimation and prediction results. A robust estimation method, based on
iteratively reweighted simple partial least squares, is introduced to improve
the prediction accuracy of the function-on-function linear regression model in
the presence of outliers. The performance of the proposed method is based on
the number of partial least squares components used to estimate the function-
on-function linear regression model. Thus, the optimum number of components is
determined via a data-driven error criterion. The finite-sample performance of
the proposed method is investigated via several Monte Carlo experiments and an
empirical data analysis. In addition, a nonparametric bootstrap method is
applied to construct pointwise prediction intervals for the response function.
The results are compared with some of the existing methods to illustrate the
improvement potentially gained by the proposed method.

Keywords: Basis function; dimension reduction; SIMPLS; weighted likelihood

##  1 Introduction

Over the last decade, the availability of functional data whose sample
elements are in the form of curves, images, shapes, or more general manifold-
valued objects sampled over a continuum measure has progressively increased in
many scientific fields. Consequently, the interest in and need for tools that
tackle the analysis of functional datasets is increasing significantly. For
example, consult Ramsay & Silverman (2002, 2006), Ferraty & Vieu (2006), Rao
(2010), Horvath & Kokoszka (2012), Cuevas (2014), Hsing & Eubank (2015),
Kokoszka & Reimherr (2017), Younse (2020), and Martines-Hernandez & Genton
(2020) for recent developments on the theory and applications of functional
data analysis tools.

This paper is devoted to studying the function-on-function linear regression
model (FFLRM), where the response and predictors are functional. FFLRM is
among the most interesting functional data analysis tools and has received
much attention in the statistical literature. This model intends to explore
the association between a functional response and one or more functional
predictors. In this context, several FFLRMs have been proposed (see, e.g.,
Ramsay & Dalzell, 1991; Yao et al., 2005; Ramsay & Silverman, 2006; Müller &
Yao, 2008; Matsui et al., 2009; He et al., 2010; Wang, 2014; Ivanescu et al.,
2015; Chiou et al., 2016). In these models, the associations between the
functional variables are characterized via the integrated functional
predictors weighted by the unknown regression parameter function. Estimating
such models is a difficult task since the regression parameter function
belongs to an infinite-dimensional space. The common practical approach to
overcoming this problem is projecting the regression parameter function onto a
finite-dimensional space using a basis expansion technique. For this purpose,
several approaches have been proposed. For example, with the general basis
expansion methods (i.e., $B$-spline, Fourier, and Gaussian basis function),
the traditional estimation methods, such as least-squares (LS) (Ramsay &
Dalzell, 1991; Ramsay & Silverman, 2006), weighted LS (Yamanishi & Tanaka,
2003), penalized maximum likelihood (Matsui et al., 2009),
expectation/conditional maximization either algorithm (Wang, 2014), and
penalized spline (Ivanescu et al., 2015) have been extended to functional data
to estimate the FFLRMs. However, the use of general basis expansion methods
may require a large number of basis functions to project the regression
parameter function, leading to overfitting of the model. In addition, when a
large number of basis functions is used, the aforementioned estimation methods
suffer from the multicollinearity problem (Preda & Saporta, 2005; Aguilera et
al., 2016). Moreover, in such a case, these methods may be computationally
intensive (Beyaztas & Shang, 2020a). Compared with general basis expansion
methods, the dimension reduction techniques such as functional principal
component (FPC) regression (Aguilera et al., 1999; Yao et al., 2005; Hall &
Hosseini-Nasab, 2006; Chiou et al., 2016) and functional partial least squares
(FPLS) regression (Preda & Schiltz, 2011; Beyaztas & Shang, 2020a) solve
multicollinearity problem and decrease the computation burden in a high-
dimensional setting. Thus, the FPC and FPLS are preferable to general basis
expansion methods in the estimation of FFLRMs.

In the FPC and FPLS regression methods, the infinite-dimensional functions of
the FFLRMs are projected onto the finite-dimensional orthonormal FPC and FPLS
bases. In doing so, both methods produce orthogonal latent components. The
regression model constructed using the latent components of the functional
response and functional predictors are then used to approximate the regression
problem of functional response on the functional predictors. The FPC considers
maximizing the covariance between the functional predictors to extract the
latent components. However, a few FPCs generally contain the most information
about the covariance between the functional predictors, which are not
necessarily important for representing FPCs. All or some of the most important
terms accounting for the interaction between the basis functions and
functional predictors might come from later principal components (Delaigle &
Hall, 2012). In contrast, the FPLS components are extracted by maximizing the
covariance between the functional response and functional predictors. Thus,
compared with FPC, the FPLS produces more informative latent components with
fewer terms, which makes FPLS more preferable than FPC (Reiss & Odgen, 2007;
Delaigle & Hall, 2012; Febrero-Bande et al., 2017). Therefore in this paper,
we consider the FPLS method to estimate the FFLRM.

For the functional regression models, various FPLS methods have been proposed
to investigate the relationship between functional variables (see, e.g., Preda
& Saporta, 2005; Hyndman & Shang, 2009; Aguilera et al., 2010, 2016; Preda &
Schiltz, 2011; Beyaztas & Shang, 2020a, 2021). In these methods, the FPLS
basis functions are iteratively extracted from squared covariance between the
functional response and vector of functional predictors to calculate the
latent components. The extracted components are then used to estimate the
underlying FFLRM. However, the direct estimation of the FPLS basis functions
is an ill-posed problem because of the infinite-dimensional structures of the
functional random variables. To overcome this problem, Preda & Schiltz (2011)
and Beyaztas & Shang (2020a) proposed approximating the FPLS regression via
the finite-dimensional PLS regression of the basis expansion coefficients. The
coefficients are obtained by projecting the infinite-dimensional functional
random variables onto a finite-dimensional space via a basis expansion method.
Preda & Schiltz (2011) and Beyaztas & Shang (2021) showed that the same PLS
components are obtained under both the FPLS and the PLS regression constructed
using the basis expansion coefficients. The numerical results provided by the
studies above showed that the FPLS produces an improved predictive performance
over the existing methods because the extracted PLS components consider the
correlation between the response and predictor variables.

Most of the existing methods, including FPLS, used to estimate the FFLRMs are
non-robust to outlying observations, which are generated by a stochastic
process with a distribution different from that of the vast majority of the
remaining observations (see, e.g., Raña et al., 2015). In functional data,
there are three types of outliers:  (1) magnitude outliers, which are points
far from the bulk of the data;  (2) shape outliers, which fall within the
range of the data but differ in shape from the bulk of the data; and  (3) the
combination of magnitude and shape outliers (e.g., see Beyaztas & Shang,
2019).  In this study, we consider the combination of magnitude and shape
outliers. In the case of outliers, the LS or maximum likelihood-based
estimation methods produce biased estimates; thus, predictions obtained from
the fitted model become unreliable (see, e.g., Zhang et al., 2018). To
overcome this issue, several robust methods have been proposed in functional
regression models. In the scalar-on-function regression model (where the
response is a scalar random variable while the predictors consist of random
curves), see for example, Maronna & Yohai (2013), Shin & Lee (2016), and
Kalogridis & Aelst (2019). In the context of FFLRM, Gervini (2012) proposed an
outlier-resistant loss function-based robust method. Also, Hullait et al.
(2021) proposed a Fisher-consistent robust regression model, and they proposed
a robust functional principal component (RFPC)-based estimator to estimate the
model parameter. However, these two methods are based on the simple FFLRM,
with only one predictor in the model. To the best of our knowledge, there is
no method to robustly estimate the FFLRM by extracting information from both
the response and multiple predictors.

This paper proposes a robust FPLS method, which allows for more than one
functional predictor in the model to estimate the FFLRM. Similar to the FPLS
methods of Preda & Schiltz (2011) and Beyaztas & Shang (2020a), in the
proposed method, the infinite-dimensional FPLS regression model is
approximated via the finite-dimensional PLS regression of the basis expansion
coefficients. In the finite-dimensional case, the PLS components are obtained
via several algorithms, such as nonlinear iterative PLS (Wold, 1974), SIMPLS
(de Jong, 1993), and improved kernel PLS (Dayal & MacGregor, 1997). Our
proposed method restricts our attention to SIMPLS, which deflates the
empirical covariance matrix to obtain the PLS components. The traditional
SIMPLS algorithm uses an LS-type estimator to extract the PLS components.
Thus, it may seriously be affected by outlying observations and lead to
unreliable estimation and prediction results. Several robust PLS algorithms
have been developed to deal with outliers. For example, Wakelinc & Macfie
(1992) and Cummins & Andrews (1995) proposed robust PLS algorithms using
robust iteratively reweighted algorithms. However, these methods use non-
robust initial weights and are not robust to outliers in the predictor space
(called leverage points). Using a robust covariance matrix with the SIMPLS
algorithm, Hubert & Branden (2003) proposed a robust PLS regression. Serrneels
et al. (2005) developed a robust PLS algorithm based on the M-regression
estimator. In addition, using the weight function developed by Markatou
(1996), Alin & Agostinelli (2017) proposed a robust iteratively reweighted
SIMPLS method (IRSIMPLS). The latter methods are robust to outliers in the
predictor space (vertical outliers) and leverage points. The numerical
analyses performed by Alin & Agostinelli (2017) demonstrate that the IRSIMPLS
produces improved results compared with its competitors. Therefore, we
consider the IRSIMPLS algorithm in the proposed method to robustly approximate
the infinite-dimensional FPLS regression model in the presence of outliers. In
the numerical analyses given in Sections 3, the predictive performance of the
proposed method is examined via several Monte Carlo experiments, and the
results are compared with those of existing methods. Our numerical results
show that the proposed method is robust to outliers and improves prediction
performance compared to the available methods when outliers are present in the
data. Also, our results show that the proposed robust FPLS approach produces
competitive or even better results with those of the RFPC method of Hullait et
al. (2021).

The remainder of this paper is organized as follows. An overview of the FPLS
approach and the methodology of the proposed method are presented in Section
2. A series of Monte Carlo experiments are conducted to evaluate the finite-
sample performance of the proposed method, and the results are given in
Section 3. The proposed method is used to analyze an empirical data analysis
and the results presented in Section 4. Section 5 concludes the paper, along
with some ideas on how the methodology can be further extended.

##  2 The functional PLS

Let $\mathcal{Y}(t)$ and $\mathcal{X}_{m}(s)$, for $m=1,\ldots,M$,
respectively, denote a response and $M$ sets of predictors given as functions
with closed and bounded intervals $t\in T$ and $s\in S$ as domains. We assume
that both the functional response and the functional predictors are second-
order stochastic processes with finite second-order moments, and are the
elements of the $\mathcal{L}_{2}$ separable Hilbert space, the space of
square-integrable and real-valued functions. Without loss of generality, we
also assume that $\mathcal{Y}(t)$ and $\mathcal{X}_{m}(s)$ for $m=1,\ldots,M$,
are zero-mean stochastic processes so that
$\mathbb{E}\left[\mathcal{Y}(t)\right]=\mathbb{E}\left[\mathcal{X}_{m}(s)\right]=0$.

We consider the following FFLRM for exploring the association between the
functional response and $M$ sets of functional predictors:

$\mathcal{Y}(t)=\int_{S}\boldsymbol{\mathcal{X}}^{\top}(s)\boldsymbol{\beta}(s,t)\textrm{d}s+\epsilon(t),$
(1)

where
$\boldsymbol{\mathcal{X}}(s)=\left[\mathcal{X}_{1}(s),\ldots,\mathcal{X}_{M}(s)\right]^{\top}$
is the matrix-valued functions comprising of $M$ sets of functional
predictors,
$\boldsymbol{\beta}(s,t)=\left[\beta_{1}(s,t),\ldots,\beta_{M}(s,t)\right]^{\top}$
denotes the vector of bivariate coefficient functions, and $\epsilon(t)$ is
the error function which is assumed to be independent of
$\mathcal{X}_{m}(t)$s. The primary aim in (1) is to estimate the regression
parameter function $\boldsymbol{\beta}(s,t)$.

The FPLS components associated with the FFLRM in (1) can be obtained as the
solutions of Tucker’s criterion (Tucker, 1938) extended to functional
variables as follows:

$\underset{\begin{subarray}{c}\kappa\in\mathcal{L}_{2},~{}\|\kappa_{m}\|_{\mathcal{L}_{2}}=1,\forall
m\in[1,\ldots,M]\\\
\zeta\in\mathcal{L}_{2},~{}\|\zeta\|_{\mathcal{L}_{2}}=1\end{subarray}}{\max}\text{Cov}^{2}\left(\int_{S}\boldsymbol{\mathcal{X}}^{\top}(s)\bm{\kappa}(s)ds,~{}\int_{T}\mathcal{Y}(t)\zeta(t)dt\right),$
(2)

where $\bm{\kappa}(s)=[\kappa_{1}(s),\ldots,\kappa_{m}(s)]^{\top}$ and
$\zeta(t)$ denote weight functions associated with the predictor and response
variables, respectively.

As noted by Preda & Saporta (2005), the FPLS components are also the
eigenvectors associated with the largest eigenvalues of the product of
Escoufier operators (Escoufier, 1970). The first FPLS component of the FFLRM
of $\mathcal{Y}(t)$ on $\boldsymbol{\mathcal{X}}(s)$, denoted by
$\bm{\eta}_{1}$, is obtained by conducting the ordinary linear regressions of
$\boldsymbol{\mathcal{X}}(s)$ on $\bm{\eta}_{1}$ and $\mathcal{Y}(t)$ on
$\bm{\eta}_{1}$ as follows:

$\bm{\eta}_{1}=\int_{S}\boldsymbol{\mathcal{X}}^{\top}(s)\bm{\kappa}_{1}(s)\textrm{d}s,$

where $\bm{\kappa}_{1}(s)$, the weight function associated with the first FPLS
component, is defined as follows:

$\bm{\kappa}_{1}(s)=\frac{\int_{T}\mathbb{E}\left[\boldsymbol{\mathcal{X}}(s)\mathcal{Y}(t)\right]\textrm{d}t}{\sqrt{\int_{S}\left(\int_{T}\mathbb{E}\left[\boldsymbol{\mathcal{X}}(s)\mathcal{Y}(t)\right]\textrm{d}t\right)^{2}\textrm{d}s}}.$

FPLS is an iterative method, and thus, the subsequent FPLS components are
computed via an iterative stepwise procedure. Let $h=1,2,\ldots$ denote the
iteration number. In addition, let
$\boldsymbol{\mathcal{X}}_{0}(s)=\boldsymbol{\mathcal{X}}(s)$ and
$\mathcal{Y}_{0}(t)=\mathcal{Y}(t)$. Denote by
$\boldsymbol{\mathcal{X}}_{h}(s)$ and $\mathcal{Y}_{h}(t)$ the residuals of
the following regression models:

$\displaystyle\boldsymbol{\mathcal{X}}_{h}(s)$
$\displaystyle=\boldsymbol{\mathcal{X}}_{h-1}(s)-\bm{p}_{h}(s)\bm{\eta}_{h},$
$\displaystyle\mathcal{Y}_{h}(t)$
$\displaystyle=\mathcal{Y}_{h-1}(t)-\zeta_{h}(t)\bm{\eta}_{h},$

where
$\bm{p}_{h}(s)=\frac{\mathbb{E}\left[\boldsymbol{\mathcal{X}}_{h-1}(s)\bm{\eta}_{h}\right]}{\mathbb{E}\left[\bm{\eta}_{h}^{2}\right]}$,
and
$\zeta_{h}(t)=\frac{\mathbb{E}\left[\mathcal{Y}_{h-1}(t)\bm{\eta}_{h}\right]}{\mathbb{E}\left[\bm{\eta}_{h}^{2}\right]}$.
Then, the $h$th FPLS component, denoted by $\bm{\eta}_{h}$, is obtained by
maximizing the Tucker’s criterion (2):

$\bm{\eta}_{h}=\int_{S}\boldsymbol{\mathcal{X}}_{h-1}^{\top}(s)\bm{\kappa}_{h}(s)\textrm{d}s,$

where the weight function $\bm{\eta}_{h}$ is given by:

$\bm{\kappa}_{h}(s)=\frac{\int_{T}\mathbb{E}\left[\boldsymbol{\mathcal{X}}_{h-1}(s)\mathcal{Y}_{h-1}(t)\right]\textrm{d}t}{\sqrt{\int_{S}\left(\int_{T}\mathbb{E}\left[\boldsymbol{\mathcal{X}}_{h-1}(s)\mathcal{Y}_{h-1}(t)\right]\textrm{d}t\right)^{2}\textrm{d}s}}.$

At the final step, the FPLS regression is completed by conducting the linear
regressions of $\boldsymbol{\mathcal{X}}_{h-1}(s)$ and $\mathcal{Y}_{h-1}(t)$
on $\bm{\eta}_{h}$.

Practically, the FPLS components of the FFLRM in (1) are obtained based on the
randomly observed sample curves of the functional response and functional
predictor variables. Although these sample curves are intrinsically infinite-
dimensional, they are observed in a set of discrete-time points in practice.
For this reason, obtaining the sample estimates of the Escoufier’s operators
becomes problematic because they need to be estimated from the discretely
observed observations (see Aguilera et al., 2010). To address this issue,
Preda & Schiltz (2011) and Beyaztas & Shang (2020a) proposed to use standard
PLS regression of the basis expansion coefficients of the functional variables
to approximate the FPLS components. In this context, several basis expansion
methods, such as $B$-spline, Fourier, wavelet, radial, and Bernstein
polynomial bases, have been proposed to project infinite-dimensional
functional variables onto the finite-dimensional space. In this study, we
consider the $B$-spline basis expansion coefficients of the functional
variables, which are then modeled using the IRSIMPLS regression of Alin &
Agostinelli (2017) to approximate the FPLS regression of $\mathcal{Y}(t)$ on
$\boldsymbol{\mathcal{X}}(s)$.

The IRSIMPLS is a weighted counterpart of the usual SIMPLS (see, e.g., Alin &
Agostinelli, 2017, for more information about the SIMPLS and IRSIMPLS
algorithms), where the weights are calculated based on the weighted likelihood
methodology of Markatou (1996). The weighted likelihood uses weighted score
equations to measure the discrepancy between the estimated and hypothesized
model densities.

Let us consider the FFLRM in (1) and let $K_{\mathcal{Y}}$ and
$K_{\boldsymbol{\mathcal{X}}}$ denote the numbers of basis functions used for
approximating the response and predictor functions, respectively. Then, they
can be represented as basis function expansion as follows:

$\displaystyle\mathcal{Y}(t)$
$\displaystyle=\sum_{k=1}^{K_{\mathcal{Y}}}c_{k}\phi_{k}(t)=\boldsymbol{C}^{\top}\boldsymbol{\Phi}(t),$
$\displaystyle\boldsymbol{\mathcal{X}}(s)$
$\displaystyle=\sum_{j=1}^{K_{\boldsymbol{\mathcal{X}}}}\bm{d}_{j}\bm{\psi}_{j}(s)=\boldsymbol{D}^{\top}\boldsymbol{\Psi}(s),$

where
$\boldsymbol{\Phi}(t)=\left[\phi_{1}(t),\ldots,\phi_{K_{\mathcal{Y}}}(t)\right]^{\top}$
and
$\boldsymbol{\Psi}(s)=\left[\bm{\psi}_{1}(s),\ldots,\bm{\psi}_{M}(s)\right]^{\top}$
with
$\bm{\psi}_{m}(s)=\left[\psi_{m1}(s),\ldots,\psi_{mK_{\boldsymbol{\mathcal{X}}}}(s)\right]^{\top}$
respectively denote the basis functions, and
$\boldsymbol{C}=\left[c_{1},\ldots,c_{K_{\mathcal{Y}}}\right]^{\top}$ and
$\boldsymbol{D}=\left[\bm{d}_{1},\ldots,\bm{d}_{M}\right]^{\top}$ with
$\bm{d}_{m}=\left[d_{m1},\ldots,d_{mK_{\boldsymbol{\mathcal{X}}}}\right]^{\top}$
are the corresponding coefficient matrices.

Let
$\boldsymbol{\Phi}=\int_{T}\boldsymbol{\Phi}(t)\boldsymbol{\Phi}^{\top}(t)dt$
and
$\boldsymbol{\Psi}=\int_{S}\boldsymbol{\Psi}(s)\boldsymbol{\Psi}^{\top}(s)ds$,
respectively, denote the $K_{\mathcal{Y}}\times K_{\mathcal{Y}}$ and
$K_{\boldsymbol{\mathcal{X}}}\times K_{\boldsymbol{\mathcal{X}}}$ dimensional
symmetric block-diagonal matrices of the inner products of the basis
functions. Denote by $\boldsymbol{\Phi}^{1/2}$ and $\boldsymbol{\Psi}^{1/2}$
the square roots of $\boldsymbol{\Phi}$ and $\boldsymbol{\Psi}$, respectively.
Then, we consider the following regression model to approximate the FPLS
regression of $\mathcal{Y}(t)$ on $\boldsymbol{\mathcal{X}}(s)$:

$\boldsymbol{\Lambda}=\boldsymbol{\Pi}\boldsymbol{\Omega}+\boldsymbol{\epsilon},$
(3)

where $\boldsymbol{\Lambda}=\boldsymbol{C}^{\top}\boldsymbol{\Phi}^{1/2}$,
$\boldsymbol{\Pi}=\boldsymbol{D}^{\top}\boldsymbol{\Psi}^{1/2}$, and
$\boldsymbol{\Omega}$ and $\boldsymbol{\epsilon}$ denote the coefficient and
residual matrices, respectively. In (3), the between columns of the obtained
design matrices ($\boldsymbol{\Lambda}$ and $\boldsymbol{\Pi}$) become highly
correlated (multicollinearity) due to the nature of functional data. In
addition, even though the functional predictors are independent, because of
the high-dimensional structures of the design matrices, there may be a strong
correlation between $\boldsymbol{\Lambda}$ and $\boldsymbol{\Pi}$. The
traditional estimation techniques such as LS and maximum likelihood extended
to functional data suffer from the multicollinearity problem. They may not
produce a valid estimate for the coefficient matrix $\boldsymbol{\Omega}$
(see, e.g., Beyaztas & Shang, 2020a). On the other hand, the PLS approach can
produce an efficient estimate for $\boldsymbol{\Omega}$ by eliminating the
multicollinearity problems. It has been showed by Preda & Schiltz (2011) and
Beyaztas & Shang (2021) that the FPLS regression of $\mathcal{Y}(t)$ on
$\boldsymbol{\mathcal{X}}(s)$ is equivalent to the PLS regression of
$\boldsymbol{\Lambda}$ on $\boldsymbol{\Pi}$ in the sense that at each step
$h$ of the PLS algorithm $1\leq h\leq H$, the same PLS components are obtained
for both PLS regressions. $H$ denotes the retained number of PLS components.

We consider the IRSIMPLS of Alin & Agostinelli (2017) to approximate the
coefficient function $\boldsymbol{\beta}(s,t)$ in (1). Similar to the
traditional PLS algorithms, the IRSIMPLS assumes a linear correlation between
predictors and the response and bypasses the multicollinearity problem. In the
IRSIMPLS, we assume that the theoretical residuals
$(\boldsymbol{\epsilon}_{1},\dots,\boldsymbol{\epsilon}_{n})$ belong to a
parametric family of densities
$\mathcal{F}=\left\\{f(\cdot,\sigma);\sigma\in\mathbb{R}^{+}\right\\}$ where
$f(\cdot,\sigma)$ is a random variable with mean zero corresponding to given
$\sigma$. Let $\widehat{\boldsymbol{\epsilon}}_{i}$ denote the calculated
residuals for the given estimates $\widehat{\boldsymbol{\Omega}}$ and
$\widehat{\sigma}$, appropriately determined by the density
$f(\cdot,\widehat{\sigma})$. Denote by
$g^{*}(\widehat{\boldsymbol{\epsilon}}_{i})$ and
$f^{*}(\widehat{\boldsymbol{\epsilon}}_{i};\sigma)$ the kernel density
estimator and the smoothed model density, respectively, as follows:

$\displaystyle g^{*}(\widehat{\boldsymbol{\epsilon}}_{i})$
$\displaystyle=\int\tau\left(\widehat{\boldsymbol{\epsilon}}_{i};u,v\right)\textrm{d}\widehat{F}(u),$
$\displaystyle f^{*}(\widehat{\boldsymbol{\epsilon}}_{i};\sigma)$
$\displaystyle=\int\tau\left(\widehat{\boldsymbol{\epsilon}}_{i};u,v\right)\textrm{d}\mathcal{F}(u;\sigma),$

where $\widehat{F}(u)$ denotes the empirical distribution function of the
calculated residuals $\widehat{\boldsymbol{\epsilon}}_{i}$ for $i=1,\ldots,n$
and $\tau\left(\widehat{\boldsymbol{\epsilon}}_{i};u,v\right)$ is the Gaussian
kernel density as follows:

$\tau\left(\widehat{\boldsymbol{\epsilon}}_{i};u,v\right)=\frac{e^{-(\widehat{\boldsymbol{\epsilon}}_{i}-u)^{2}/(2v^{2})}}{\sqrt{2\pi}v}.$

Herein, $v$ denotes the smoothing parameter, and in this study, it is chosen
as $v=\sqrt{\gamma\sigma^{2}}$, where $\gamma$ is a constant, to assign small
weights to outliers (Markatou, 1996; Alin & Agostinelli, 2017). Then, the
Pearson residuals $\delta\left(\widehat{\boldsymbol{\epsilon}}_{i}\right)$ for
$i=1,\ldots,n$, which are used to calculate weights, are defined as follows:

$\delta\left(\widehat{\boldsymbol{\epsilon}}_{i}\right)=\frac{g^{*}(\widehat{\boldsymbol{\epsilon}}_{i})}{f^{*}(\widehat{\boldsymbol{\epsilon}}_{i};\sigma)}.$

Finally, the weight for $i$th observation, $\omega_{i}$, is defined as
follows:

$\displaystyle\omega_{i}$
$\displaystyle=\omega\left(\widehat{\boldsymbol{\epsilon}}_{i};f^{*}(\widehat{\boldsymbol{\epsilon}}_{i};\sigma),\widehat{F}\right)$
$\displaystyle=\min\left\\{1,\frac{\left[A\left(\delta\left(\widehat{\boldsymbol{\epsilon}}_{i}\right)\right)+1\right]^{+}}{\delta\left(\widehat{\boldsymbol{\epsilon}}_{i}\right)+1}\right\\},$

where $A(\cdot)$ denotes the residual adjustment function. In this study, the
weights are calculated based on the Hellinger residual adjustment function of
Lindsay (1994): $A[\delta(\cdot)]=2[\delta(\cdot)+1]^{1/2}-1$.

The IRSIMPLS produces PLS components by deflating a weighted covariance matrix
$\boldsymbol{S}^{\omega}=\left(\boldsymbol{\Pi}^{\omega}\right)^{\top}\boldsymbol{\Lambda}^{\omega}$.
The $\boldsymbol{\Pi}^{\omega}$ and $\boldsymbol{\Lambda}^{\omega}$ are the
weighted versions of $\boldsymbol{\Pi}$ and $\boldsymbol{\Lambda}$,
respectively, and obtained by multiplying each row of $\boldsymbol{\Pi}$ and
$\boldsymbol{\Lambda}$ by the square root of weights:

$\omega_{i}=\underset{k}{\text{median}}\\{\omega_{ik}\\},\qquad
i=1,\ldots,n,\quad k=1,\ldots,K_{\mathcal{Y}}.$

Herein, $\omega_{ik}$ denotes the the weight for $i$th residual corresponding
to the $k$th column. The algorithm of the IRSIMPLS is described in Algorithm
1.

  * Step 1.

Center the matrices $\boldsymbol{\Pi}$ and $\boldsymbol{\Lambda}$ using
$L_{1}$ median as defined by Serrneels et al. (2005). In addition, scale the
centered $\boldsymbol{\Pi}$ using median absolute deviation. Let
$\boldsymbol{\Pi}^{*}$ denote the centered and scaled version of
$\boldsymbol{\Pi}$.

  * Step 2.

For $i=1,\ldots,n$, compute weights $\omega_{i}$ for each row of
$\boldsymbol{\Pi}^{*}$, and obtain $\boldsymbol{\Pi}^{\omega}$ and
$\boldsymbol{\Lambda}^{\omega}$ by multiplying each row of $\boldsymbol{\Pi}$
and $\boldsymbol{\Lambda}$ by the square roots of $\omega_{i}$.

  * Step 3.

Perform ordinary SIMPLS on a sample of $n_{r}$ ($n_{r}\ll n$) randomly drawn
pairs without replacement from
$\left(\boldsymbol{\Pi},\boldsymbol{\Lambda}\right)$ to obtain starting
coefficient estimate $\widehat{\boldsymbol{\Omega}}^{*}$. Note that the
starting values to calculate $\widehat{\boldsymbol{\Omega}}^{*}$ are
determined by applying ordinary SIMPLS on the centered and scaled
$\boldsymbol{\Pi}$ and $\boldsymbol{\Lambda}$.

  * Step 4.

Compute residuals using the starting coefficient estimate obtained in the
previous step as follows:

$\widehat{\boldsymbol{\epsilon}}_{i}=\boldsymbol{\Lambda}-\boldsymbol{\Pi}\widehat{\boldsymbol{\Omega}}^{*}.$

  * Step 5.

Obtain the centered and scaled residuals
$\widehat{\boldsymbol{\epsilon}}^{*}_{i}$ using $L_{1}$ median and median
absolute deviation.

  * Step 6.

Calculate new weights, $\omega^{*}_{i}$, using the residuals
$\widehat{\boldsymbol{\epsilon}}^{*}_{i}$ obtained in Step 5 to reweight
$\boldsymbol{\Pi}$ and $\boldsymbol{\Lambda}$ multiplying each row of them by
the square roots of $\omega^{*}_{i}$. Let $\boldsymbol{\Pi}^{\omega*}$ and
$\boldsymbol{\Lambda}^{\omega*}$ denote the reweighted versions of
$\boldsymbol{\Pi}$ and $\boldsymbol{\Lambda}$, respectively.

  * Step 7.

Obtain the reweighted covariance matrix
$\boldsymbol{S}^{\omega*}=\left(\boldsymbol{\Pi}^{\omega*}\right)^{\top}\boldsymbol{\Lambda}^{\omega*}$.
Then, compute the new coefficient estimate
$\widehat{\boldsymbol{\Omega}}^{**}$ performing the ordinary SIMPLS on
$\boldsymbol{S}^{\omega*}$.

  * Step 8.

With the coefficient estimate in Step 7, repeat steps 4-7 until convergence is
achieved. Convergence is reached if the maximum of the absolute deviations
between two $\widehat{\boldsymbol{\Omega}}^{**}$s obtained consecutive repeats
less than a predetermined threshold, such as $10^{-4}$.

  * Step 9.

Repeat Steps 3-8 a few times to obtain the final estimate of
$\boldsymbol{\Omega}$, given by $\boldsymbol{\Omega}^{\omega}$, which is equal
to one having the minimum absolute deviation among others.

Algorithm 1 IRSIMPLS algorithm

Let $\boldsymbol{\Omega}^{\omega,h}$ denote the estimated coefficient matrix
of $\boldsymbol{\Omega}$, obtained after $h$ iterations from the IRSIMPLS
regression. Then, the IRSIMPLS approximation of $\boldsymbol{\beta}(s,t)$,
given by $\widehat{\boldsymbol{\beta}}^{\omega}(s,t)$, is obtained using the
bases $\boldsymbol{\Phi}(t)$ and $\boldsymbol{\Psi}(s)$ as follows:

$\widehat{\boldsymbol{\beta}}^{\omega}(s,t)=\left[\left(\boldsymbol{\Psi}^{1/2}\right)^{-1}\boldsymbol{\Omega}^{\omega,h}\left(\boldsymbol{\Phi}^{1/2}\right)^{-1}\right]\boldsymbol{\Psi}(s)\boldsymbol{\Phi}(t).$
(4)

In what follows, the robust estimation of the functional response,
$\widehat{\mathcal{Y}}(t)$, can be calculated as follows:

$\widehat{\mathcal{Y}}(t)=\int_{S}\bm{\mathcal{X}}(s)\widehat{\boldsymbol{\beta}}^{\omega}(s,t)\textrm{d}s.$

Since the IRSIMPLS approximation of $\boldsymbol{\beta}(s,t)$ given in (4) is
based on the weighted residuals, it is more robust to the outliers than the
ordinary LS-based PLS methods.

The performance of the proposed method depends on the number of PLS components
$h$ used to estimate Model (1). Thus, the optimum number $h$ should be
determined based on an error metric. The outliers may affect the usual error
metrics, leading to incorrectly determining the optimum number $h$. Therefore,
we consider a trimmed mean absolute prediction percentage error (TMAPE) to
determine the optimum $h$ robustly. To this end, we consider the following
procedure. First, the entire dataset is randomly divided into two parts
roughly the same size. The first part of the data is used to construct the
FFLRM, and the PLS approximation of $\boldsymbol{\beta}(s,t)$ is obtained for
$h=1,2,\ldots,H$ number of PLS components. With the second part of the data,
the absolute prediction percentage errors (APE) are calculated as follows:

$\text{APE}_{i}(h)=\left\|\frac{\left|\widehat{\mathcal{Y}}_{i}(t)-\mathcal{Y}_{i}(t)\right|}{\left|\mathcal{Y}_{i}(t)\right|}\right\|_{\mathcal{L}_{2}}\qquad
i=1,2,\ldots,\lfloor n/2\rfloor,$

where $\left\|\cdot\right\|_{\mathcal{L}_{2}}$ denotes $\mathcal{L}_{2}$ norm,
which is approximated by the Riemann sum (Luo & Qi, 2019),
$\lfloor\cdot\rfloor$ denotes an integer value that rounds down to $n/2$, and
$\mathcal{Y}_{i}(t)$ and $\widehat{\mathcal{Y}}_{i}(t)$ are the observed
response functions for $i^{\text{th}}$ individual and its prediction,
respectively. Then, the following TMAPE is computed for each $h$:

$\text{TMAPE}(h)=\frac{1}{n*}\sum_{i=1}^{n*}\text{APE}_{i}(h),$

where $\left\\{n*\subset\\{1,2,\ldots,\lfloor
n/2\rfloor\\},~{}|n*|=\left[q\lfloor n/2\rfloor\right]\right\\}$. Herein, we
consider 20% trimming proportion as suggested by Wilcox (2012) so that $q=0.8$
(other trimming proportions such as 5% and 10% can also be used in the
proposed method). In this way, the optimum $h$ is robustly determined by
ignoring the APE values corresponding to the outliers. Finally, the optimum
number of PLS components is determined according to the minimum TMAPE.

##  3 Monte Carlo experiments

Computationally, the proposed method works with a similar structure with the
FPLS methods of Preda & Schiltz (2011) and Beyaztas & Shang (2020a). We note
that all the numerical calculations in this study are performed using R 4.1.1.
on an Intel Core i7 6700HQ 2.6 GHz PC. We summarize the steps and R packages
used in the computation of the proposed method as follows. In the first step
of the proposed method, the $B$-spline basis expansion method is used to
project the functional random variables onto the finite-dimensional space,
which can easily be done using create.bspline.basis function in the R package
fda (Ramsay et al., 2020). In addition, the inprod and sqrtm functions in the
R packages fda and expm (Goulet et al., 2021) are used to obtain the inner
product matrices and their square roots, respectively. In the second step, the
IRSIMPLS algorithm proposed by (Alin & Agostinelli, 2017) uses the R package
wle (Agostinelli & Library, 2015) to obtain the weights. These weights are
used to robustly estimate the regression parameter matrix of the PLS
regression constructed using basis expansion coefficients. The estimated
parameter matrix and $B$-spline basis functions are then used to obtain the
approximate form of the regression parameter function $\bm{\beta}(s,t)$.

Various Monte Carlo experiments under different data generating processes are
conducted to evaluate the finite-sample performance of the proposed method.
The results are compared with those obtained via the LS, naïve SIMPLS, and the
RFPC of Hullait et al. (2021). Throughout the experiments, the following FFLRM
is considered:

$\mathcal{Y}(t)=\sum_{m=1}^{M=5}\int_{S}\mathcal{X}_{m}(s)\beta_{m}(s,t)\textrm{d}s+\epsilon(t),$
(5)

where $s,t\in[0,1]$. Two cases are considered for generating the data:

Case-1:

The functions of both the response and predictor variables are generated from
smooth processes with no outliers. In this case, the proposed method is
expected to produce similar results with available traditional methods. Our
aim, in this case, is to show the correctness of the proposed method.

Case-2:

The functional variables of the response and predictor are contaminated by
inserting outliers at 20% contamination level. In this case, the proposed
method is expected to produce better prediction performance than the
traditional methods.

We consider two scenarios to generate the functions of the predictor
variables. In the first scenario (Scenario-1), the functional predictors are
assumed to be independent (i.e.,
$\text{Cov}(\mathcal{X}_{m}(s),\mathcal{X}_{m^{\prime}}(s))=0$ for $m\neq
m^{\prime}$ and $m,m^{\prime}=1,\ldots,5$) and are generated as follows:

$\mathcal{X}_{m}(s)=10+V_{m}(s),$

where $V_{m}(s)$ are generated from the Gaussian process with mean zero and
variance-covariance function
$\boldsymbol{\Sigma}_{V}(s,s^{\prime})=e^{-100(s-s^{\prime})^{2}}$. In the
second scenario (Scenario-2), on the other hand, the functional predictors are
generated as follows:

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

where $V_{m}(s)$ ($m=1,\ldots,9$) are generated as in the first scenario and
Lag controls the correlation between the functional predictors (Luo & Qi,
2019). In the experiments, similar to Luo & Qi (2019), we consider Lag = 4 to
generate correlated functional predictors. The functions of the response
variable are obtained using (5) with the following coefficient functions:

$\displaystyle\beta_{1}(s,t)$ $\displaystyle=(1-s)^{2}(t-0.5)^{2},$
$\displaystyle\beta_{2}(s,t)$
$\displaystyle=e^{-3(s-1)^{2}}e^{-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}.$

For each variable, $n=500$ functions are generated at 100 equally spaced
points in the interval $[0,1]$. Then, to evaluate the out-of-sample predictive
performance of the methods, the generated data are divided into the following
training and test samples with sample sizes $n_{\text{train}}$ and
$n_{\text{test}}$, respectively:  1) the first $n_{\text{train}}=200$
functions are used to construct the FFLRM, and  2) the last
$n_{\text{test}}=300$ functions are used to predict the response based on the
constructed model and a new set of predictors.  To generate the outlier-
contaminated data, $\left(0.2\times n_{\text{train}}\right)$ randomly selected
functions of the predictors are replaced by
$\widetilde{\mathcal{X}}_{m}(s)=\mathcal{X}_{m}(s)+\varepsilon_{m}(s)$, where
$\varepsilon_{m}(s)$ is generated from the Ornstein-Uhlenbeck process:

$\varepsilon_{m}(s)=\gamma+[\varepsilon_{m0}(s)-\rho]e^{-\theta
s}\sigma\int_{0}^{s}e^{-\theta(s-u)}\textrm{d}W_{u},$

where $\rho$, $\theta$, and $\sigma>0$ are constants, and $W_{u}$ is the
Wiener process. Note that the initial values $\varepsilon_{m0}(s)$ are
independently taken from $W_{u}$. Then, the outlying functions of the response
variable are obtained using (5) but with $\widetilde{\mathcal{X}}_{m}(s)$ and
the following modified regression coefficient functions (only
$\beta_{2}(s,t)$, $\beta_{3}(s,t)$, and $\beta_{4}(s,t)$ are modified):

$\displaystyle\widetilde{\beta}_{2}(s,t)$
$\displaystyle=e^{-3(s-1)^{2}}e^{-3(t-0.5)^{2}},$
$\displaystyle\widetilde{\beta}_{3}(s,t)$
$\displaystyle=6e^{-5(s+0.5)^{2}-5(t-0.5)^{2}}+4e^{-5(s-1.5)^{2}-5(t-0.5)^{2}},$
$\displaystyle\widetilde{\beta}_{4}(s,t)$ $\displaystyle=6\cos(6\pi s)\cos(\pi
t).$

In this way, the functional predictors are contaminated by the magnitude
outliers; the response variable is contaminated by the mixture of the
magnitude and shape outliers. Examples of the generated functions and the
outliers are presented in Figure 1.

Figure 1: Plots of the generated functions of the response variable (left
panel) and predictor variables (right panel). Curves in black color present
the non-outlying functions, while gray curves present the generated outlying
functions

The functions of both the response and predictor variables in the training
sample are distorted by the Gaussian random noise
$\upsilon(s)\sim\mathcal{N}(0,4)$ before fitting the model. Also,
$K_{\mathcal{Y}}=K_{\boldsymbol{\mathcal{X}}}=10$ number of basis functions
are used to smooth the noisy data. Note that in the Monte Carlo experiments,
the number of basis functions
$K_{\mathcal{Y}}=K_{\boldsymbol{\mathcal{X}}}=10$ is chosen arbitrarily to
evaluate the predictive performance of the methods under the same conditions.
For each scenario, MC = 250 Monte Carlo experiments are performed. In each
experiment, three performance metrics; mean absolute prediction percentage
error (MAPE), median absolute prediction error (MdAPE), and coefficient of
determination $R^{2}$ are computed to evaluate the out-of-sample prediction
performance of the methods:

MAPE
$\displaystyle=\frac{1}{n_{\text{test}}}\sum_{i=1}^{n_{\text{test}}}\left\|\frac{\left|\mathcal{Y}_{i}(t)-\widehat{\mathcal{Y}}_{i}(t)\right|}{\left|\mathcal{Y}_{i}(t)\right|}\right\|_{\mathcal{L}_{2}},$
MdAPE
$\displaystyle=\text{median}\left\\{\text{AE}_{1}^{p},\ldots,\text{AE}_{n_{\text{test}}}^{p}\right\\},$
$\displaystyle R^{2}$
$\displaystyle=1-\frac{\sum_{i=1}^{n_{\text{test}}}\left\|\mathcal{Y}_{i}(t)-\widehat{\mathcal{Y}}_{i}(t)\right\|_{\mathcal{L}_{2}}^{2}}{\sum_{i=1}^{n_{\text{test}}}\left\|\mathcal{Y}_{i}(t)-\bar{\mathcal{Y}}(t)\right\|_{\mathcal{L}_{2}}^{2}},$

where $\mathcal{Y}_{i}(t)$ and $\widehat{\mathcal{Y}}_{i}(t)$ respectively
denote the $i^{\textsuperscript{th}}$ generated and predicted functions of the
response variable in the test sample and
$\text{AE}_{i}=|\mathcal{Y}_{i}(t)-\widehat{\mathcal{Y}}_{i}(t)|$ denotes the
absolute prediction error.

For further investigation of the uncertainty of predictions, the nonparametric
bootstrap method is applied to construct pointwise prediction intervals for
the functions of the response variable in the test sample. Let
$\widehat{\epsilon}(t)=\mathcal{Y}(t)-\widehat{\mathcal{Y}}(t)$ denote the
estimated model error. The following algorithm then describes the bootstrap
method used in this study to construct pointwise prediction intervals for the
functions of the response variable.

  * Step 1.

Draw a pair of bootstrap samples
$\left(\mathcal{Y}^{*}(t),\boldsymbol{\mathcal{X}}^{*}(s)\right)$ of size
$n_{\text{train}}$ by sampling with replacement from the original sample
$\left(\mathcal{Y}(t),\boldsymbol{\mathcal{X}}(s)\right)$.

  * Step 2.

Using the bootstrap sample
$\left(\mathcal{Y}^{*}(t),\boldsymbol{\mathcal{X}}^{*}(s)\right)$, obtain the
IRSIMPLS approximation of $\boldsymbol{\beta}(s,t)$, namely
$\widehat{\boldsymbol{\beta}}^{\omega*}(s,t)$ as described in Section 2.

  * Step 3.

Compute the bootstrap predictions of $\mathcal{Y}(t)$ in the test sample as
follows:

$\widehat{\mathcal{Y}}^{*}(t)=\int_{S}\boldsymbol{\mathcal{X}}(s)\widehat{\boldsymbol{\beta}}^{\omega*}(s,t)\textrm{d}s+\epsilon^{*}(t),$

where $\boldsymbol{\mathcal{X}}(s)$ denotes the matrix of predictors in the
test sample, and $\epsilon^{*}(t)$ is a random drawn from
$\widehat{\epsilon}(t)$.

  * Step 4.

Repeat Steps 1-3 $B$ times, where $B=200$ denotes the number of bootstrap
simulations, to obtain $B$ sets of bootstrap replicates of the predicted
response functions
$\left\\{\widehat{\mathcal{Y}}^{*,1}(t),\ldots,\widehat{\mathcal{Y}}^{*,B}(t)\right\\}$.

Denote by $Q_{\alpha}^{i}(t)$ the $\alpha$th quantile of the generated $B$
sets of bootstrap replicates of the $i$th predicted response function, where
$\alpha$ denotes a level of significance. In our numerical analyses, the
significance level is set to $\alpha=0.05$ to construct 95% prediction
intervals for the functions of the response variable in the test sample. Then,
$\left[Q_{\alpha/2}^{i}(t),Q_{1-\alpha/2}^{i}(t)\right]$ defines the
$100(1-\alpha)\%$ bootstrap prediction interval for $\mathcal{Y}_{i}(t)$. For
evaluating the performance of the methods, two bootstrap performance metrics;
the coverage probability deviance (CPD), which is the absolute difference
between the nominal and empirical coverage probabilities, and the interval
score (score), are calculated as follows:

CPD
$\displaystyle=\left|\alpha-\frac{1}{n_{\text{test}}}\sum_{i=1}^{n_{\text{test}}}\mathbb{1}\left\\{Q_{\alpha/2}^{i}(t)>\mathcal{Y}_{i}(t)\right\\}+\mathbb{1}\left\\{Q_{1-\alpha/2}^{i}(t)<\mathcal{Y}_{i}(t)\right\\}\right|,$
score
$\displaystyle=\frac{1}{n_{\text{test}}}\sum_{i=1}^{n_{\text{test}}}\bigg{\|}\left[\left\\{Q_{1-\alpha/2}^{i}(t)-Q_{\alpha/2}^{i}(t)\right\\}\right.$
$\displaystyle+\frac{2}{\alpha}\left(Q_{\alpha/2}^{i}(t)-\mathcal{Y}_{i}(t)\right)\mathbb{1}\left\\{\mathcal{Y}_{i}(t)<Q_{\alpha/2}^{i}(t)\right\\}$
$\displaystyle+\left.\frac{2}{\alpha}\left(\mathcal{Y}_{i}(t)-Q_{1-\alpha/2}^{i}(t)\right)\mathbb{1}\left\\{\mathcal{Y}_{i}(t)>Q_{1-\alpha/2}^{i}(t)\right\\}\right]\bigg{\|}_{\mathcal{L}_{2}},$

where $\mathbb{1}\\{\cdot\\}$ denotes the binary indicator function.

For the aforementioned data generation process, we compare the out-of-sample
predictive performance of our proposed method only with the LS and naïve
SIMPLS-based functional PLS. The error measures for the RFPC are not
calculated since it allows only one functional predictor in the model. The
finite-sample performance of the proposed method will be compared with the
RFPC using a different data generation process.

Our findings obtained from the Monte Carlo experiments are presented in
Figures 2 and 3. From these figures, it is obvious that all the methods
produce slightly better results when the functional predictors are independent
(Scenario-1) than when the functional predictors are correlated (Scenario-2).
Figure 2 demonstrates that when no outliers are present in the data, all the
methods produce similar MAPE, MdAPE, and $R^{2}$ values. The traditional
SIMPLS method produces slightly better performance metrics than the proposed
method, but the observed differences are not significant. In addition, from
Figure 3, when outliers do not contaminate the data, the proposed method
produces similar bootstrap-based CPD values with slightly larger score values
compared with LS and SIMPLS.

Figure 2: Boxplots of the calculated MAPE (first row), MdAPE (second row), and
$R^{2}$ (third row) values for the LS, naïve SIMPLS, and proposed IRSIMPLS
methods under Scenario-1 (first column) and Scenario-2 (second column)

Figure 3: Boxplots of the calculated CPD (first row) and score (second row)
values for the LS, naïve SIMPLS, and proposed IRSIMPLS methods under
Scenario-1 (first column) and Scenario-2 (second column)

When the generated data include outliers, our proposed method produces
considerably smaller MAPE, MdAPE, and larger $R^{2}$ values under both cases
and scenarios compared with the LS and naïve SIMPLS. From Figure 2, the LS and
naïve SIMPLS are affected by the outliers and produce larger error values
compared with the case where outliers do not contaminate the data. On the
other hand, the proposed method down-weights the effects of outliers and
produces almost similar error values under both cases. Moreover, compared with
LS and traditional SIMPLS, the proposed method produces significantly smaller
CPD values with slightly larger score values when outliers contaminate the
data (see Figure 3). In other words, compared with LS and traditional SIMPLS,
much accurate prediction intervals are obtained by the proposed method when
the data have outliers.

Moreover, we compare the proposed method with LS and naïve SIMPLS in terms of
their computing times. While doing so, we consider Case-1 and Scenario-1 as
data generation processes and record the computing times of the methods using
the R function system.time when the numbers of basis function
$K_{\mathcal{Y}}$ and $K_{\bm{\mathcal{X}}}$ are increased from 10 to 80. The
results are presented in Figure 4. From this figure, the proposed method
requires more computing time compared with other methods. For example, the
proposed method respectively requires 1.65 and 1.18 times more computing time
than the LS and naïve SIMPLS when $K_{\mathcal{Y}}=K_{\bm{\mathcal{X}}}=10$
while it requires 2.45 and 1.71 times more computing time when
$K_{\mathcal{Y}}=K_{\bm{\mathcal{X}}}=80$.

Figure 4: Estimated computing times for the LS, SIMPLS, and IRSIMPLS methods

All in all, the numerical results produced by our Monte Carlo experiments have
demonstrated that:  1) the proposed IRSIMPLS-based FFLRM has a similar
prediction performance to that of LS and traditional SIMPLS-based regression
when no outliers are present in the data, and  2) the proposed method has
improved prediction performance than the LS and traditional SIMPLS when the
data are contaminated by magnitude or shape outliers.  This is because the
traditional approximation of $\boldsymbol{\beta}(s,t)$ is significantly
affected by outliers, leading to poor prediction performance. However,
compared with LS and traditional SIMPLS, the proposed method increases the
prediction performance of FFLRM by down-weighting the effects of outliers.

To compare the performance of our proposed method with the RFPC, we use a
similar data generation process, which is a simple FFLRM, as used by Hullait
et al. (2021). The predictor variable $\mathcal{X}(s)$ is generated based on a
FPC analysis with the following mean function ($\mu_{\mathcal{X}}(s)$) and
eigenfunctions ($\iota_{i}^{\mathcal{X}}$ for $i=1,2,3$):

$\displaystyle\mu_{\mathcal{X}}(s)$ $\displaystyle=-10(s-0.5)^{2}+2,$
$\displaystyle\iota_{1}^{\mathcal{X}}$ $\displaystyle=\sqrt{2}\sin(\pi s),$
$\displaystyle\iota_{2}^{\mathcal{X}}$ $\displaystyle=\sqrt{2}\sin(7\pi s),$
$\displaystyle\iota_{3}^{\mathcal{X}}$ $\displaystyle=\sqrt{2}\cos(7\pi s),$

where $s\in[0,1]$ and the principal component scores are generated by sampling
from the Gaussian distributions with mean-zero and variances 40, 10, and 1,
respectively. The response variable $\mathcal{Y}(t)$, on the other hand, is
generated using the following mean function ($\mu_{\mathcal{Y}}(t)$) and
eigenfunctions ($\iota_{i}^{\mathcal{Y}}$ for $i=1,2,3$):

$\displaystyle\mu_{\mathcal{Y}}(t)$ $\displaystyle=60\exp(-(t-1)^{2}),$
$\displaystyle\iota_{1}^{\mathcal{Y}}$ $\displaystyle=\sqrt{2}\sin(12\pi t),$
$\displaystyle\iota_{2}^{\mathcal{Y}}$ $\displaystyle=\sqrt{2}\sin(5\pi t),$
$\displaystyle\iota_{3}^{\mathcal{Y}}$ $\displaystyle=\sqrt{2}\cos(2\pi t),$

where $t\in[0,1]$. The following simple FFLRM is used to generate functional
observations:

$\mathcal{Y}(t)=\mu_{\mathcal{Y}}(t)+\int_{0}^{1}\beta(s,t)[\mathcal{X}(s)-\mu_{\mathcal{X}}(s)]\textrm{d}s+\epsilon(t),$

where
$\beta(s,t)=\bm{\iota}^{\mathcal{X}}(s)\bm{B}\bm{\iota}^{\mathcal{Y}}(t)$ with
$\bm{B}$ having random entries between uniform $[-3,3]$,
$\epsilon(t)=\bm{q}^{\top}\bm{\iota}^{\mathcal{Y}}(t)+d$, and $\bm{q}^{\top}$
and $d$ are sampled from $\mathcal{N}(0,0.1)$. We generate $n=200$ functions
for each functional variable at 200 equally spaced points in the interval
$[0,1]$. The generated data are contaminated at 20% contamination level, and
two scenarios are considered to generate outliers:

Scenario-1

where the outliers are generated by replacing $\bm{B}$ with
$\bm{B}_{1}=\bm{B}+\bm{R}$ with $\bm{R}\sim\mathcal{N}(0,0.5)$.

Scenario-2

where the outliers are generated using
$\beta(s,t)=\bm{\iota}^{\mathcal{X}}(s)\bm{B}_{2}[\bm{\iota}^{\mathcal{Y}}(t),p(t)]$
where $\bm{B}_{2}$ is an $3\times 4$ matrix $\bm{B}_{2}=[\bm{B}_{2},\bm{I}]$,
$\bm{I}$ is sampled from $\mathcal{N}(2,1)$, and $p(t)$ is a random $B$-spline
function defined on an interval of length $1/10$.

As noted by Hullait et al. (2021), the outliers generated under Scenario-1
affect the regression function across the entire interval, whereas the
outliers generated under Scenario-2 only affect a small interval of the
curves. A graphical display of the generated data is presented in Figure 5.

Figure 5: Plots of the generated 50 functions for the functional response and
predictor variables

The simple FFLRM is constructed for both scenarios, and the bivariate
coefficient is estimated using the proposed method and RFPC based on $K=15$
$B$-spline basis functions. Similar to Hullait et al. (2021), the following
mean squared error (MSE) is considered to compare the performance of the
methods:

$\text{MSE}=\frac{1}{(1-0.20)\times
n}\sum_{i=1}^{n}\left\|(1-u_{i})\widehat{\mathcal{Y}}_{i}(t)-\mathcal{Y}_{i}(t)\right\|^{2}_{\mathcal{L}_{2}},$
(6)

where $\widehat{\mathcal{Y}}_{i}(t)$ is the $i\textsuperscript{th}$ fitted
response function and $u_{i}=1$ if $i\textsuperscript{th}$ response function
is outlier and $u_{i}=0$ otherwise.

The MSE values obtained from 100 replications are presented in Figure 6. From
this figure, our proposed method produces smaller MSE values compared with the
RFPC of Hullait et al. (2021).

Figure 6: Boxplots of the calculated MSE values for the RFPC and proposed
IRSIMPLS methods, when the data are generated under Scenario-1 (left panel)
and Scenario-2 (right panel)

##  4 Empirical data example: Oman weather data

The finite-sample performance of the proposed method is also evaluated using
the Oman weather empirical dataset (dataset is available from the National
Center for Statistics & Information: https://data.gov.om). The dataset
contains three variables: maximum monthly evaporation (mm), humidity (%), and
temperature (∘C), which were collected from 66 weather stations across Oman
from January 2018 to December 2018. The observations are considered as the
functions of months ($1\leq s,t\leq 12$). Figure 7 presents the graphical
display of the variables, as well as the functional highest density region
(HDR) boxplots (Hyndman & Shang, 2010) for these variables. In the functional
HDR boxplots, the inner and outer regions are presented with dark and light
gray colors for each variable, respectively. The inner region corresponds to
the region bounded by all the curves corresponding to the points inside the
50% bivariate HDR. In contrast, the outer region corresponds to the region
bounded by all the curves corresponding to the points within the outer
bivariate HDR (Hyndman & Shang, 2010). Functions that are outside of these
regions are called outliers. From this figure, all the variables have clear
outlying curves.

Figure 7: Graphical display of the temperature, humidity, evaporation
variables (first row), and their HDR boxplots (second row) for the Oman
weather data. The observations are the functions of months and $1\leq s,t\leq
12$

For this dataset, we consider predicting maximum monthly evaporation using the
maximum temperature and humidity variables. The following procedure is
repeated 100 times to evaluate the predictive performance of the methods. The
dataset is randomly divided into the training and test samples with sizes 40
and 26, respectively. The model is constructed based on the curves in the
training sample using $K_{\mathcal{Y}}=K_{\boldsymbol{\mathcal{X}}}=6$ number
of basis functions, which are determined based on the generalized cross-
validation procedure of Craven & Wahba (1978), to predict 26 curves in the
test sample. For each replication, the MAPE, MdAPE, and $R^{2}$, as well as
the bootstrap-based score and CPD values, are calculated.

Our results are presented in Figure 8, which indicates that the proposed
IRSIMPLS-based method produces better prediction performance compared with the
LS and naïve SIMPLS. The proposed method has smaller MAPE and MdAPE and larger
$R^{2}$ values than LS and classical SIMPLS. In addition, the proposed method
produces smaller bootstrap-based CPD values than the LS and classical SIMPLS.
However, it produces slightly larger score values than classical SIMPLS, as
presented in Figure 8.

Figure 8: Computed MAPE, $R^{2}$, bootstrap-based CPD, and score values for
the LS, traditional SIMPLS, and proposed IRSIMPLS methods for the Oman weather
data

To compare the performance of our method with the RFPC, we consider two simple
FFLRMs. In the first model, we consider temperature as the predictor and
evaporation as the response variable, whereas humidity is the predictor
variable in the second model. We construct models based on the entire data and
calculate the MSE values using the observed and fitted curves as in (6). The
computed MSE values for the first model are 21.75 and 17.87 for the RFPC and
IRSIMPLS, respectively. For the second model, the computed MSE values are
observed as 20.43 and 21.98. Our proposed method can produce competitive
performance compared with the RFPC for the Oman weather data.

##  5 Conclusions, limitations, and future research

FFLRMs have become an important analytical tool to explore the relationships
between complex and high-dimensional functional variables. Accordingly,
several FFLRMs have been proposed. The traditional estimation techniques,
including PLS regression, have been extended to $\mathcal{L}_{2}$ Hilbert
space to estimate the coefficient function of the considered regression
models. However, these techniques may be significantly affected by the
presence of outliers. In this case, traditional methods produce biased
estimates for the regression parameter, leading to poor prediction results.

In this study, we propose a robust PLS method based on the IRSIMPLS approach
of Alin & Agostinelli (2017) to obtain an outlier-resistant estimate for the
coefficient function of the FFLRM. Its predictive performance is evaluated
using various Monte Carlo experiments and an empirical data example. The
results are compared with the LS, naïve SIMPLS, and the RFPC method of Hullait
et al. (2021). In addition, we applied a nonparametric bootstrap to obtain
pointwise prediction intervals for the functions of the response variable. Our
findings demonstrate that the proposed method produces similar results to LS
and naïve SIMPLS when the data do not include outliers. On the other hand, our
proposed method produces better prediction performance than the LS and
traditional SIMPLS when the data contain outliers. The numerical results
obtained from the Monte Carlo experiments performed in this study have shown
that the proposed method produces more accurate bootstrap prediction intervals
with wider prediction interval lengths than the traditional methods in the
presence of outliers. Another interesting result produced by our Monte Carlo
experiments is that, similar to the traditional methods, the predictive
performance of the proposed method is affected by the correlation level
between the predictors. Compared to the case where the functional predictors
are correlated, our proposed method produces better predictive performance
when the predictors are independent. However, it still produces improved
performance over competitors when the predictors are correlated, and outliers
contaminate the data. Contrary to the superiority of the proposed method over
the traditional methods, it requires more computing time than the LS and naïve
SIMPLS methods since the proposed method uses an iterative weighted likelihood
algorithm. The results of our Monte Carlo experiments and empirical data
analysis have also shown that the proposed IRSIMPLS method produces improved
estimation and prediction accuracy than RFPC.

There are several directions that our proposed method can be further extended.
For example:  1) When a relatively large number of functional predictors are
used in the FFLRM, not all of them may significantly affect the dependent
variable. In such a case, the IRSIMPLS-based FPLS method may produce poor
predictive performance. Thus, similar to Luo & Qi (2019) and Beyaztas & Shang
(2021), a variable selection procedure along with the IRSIMPLS can be used to
improve the predictive performance of the proposed method.  2) Our numerical
analyses consider only the $B$-spline basis function expansion to project the
infinite-dimensional random functions into the finite-dimensional space. The
predictive performance of the proposed method may be dependent on the type of
basis functions under consideration. As a future study, the finite-sample
performance of our proposed method can be explored using other available basis
function expansion methods, such as Fourier, wavelet, radial, and Bernstein
polynomial bases.  3) In this study, we only consider the main effects of the
predictors on the response variable. However, recent studies (see, e.g., Luo &
Qi, 2019; Matsui, 2020; Beyaztas & Shang, 2021) have demonstrated that the
FFLRM with quadratic and interaction effects of the functional predictors
perform better than the standard FFLERM. On the other hand, the outliers in
the quadratic and interaction effects of the predictors may be more erroneous
than the effects of outliers in the main effects. The proposed method can also
be extended to the FFLRM with quadratic and interaction effects to estimate
this model robustly.  4) Our numerical analyses consider a standard bootstrap
method to construct a pointwise prediction interval for the response variable.
While our Monte Carlo experiments show that the proposed method produces
accurate prediction intervals in the presence of outliers, it produces wider
prediction intervals than traditional methods. A robust bootstrap method such
as the one proposed by Beyaztas & Shang (2020b) can be used along with the
proposed method to construct accurate prediction intervals with narrower
prediction interval lengths.  5) One of the recent studies in the functional
data analysis literature is about the misaligned observations (Olsen et al.,
2018). When misaligned functional observations are present in the data, they
may be considered as shape outliers, and such observations may lead to poor
modeling results. The proposed method can also be extended to the analysis of
misaligned functional observations to reduce the effects of such observations
on the analysis results.

## Acknowledgements

This work was supported by The Scientific and Technological Research Council
of Turkey (TUBITAK) (grant no: 120F270). We thank Dr. Aylin Alin for providing
some R code.

## References

  * (1)
  * Agostinelli & Library (2015) Agostinelli, C. & Library, S. C. M. (2015), wle: Weighted Likelihood Estimation.  R package version 0.9-91.
https://CRAN.R-project.org/package=wle

  * Aguilera et al. (2010) Aguilera, A. M., Escabias, M., Preda, C. & 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. & 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. & Valderrama, M. J. (1999), ‘Forecasting with unequally spaced data by a functional principal component approach’, Test 8(1), 233–254.
  * Alin & Agostinelli (2017) Alin, A. & Agostinelli, C. (2017), ‘Robust iteratively reweighted SIMPLS’, Journal of Chemometrics 31(3), e2881.
  * Beyaztas & Shang (2019) Beyaztas, U. & Shang, H. L. (2019), ‘Forecasting functional time series using weighted likelihood methodology’, Journal of Statistical Computation and Simulation 89(16), 3046–3060.
  * Beyaztas & Shang (2020a) Beyaztas, U. & Shang, H. L. (2020a), ‘On function-on-function regression: Partial least squares approach’, Environmental and Ecological Statistics 27(1), 95–114.
  * Beyaztas & Shang (2020b) Beyaztas, U. & Shang, H. L. (2020b), ‘Robust bootstrap prediction intervals for univariate and multivariate autoregressive time series models’, Journal of Applied Statistics in-press.
  * Beyaztas & Shang (2021) Beyaztas, U. & Shang, H. L. (2021), ‘A partial least squares approach for function-on-function interaction regression’, Cmputational Statistics 36, 911–939.
  * Chiou et al. (2016) Chiou, J.-M., Yang, Y.-F. & Chen, Y.-T. (2016), ‘Multivariate functional linear regression and prediction’, Journal of Multivariate Analysis 146, 301–312.
  * Craven & Wahba (1978) Craven, P. & Wahba, G. (1978), ‘Smoothing noisy data with spline functions’, Numerische Mathematik 31(4), 377–703.
  * 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.
  * Cummins & Andrews (1995) Cummins, D. J. & Andrews, C. W. (1995), ‘Iteratively reweighted partial least squares: A performance analysis by Monte Carlo simulation’, Journal of Chemometrics 9(6), 489–507.
  * Dayal & MacGregor (1997) Dayal, B. S. & 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.
  * Delaigle & Hall (2012) Delaigle, A. & Hall, P. (2012), ‘Methodology and theory for partial least squares applied to functional data’, The Annals of Statistics 40(1), 322–352.
  * 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. & Gonzalez‐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 & Vieu (2006) Ferraty, F. & Vieu, P. (2006), Nonparametric Functional Data Analysis, Springer, New York.
  * Gervini (2012) Gervini, D. (2012), Functional robust regression for longitudinal data, Working paper, University of Wisconsin–Milwaukee.
https://arxiv.org/abs/1211.7332

  * Goulet et al. (2021) Goulet, V., Dutang, C., Maechler, M., Firth, D., Shapira, M. & Stadelmann, M. (2021), expm: Matrix Exponential, Log, ’etc’.  R package version 0.999-6.
https://CRAN.R-project.org/package=expm

  * Hall & Hosseini-Nasab (2006) Hall, P. & Hosseini-Nasab, M. (2006), ‘On properties of functional principal components analysis’, Journal of the Royal Statistical Society, Series B 68(1), 109–126.
  * He et al. (2010) He, G., Müller, H. G., Wang, J. L. & Yang, W. (2010), ‘Functional linear regression via canonical analysis’, Bernoulli 16(3), 705–729.
  * Horvath & Kokoszka (2012) Horvath, L. & Kokoszka, P. (2012), Inference for Functional Data with Applications, Springer, New York.
  * Hsing & Eubank (2015) Hsing, T. & Eubank, R. (2015), Theoretical Foundations of Functional Data Analysis, with an Introduction to Linear Operators, John Wiley & Sons, Chennai, India.
  * Hubert & Branden (2003) Hubert, M. & Branden, K. V. (2003), ‘Robust methods for partial least squares regression’, Journal of Chemometrics 17(10), 537–549.
  * Hullait et al. (2021) Hullait, H., Leslie, D. S., Pavlidis, N. G. & King, S. (2021), ‘Robust function-on-function regression’, Technometrics 63(3), 396–409.
  * Hyndman & Shang (2009) Hyndman, R. J. & Shang, H. L. (2009), ‘Forecasting functional time series’, Journal of the Korean Statistical Society 38(3), 199–211.
  * Hyndman & Shang (2010) Hyndman, R. J. & Shang, H. L. (2010), ‘Rainbow plots, bagplots, and boxplots for functional data’, Journal of Computational and Graphical Statistics 19(1), 29–45.
  * Ivanescu et al. (2015) Ivanescu, A. E., Staicu, A.-M., Scheipl, F. & Greven, S. (2015), ‘Penalized function-on-function regression’, Computational Statistics 30(2), 539–568.
  * Kalogridis & Aelst (2019) Kalogridis, I. & Aelst, S. V. (2019), ‘Robust functional regression based on principal components’, Journal of Multivariate Analysis 173, 393–415.
  * Kokoszka & Reimherr (2017) Kokoszka, P. & Reimherr, M. (2017), Introduction to Functional Data Analysis, CRC Press, Boca Raton.
  * Lindsay (1994) Lindsay, B. G. (1994), ‘Efficiency versus robustness: the case for minimum hellinger distance and related methods’, Annals of Statistics 22(2), 1018–1114.
  * Luo & Qi (2019) Luo, R. & Qi, X. (2019), ‘Interaction model and model selection for function-on-function regression’, Journal of Computational and Graphical Statistics 28(2), 309–322.
  * Markatou (1996) Markatou, M. (1996), ‘Robust statistical inference: weighted likelihoods or usual M-estimation?’, Communications in Statistics - Theory and Methods 25(11), 2597–2613.
  * Maronna & Yohai (2013) Maronna, R. A. & Yohai, V. J. (2013), ‘Robust functional linear regression based on splines’, Computational Statistics and Data Analysis 65, 46–55.
  * Martines-Hernandez & Genton (2020) Martines-Hernandez, I. & Genton, M. G. (2020), ‘Recent developments in complex and spatially correlated functional data’, Brazilian Journal of Probability and Statistics 34(2), 204–229.
  * 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. & Konishi, S. (2009), ‘Regularized functional regression modeling for functional response and predictors’, Journal of Math-for-Industry 1(A3), 17–25.
  * Müller & Yao (2008) Müller, H.-G. & Yao, F. (2008), ‘Functional additive models’, Journal of the American Statistical Association, Theory and Methods 103(484), 1534–1544.
  * Olsen et al. (2018) Olsen, N. L., Markussen, B. & Raket, L. L. (2018), ‘Simultaneous inference for misaligned multivariate functional data’, Journal of the Royal Statistical Society, Series C 67(5), 1147–1176.
  * Preda & Saporta (2005) Preda, C. & Saporta, G. (2005), ‘PLS regression on a stochastic process’, Computational Statistics & Data Analysis 48(1), 149–158.
  * Preda & Schiltz (2011) Preda, C. & Schiltz, J. (2011), Functional PLS regression with functional response: the basis expansion approach, in ‘Proceedings of the 14th Applied Stochastic Models and Data Analysis Conference’, Universita di Roma La Spienza, pp. 1126–1133.
  * Raña et al. (2015) Raña, P., Aneiros, G. & Vilar, J. M. (2015), ‘Detection of outliers in functional time series’, Environmetrics 26(3), 178–191.
  * Ramsay & Dalzell (1991) Ramsay, J. O. & Dalzell, C. J. (1991), ‘Some tools for functional data analysis’, Journal of the Royal Statistical Society, Series B 53(3), 539–572.
  * Ramsay et al. (2020) Ramsay, J. O., Graves, S. & Hooker, G. (2020), fda: Functional Data Analysis.  R package version 5.1.9.
https://CRAN.R-project.org/package=fda

  * Ramsay & Silverman (2002) Ramsay, J. O. & Silverman, B. W. (2002), Applied Functional Data Analysis, Springer, New York.
  * Ramsay & Silverman (2006) Ramsay, J. O. & Silverman, B. W. (2006), Functional Data Analysis, Springer, New York.
  * Rao (2010) Rao, B. L. S. P. (2010), ‘Nonparametric density estimation for functional data by delta sequences’, Brazilian Journal of Probability and Statistics 24(3), 468–478.
  * Reiss & Odgen (2007) Reiss, P. T. & Odgen, R. T. (2007), ‘Functional principal component regression and functional partial least squares’, Journal of the American Statistical Association: Theory and Methods 102(479), 984–996.
  * Serrneels et al. (2005) Serrneels, S., Croux, C., Filzmoser, P. & Espen, P. J. V. (2005), ‘Partial robust M-regression’, Chemometrics and Intelligent Laboratory Systems 79(1-2), 55–64.
  * Shin & Lee (2016) Shin, H. & Lee, S. (2016), ‘An RKHS approach to robust functional linear regression’, Statistica Sinica 26, 255–272.
  * Tucker (1938) Tucker, R. S. (1938), ‘The reasons for price rigidity’, The American Economic Review 28(1), 41–54.
  * Wakelinc & Macfie (1992) Wakelinc, I. N. & Macfie, H. J. H. (1992), ‘A robust pls procedure’, Journal of Chemometrics 6(4), 189–198.
  * Wang (2014) Wang, W. (2014), ‘Linear mixed function‐on‐function regression models’, Biometrics 70(4), 794–801.
  * Wilcox (2012) Wilcox, R. (2012), Introduction to Robust Estimation and Hypothesis Testing, Elsevier, Waltham, MA.
  * 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.
  * Yamanishi & Tanaka (2003) Yamanishi, Y. & Tanaka, Y. (2003), ‘Geographically weighted functional multiple regression analysis: A numerical investigation’, Journal of the Japanese Society of Computational Statistics 15(2), 307–317.
  * Yao et al. (2005) Yao, F., Müller, H.-G. & Wang, J.-L. (2005), ‘Functional linear regression analysis for longitudinal data’, Annals of Statistics 33(6), 2873–2903.
  * Younse (2020) Younse, A. (2020), ‘Nonparametric discrimination of areal functional data’, Brazilian Journal of Probability and Statistics 34(1), 112–126.
  * Zhang et al. (2018) Zhang, J., Zhou, Y., Cui, X. & Xu, W. (2018), ‘Semiparametric quantile estimation for varying coefficient partially linear measurement errors models’, Brazilian Journal of Probability and Statistics 32(3), 616–656.