Introduction

The interest and need for developing statistical methods for analyzing functional data have increased in the last few decades. There have been many recent theoretical and applied developments in functional data analysis tools (see Ramsay and Dalzell 1991; Ramsay and Silverman 2002, 2006; Ferraty and Vieu 2006; Horváth and Kokoszka 2012; Cuevas 2014; Hsing and Eubank 2015; Marron et al. 2015; Srivastava and Klassen 2016; Dryden and Mardia 2016; Kokoszka and Reimherr 2017). Among many others, the scalar-on-function linear regression model (SFLRM), where the response is scalar-valued and predictor(s) consist of random functions (see Cardot, Ferraty, and Sarda 1991, 2003; James 2002; Reiss and Ogden 2007; Chen, Hall, and Müller 2011; Jiang and Wang 2011; J. Goldsmith et al. 2011; Dou, Pollard, and Zhou 2012; Tucker, Lewis, and Srivastava 2019; Ahn et al. 2020; Beyaztas and Shang 2022), and the function-on-function linear regression model (FFLRM), where both the response and predictor(s) consist of random curves (see Yao, Müller, and Wang 2005; Harezlak et al. 2007; Şentürk and Müller 2008; Matsui, Kawano, and Konishi 2009; Ivanescu et al. 2015; Scheipl, Staicu, and Greven 2015; Chiou, Yang, and Chen 2016; Beyaztas and Shang 2020), have received considerable attention among researchers as a tool for exploring the functional relationship between a scalar response-functional predictors and functional response-functional predictors, respectively. In addition, please see the available R packages fda (Ramsay, Graves, and Hooker 2022) and refund (Jeff Goldsmith et al. 2022) for the implementation of many functional data analysis methods, including SFLRM and FFLRM.

Most of the existing methods developed to estimate the SFLRM and FFLRM are non-robust to outlying observations, i.e. observations which are generated by a stochastic process with a distribution different from that of the rest of the observations (see, e.g., Raña, Aneiros, and Vilar 2015). In the case of outliers, the non-robust methods produce biased estimates; thus, predictions obtained from the fitted models become unreliable (see, e.g., Zhu, Brown, and Morris 2011; Maronna and Yohai 2013; Shin and Lee 2016; Kalogridis and Aelst 2019; Boente, Salibian-Barrera, and Vena 2020; Hullait et al. 2021; Beyaztas and Shang 2022). These methods may also lack robustness because they are coss-sectional and the fitted mean of the response variable may not be representative of the underlying data generating process. In this paper, we provide a hands-on tutorial for the implementation of functional principal component regression based on several robust approaches ( available in the R package robflreg), for robustly modeling and predicting the SFLRM and FFLRM in the presence of outliers.

The methods available in the robflreg package are centered on the robust functional principal component analysis (RFPCA) approach of (Bali et al. 2011). It uses the robust projection pursuit approach of (Croux and Ruiz-Gazen 1996) combined with a robust scale estimator to produce functional principal components and the corresponding principal component scores. With this approach, the infinite-dimensional SFLRM and FFLRM are projected onto a finite-dimensional space of RFPCA bases. Then, for the SFLRM, the robust estimation methods, including the least trimmed squares (LTS) of (Rousseeuw 1984), MM-type regression estimator (MM) of (Yohai 1987) and (Koller and Stahel 2011), S estimator, and the tau estimator of (Salibian-Barrera, Willems, and Zamar 2008), are used to estimate the parameter vector of regression model, where the scalar-valued response is predicted by the robust principal component scores of the functional predictors. For the FFLRM, on the other hand, the robust estimation methods, including the minimum covariance determinant estimator (MCD) of (Rousseeuw et al. 2004), multivariate least trimmed squares estimator (MLTS) of (Bali et al. 2008), MM estimator of (Kudraszow and Moronna 2011), S estimator of (Bilodeau and Duchesne 2000), and the tau estimator of (Ben, Martinez, and Yohai 2006), are used to estimate the parameter matrix of the regression model between the robust principal component scores of the functional response and the functional predictor variables. Besides the robust procedures, the package robflreg allows for non-robust estimation of the functional principal component regression model using the classical functional principal component analysis (FPCA) of (Ramsay and Silverman 2006) and the least-squares estimator.

The remainder of this paper is organized as follows. The SFLRM and FFLRM, as well as the techniques used for modeling and predicting these regression models, are reviewed and implemented using the robflreg package. Some ideas on how the available R package robflreg can be further improved are given at the end.

Functional linear regression models

We present the SFLRM and FFLRM, respectively.

The SFLRM

We consider a random sample \(\left\lbrace Y_i, \bm{\mathcal{X}}_{i}(s): i = 1, \ldots, n \right\rbrace\) from the pair \(( Y, \bm{\mathcal{X}} )\), where \(Y \in \mathbb{R}\) is the scalar response and \(\bm{\mathcal{X}} = [\mathcal{X}_1(s), \ldots, \mathcal{X}_P(s)]^\top\) with \(\mathcal{X}_p(s) \in \mathcal{L}_2[0,\mathcal{I}]\), \(\forall~p = 1, \ldots, P\) is the vector of \(P\) set of functional predictors whose sample elements are denoted by curves belonging to \(\mathcal{L}_2\) Hilbert space, denoted by \(\mathcal{H}\), with bounded and closed interval \(s \in \mathcal{I}\). We assume that the functional predictors \(\mathcal{X}_p(s)\), for \(p = 1, \ldots, P\), have finite second-order moments, i.e., \(\text{E}[\Vert \mathcal{X}_p(s) \Vert] < \infty\). Without loss of generality, we also assume that \(Y\) and \(\mathcal{X}_p(s)\), for \(p = 1, \ldots, P\), are mean-zero processes, so that \(\text{E}[Y] = \text{E}[\mathcal{X}_p(s)] = 0\) and \(s \in [0,1]\). Then, the SFRM is defined as follows: \[\label{eq:sof} Y_i = \int_0^1 \bm{\mathcal{X}}_i^\top(s) \bm{\beta}(s) ds + \epsilon_i,\] where \(\beta_p(s) \in \mathcal{L}_2[0,1]\) is the regression coefficient function linking \(Y\) with \(\mathcal{X}_p(s)\), and \(\bm{\beta}(s) = [ \beta_1(s), \ldots, \beta_P(s) ]^\top \in \mathcal{L}_2^P[0,1]\), and \(\epsilon_i\) is the error term which is assumed to follow a Gaussian distribution with mean-zero and variance \(\sigma^2\).

Simulation of a dataset for the SFLRM

The function generate.sf.data() in the package robflreg allows to simulate a dataset for the SFRM \[eq:sof\] as follows:

0)

Here, the argument n denotes the number of observations for each variable to be generated, n.pred denotes the number of functional predictors to be generated, n.gp denotes the number of grid points (i.e. the number of breaks in a fine grid on the interval \([0, 1]\)), and out.p is a number between 0 and 1, denoting percentage of outliers in the generated data. In the data generation process, first, generate.sf.data() simulates the functional predictors based on the following process: \[\mathcal{X}(s) = \sum_{j=1}^5 \kappa_j \nu_j(s),\] where \(\kappa_j\) is a vector generated from a Normal distribution with mean one and variance \(\sqrt{a}j^{-3/2}\), where \(a\) is is a uniformly generated random number between 1 and 4, and \[\nu_j(s) = \sin (j \pi s) - \cos (j \pi s).\] The regression coefficient functions are generated from a coefficient space that includes ten different functions, such as \(b \sin(2 \pi t)\) and \(b \cos(2 \pi t)\), where \(b\) is generated from a uniform distribution between 1 and 3. The error process is generated from the standard normal distribution. Finally, the scalar response is obtained using \[eq:sof\]. Suppose outliers are allowed in the generated data, i.e., \(out.p > 0\). Then, the randomly selected \(n \times out.p\) cases are generated differently from the process mentioned above. In more detail, if \(out.p > 0\), the regression coefficient functions (possibly different from the previously generated coefficient functions) generated from the coefficient space with \(b^*\) (instead of \(b\)), where \(b^*\) is generated from a uniform distribution between 3 and 5, are used to generate the outlying observations. Further, in this case, the following process is used to generate functional predictors: \[\mathcal{X}^*(s) = \sum_{j=1}^5 \kappa_j^* \nu_j^*(s),\] where \(\kappa_j^*\) is a vector generated from a Normal distribution with mean one and variance \(\sqrt{a}j^{-1/2}\) and \[\nu_j^*(s) = 2 \sin (j \pi s) - \cos (j \pi s).\] Moreover, the error process is generated from a normal distribution with mean of zero and a variance of one. A graphical display of the generated dataset with five functional predictors and n = 400 observations at 101 equally spaced points in the interval \([0, 1]\) obtained by generate.sf.data() is presented in Figure 1.

image image image
image image image

Plots of the simulated scalar response and the functional predictor variables. a) Observations generated under the main data generating process (black) and outlier points (grey). b-f) Curves of the functional predictors generated under the main data generating process (black) and outlier curves (grey).

Figure 1 can be produced by the following code (refer to the supplement code file for all the reproducible code):

library(robflreg)
library(fda.usc)
set.seed(202)

400
1]
model
0.1)

variable
sim.data\(Y # Predictors X <- sim.data\)X
functions
sim.data\(f.coef # Plot the scalar response out.indx <- sim.data\)out.indx
"",
range(Y))
"blue")
predictor
101))
"black",
range(fX1))
"grey")

The FFLRM

Let us consider a random sample \(\lbrace \mathcal{Y}_i(t), \bm{\mathcal{X}}_i(s)\): \(i = 1, 2, \ldots n \rbrace\) from the pair \((\mathcal{Y}, \bm{\mathcal{X}})\), where \(\mathcal{Y} \in \mathcal{L}_2 [0,1]\) is the functional response and \(\bm{\mathcal{X}} = [\mathcal{X}_1(s), \ldots, \mathcal{X}_P(s)]^\top\) with \(\mathcal{X}_p(s) \in \mathcal{L}_2[0,1]\), \(\forall~p = 1, \ldots, P\) is the vector of \(P\) set of functional predictors. We assume that the functional response and functional predictors have finite second-order moments, i.e., \(\text{E}[\Vert \mathcal{Y}(t) \Vert] = \text{E}[\Vert \mathcal{X}_p(s) \Vert] < \infty\), for \(p = 1, \ldots, P\). Without loss of generality, we also assume that both \(\mathcal{Y}(t)\) and \(\mathcal{X}_p(s)\), for \(p = 1, \ldots, P\), are mean-zero processes, so that \(\text{E}[Y(t)] = \text{E}[\mathcal{X}_p(s)] = 0\). Then, the FFRM is defined as follows: \[\label{eq:fof} Y_i(t) = \int_0^1 \bm{\mathcal{X}}_i^\top(s) \bm{\beta}(s,t) ds dt + \epsilon_i(t),\] where \(\beta_p(s,t) \in \mathcal{L}_2[0,1]\) is the bivariate regression coefficient function linking \(\mathcal{Y}(t)\) with \(\mathcal{X}_p(s)\), and \(\bm{\beta}(s,t) = [ \beta_1(s,t), \ldots, \beta_P(s,t) ]^\top \in \mathcal{L}_2^P[0,1]\), and \(\epsilon_i(t) \in \mathcal{L}_2[0,1]\) is the error term which is assumed to be independent of \(\mathcal{X}_p(s)\), for \(p = 1, \ldots, P\) and \(\text{E}[\epsilon_i(t)] = 0\).

Simulation of a dataset for the FFLRM

The function generate.ff.data() allows for the simulation of a dataset for the FFLRM as follows:

0)

Similar to the generate.sf.data(), n.pred denotes the number of functional predictors to be generated, n.curve denotes the number of observations for each functional variable to be generated, n.gp denotes the number of grid points, and out.p is an integer between 0 and 1, denoting the outlier percentage in the generated data. When generating a dataset, first, the function generate.ff.data() first simulates the functional predictors via the following process: \[\mathcal{X}(s) = \sum_{j=1}^5 \kappa_j \nu_j(s),\] where \(\kappa_j\) is a vector generated from a Normal distribution with mean one and variance \(\sqrt{a}j^{-1/2}\), where \(a\) is is a uniformly generated random number between 1 and 4, and \[\nu_j(s) = \sin (j \pi s) - \cos (j \pi s).\] The bivariate regression coefficient functions are generated from a coefficient space that includes ten different functions such as \(b \sin(2 \pi s) \sin(\pi t)\) and \(b e^{-3(s-0.5)^2} e^{-4(t-1)^2}\), where \(b\) is generated from a uniform distribution between 1 and 3. The error process \(\epsilon(t)\), on the other hand, is generated from the Ornstein-Uhlenbeck process: \[\epsilon(t) = l + [\epsilon_0(t) - l] e^{-\theta t} \sigma \int_0^t e^{-\theta (t-u)} d W_u,\] where \(l, \theta > 0, \sigma > 0\) are real constants, \(\epsilon_0(t)\) is the initial value of \(\epsilon(t)\) taken from \(W_u\), and \(W_u\) is the Wiener process. If outliers are allowed in the generated data, i.e., \(out.p > 0\), then, the randomly selected \(n \times out.p\) cases are generated differently from the process mentioned above. In more detail, if \(out.p > 0\), the regression coefficient functions (possibly different from the previously generated coefficient functions) generated from the coefficient space with \(b^*\) (instead of \(b\)), where \(b^*\) is generated from a uniform distribution between 1 and 2, are used to generate the outlying observations. In addition, in this case, the following process is used to generate functional predictors: \[\mathcal{X}^*(s) = \sum_{j=1}^5 \kappa_j^* \nu_j^*(s),\] where \(\kappa_j^*\) is a vector generated from a Normal distribution with mean one and variance \(\sqrt{a}j^{-3/2}\) and \[\nu_j^*(s) = 2 \sin (j \pi s) - \cos (j \pi s).\] A graphical display of the generated dataset with five functional predictors and n = 200 observations at 101 equally spaced points in the interval \([0, 1]\) obtained by generate.ff.data() is presented in Figure 2.

image image image
image image image

Plots of the simulated functional response and functional predictor variables: In the plots, the black curves denote the observations of the response and functional predictors generated under the smooth data generation process. The grey curves denote the outlying observations in the functional variables.

Figure 2 can be produced by the following code (refer to the supplement code file for all the reproducible code):

library(robflreg)
library(fda.usc)
set.seed(202)

200
1]
model
0.1)

variable
sim.data\(Y # Predictors X <- sim.data\)X
functions
sim.data\(f.coef # Plot the scalar response out.indx <- sim.data\)out.indx
101))
"",
range(fY))
"grey")
predictor
101))
"black",
range(fX1))
"grey")

Estimation

We first review the classical and robust FPCA methods. Then, we will focus on the robust estimation of the SFLRM and FFLRM by the functional principal component regression.

Functional principal component analysis (FPCA)

For a functional random variable \(\mathcal{X}(s)\), let us denote its covariance function by \(\mathcal{C}(s_1, s_2) = \text{Cov}[\mathcal{X}(s_1),~\mathcal{X}(s_2)]\) satisfying \(\int_0^1 \int_0^1 \mathcal{C}(s_1, s_2) ds_1 ds_2 < \infty\). Then, by Mercer’s Theorem, the following representation holds: \[\mathcal{C} = \sum_{k=1}^{\infty} \kappa_k \psi_k(s_1) \psi_k(s_2), \quad \forall s_1, s_2 \in [0,1],\] where \(\left \lbrace \psi_k(s): k = 1, 2, \ldots \right \rbrace\) are orthonormal bases of eigenfunctions in \(\mathcal{L}_2[0,1]\) corresponding to the non-negative eigenvalues \(\left \lbrace \kappa_k: k = 1, 2, \ldots \right \rbrace\) with \(\kappa_k \geq \kappa_{k+1}\). In practice, most of the variability in functional variables can be captured via a finite number of the first \(K\) eigenfunctions. Thus, the covariance function of a functional variable is estimated using a pre-determined truncation constant \(K\). In addition, the orthonormal bases of eigenfunctions are unknown in practice. Thus, they are approximated via a suitable basis expansion method such as B-spline, which is used in the robflreg package.

The RFPCA of (Bali et al. 2011) follows a similar structure as the classical FPCA, but it uses a robust scale functional instead of variance. Now let \(\Vert \alpha \Vert^2 = \langle \alpha, \alpha \rangle\) denote the norm generated by the inner product \(\langle \cdot, \cdot \rangle\). Also, let \(\mathcal{F}[\alpha]\) denote the distribution of \(\langle \alpha, \mathcal{X} \rangle\) where \(\mathcal{F}\) is the distribution of \(\mathcal{X}\). Then, for a given M-scale functional \(\sigma_M(\mathcal{F})\), the orthonormal bases of eigenfunctions defined by (Bali et al. 2011) are as follows: \[\begin{cases} \psi_k(\mathcal{F}) = \underset{\begin{subarray}{c} \Vert \alpha \Vert^2 = 1 \end{subarray}}{\arg \max}~ \sigma_M(\mathcal{F}[\alpha]), & k = 1, \\ \psi_k(\mathcal{F}) = \underset{\begin{subarray}{c} \Vert \alpha \Vert^2 = 1, \alpha \in \mathcal{B}_k \end{subarray}}{\arg \max}~ \sigma_M(\mathcal{F}[\alpha]), & k \geq 2, \end{cases}\] where \(\mathcal{B}_k = \left\lbrace \alpha \in \mathcal{L}_2[0,1]: \langle \alpha, \psi_k ( \mathcal{F} ) \rangle = 0, ~ 1 \leq k \leq \ K - 1 \right\rbrace\). The \(k\)-th largest eigenvalue is given by: \[\kappa_k(\mathcal{F}) = \sigma_M^2(\mathcal{F}[\psi_k]) = \underset{\begin{subarray}{c} \Vert \alpha \Vert^2 = 1, \alpha \in \mathcal{B}_k \end{subarray}}{\max} \sigma_M^2(\mathcal{F}[\alpha]).\]

Denote by \(\sigma_M(\mathcal{F}_n[\alpha])\) the functional for \(\sigma_M\). Let \(s^2_n: \mathcal{L}_2[0,1] \rightarrow \mathbb{R}\) denote the function of empirical M-scale functional such that \(s^2(\alpha) = \sigma_M^2(\mathcal{F}[\alpha])\). Then, the RFPCA estimates of the orthonormal bases of eigenfunctions for \(\mathcal{X}(s)\) are given by \[\begin{cases} \widehat{\psi}_k(s) = \underset{\begin{subarray}{c} \Vert \alpha \Vert^2 = 1 \end{subarray}}{\arg \max}~ s_n(\alpha), & k = 1, \\ \widehat{\psi}_k(s) = \underset{\begin{subarray}{c} \alpha \in \widehat{\mathcal{B}}_k \end{subarray}}{\arg \max}~ s_n(\alpha), & k \geq 2, \end{cases}\] where \(\widehat{\mathcal{B}}_k = \left\lbrace \alpha \in \mathcal{L}_2 [0,1]: \Vert \alpha \Vert = 1, \langle \alpha, \widehat{\psi}_k \rangle = 0, ~ \forall~ 1 \leq k \leq \ K - 1 \right\rbrace\). The corresponding eigenvalues, on the other hand, are given by \[\widehat{\kappa}_k = s^2_n (\widehat{\psi}_k), \quad k \geq 1.\]

Main RFPCA function and its arguments

The main function to obtain the robust estimates of functional principal components and the corresponding principal component scores is called getPCA():

"robust"))

In the getPCA() function, the data set is provided in the data argument as a matrix. The grid points of the functional predictors are provided in the gp argument as a vector. nbasis denotes the number of B-spline basis expansion functions used to approximate the robust functional principal components. ncomp specifies the number of functional principal components to be computed. The argument emodel denotes the method used for functional principal component decomposition. If emodel = "classical", then the classical functional principal component decomposition is performed. On the other hand, if emodel = "robust", then the RFPCA method of (Bali et al. 2011) is used to obtain the functional principal components and the corresponding principal component scores. Figure 3 presents the plot of five functional principal components computed from simulated functional data using RFPCA and nbasis = 20 B-spline basis expansion functions.

Plot of the first five functional principal components of the simulated functional data. The principal components were obtained using the RFPCA and different colors correspond to different principal components.

Figure 3 can be produced by the following code:

library(robflreg)
200
1]
model
set.seed(202)
101)
variable
sim.data$Y gpY <- seq(0, 1, length.out = 101) # grid points

Perform robust functional principal component analysis on the response variable Y

rob.fpca <- getPCA(data = Y, nbasis = 20, ncomp = 5, gp = gpY, emodel = “robust”)

Principal components

PCs <- rob.fpca$PCAcoef

1)

Robust estimation of the SFLRM

In the robust estimation of the SFRM, we first consider the principal component decomposition of the functional predictors as follows: \[\mathcal{X}_p(s) = \sum_{k=1}^{K_p} \xi_{pk} \psi_{pk}(s),\] where \(K_p\) is the truncation constant for the \(p\)-th functional predictor \(\mathcal{X}_p(s)\), \(\psi_{pk}(s)\) is the \(k\)-th eigenfunction obtained by the RFPCA of (Bali et al. 2011), and \(\xi_{pk}\) is the corresponding principal component score, given by: \[\xi_{pk} = \int_0^1 \mathcal{X}_p(s) \psi_{pk}(s) ds.\]

In practice, the eigenfunctions are approximated via a basis expansion function such as B-spline. Let \(\varphi_p(s)\) denote the B-spline basis expansion function, and \(A_p = (a_{pl})\) be an \(n \times L\)-dimensional matrix of basis expansion coefficients for the \(p\)-th functional predictor variable. In addition, let \(\bm{\varphi} = \int_0^1 \varphi_p(s) \varphi_p^\top (s) ds\) and \(\bm{\varphi}^{1/2}\) denote the \(L \times L\) dimensional matrix of inner products between the basis expansion functions and its square root, respectively. Then, the infinite-dimensional RFPCA of \(\mathcal{X}_p(s)\) is equivalent to the multivariate principal component analysis of \(A_p \bm{\varphi}^{1/2}\) and the \(k\)-th eigenfunction is given by \(\psi_{pk}(s) = \bm{\varphi}^{-1/2} v_{pk}\), where \(v_{pk}\) denotes the \(p\)-th eigenvector of the sample covariance matrix of \(A_p \bm{\varphi}^{1/2}\) (see, e.g., Ocana, Aguilera, and Escabias 2007 for more information).

If we assume that the \(p\)-th regression coefficient function \(\beta_p(s)\) admits the similar functional principal decomposition as the functional predictors as follows: \[\beta_p(s) = \sum_{k=1}^{K_p} b_{pk} \psi_{pk}(s),\] where \(b_{pk} = \int_0^1 \beta_p(s) \psi_{pk}(s) ds\). Then, the infinite-dimensional SFRM in \[eq:sof\] is approximated by the finite-dimensional regression model of scalar response on all the functional principal component scores as follows: \[Y = \sum_{p=1}^P \sum_{k=1}^{K_p} b_{pk} \xi_{pk}.\]

Main functions for the robust estimation of a SFRM and their arguments

The main function for robust estimation of SFRMs is called rob.sf.reg():

"robust"),
NULL)

In the rob.sf.reg() function, the scalar response is provided in the Y argument as an \(n \times 1\)-dimensional column vector, where \(n\) is the sample size. The functional predictors, on the other hand, are provided as a list object in the X argument. Each element of X is an \(n \times L_p\)-dimensional matrix containing the observations of \(p\)-th functional predictor \(\mathcal{X}_p(s)\), where \(L_p\) is the number of grid points for \(\mathcal{X}_p(s)\). The rob.sf.reg() function also allows for scalar predictors, which can be provided in the X.scl as an \(n \times R\)-dimensional matrix, where \(R\) denotes the number of scalar predictors. In this case, the following SFRM is considered: \[Y_i = \int_0^1 \bm{\mathcal{X}}_i^\top(s) \bm{\beta}(s) ds + \text{X.scl}_i \bm{\gamma} + \epsilon_i,\] where \(\bm{\gamma}\) denotes the vector of coefficients for the scalar predictors’ matrix. The functional principal component decomposition method is provided in the emodel argument. If emodel = "classical", then the classical functional principal component decomposition is performed to obtain principal components and the corresponding principal component scores. The coefficient vector of the regression problem of scalar response on the principal component scores is estimated via the least-squares method. If emodel = "robust", then the RFPCA of (Bali et al. 2011) is performed to obtain the principal components and the corresponding principal component scores. In this case, the method used to estimate the coefficient vector of the regression problem constructed by the scalar response and principal component scores is provided in the fmodel argument. One of the methods among LTS, MM, S, and tau can be chosen to estimate the parameter vector (by specifying, for example, fmodel = tau"). The number of B-spline basis expansion functions used to approximate the functional principal components are provided in the nbasis argument as a vector with length \(p\). Suppose nbasis = NULL, then \(\min(20, L_p / 4)\) number of B-spline basis expansion functions are used for each functional predictor. The grid points for the functional predictors are provided in the gp argument as a list object. The \(p\)-th element of gp is a vector containing the grid points of the \(p\)-th functional predictor \(\mathcal{X}_p(s)\). If gp = NULL, then \(L_p\) equally spaced time points in the interval \([0,1]\) are used for the \(p\)-th functional predictor. The number of principal components to be computed for the functional predictors are provided in the ncomp argument as a vector with length \(P\). If ncomp = NULL, then, for each functional predictor, the number whose usage results in at least 95% explained variation is used as the number of principal components.

The function get.sf.coeffs() can be used to obtain the estimated regression coefficient functions from a fitted functional principal component regression:

get.sf.coeffs(object)

In this function, the argument object is the output object obtained using the function rob.sf.reg(). The function get.sf.coeffs() produces a list object whose \(p\)-th element is a vector with length \(L_p\) containing the \(p\)-th regression coefficient function \(\beta_p(s)\).

The plots of the estimated regression coefficient functions can be obtained using the function plot_sf_coeffs():

b)

In Figure 4, the plots of the regression coefficient functions obtained from simulated data (outlier-contaminated) using RFPCA and MM estimator, as well as the classical functional principal component regression, are presented. In this function, the argument object is the output object obtained by the function get.sf.coeffs(). On the other hand, the argument b is an integer value indicating which regression parameter function to plot. It is clear from this figure that, compared with the classical functional principal component regression, the robust approach produces estimated parameter functions which are more similar to the true parameter functions.

image image image

A plot of the estimated regression coefficient functions obtained from simulated data (outlier-contaminated) using RFPCA and MM estimator and the classical functional principal component regression. Blue lines denote the true parameter functions. On the other hand, black and red lines denote the estimated parameter functions by the classical functional principal component regression and robust MM-based functional principal component regression, respectively.

The following code can be used to produce this figure and the results (refer to the supplement code file for all of the reproducible code):

library(robflreg)
400
1]
model
set.seed(202)
0.1)

functions
sim.data$f.coef[[1]]
sim.data$f.coef[[2]]
sim.data$f.coef[[3]]

variable
sim.data$Y
Predictors
sim.data$X

Xs

data
method:
gp)

data
estimator:
gp)

functions
get.sf.coeffs(classical.fit)
get.sf.coeffs(robust.fit)

function
1)
"red")
2)

Robust estimation of the FFLRM

Let us consider the functional principal decompositions of both the functional response and functional predictor variables as follows: \[\mathcal{Y}(t) = \sum_{k=1}^K \zeta_k \phi_k(t), \quad \mathcal{X}_p(s) = \sum_{j=1}^{K_p} \xi_{pj} \psi_{pj}(s),\] where \(\phi_k(t)\) and \(\psi_{pj}(s)\) respectively are the \(k\)-th and \(j\)-th eigenfunctions of \(\mathcal{Y}(t)\) and \(\mathcal{X}_p(s)\) obtained by the RFPCA and \(\zeta_k\) and \(\xi_{pj}\) are the corresponding principal component scores given by \[\zeta_k = \int_0^1 \mathcal{Y}(t) \phi_k(t) dt, \quad \xi_{pj} = \int_0^1 \mathcal{X}_p(s) \psi_{pj}(s) ds.\] If we assume that the \(p\)-th bivariate regression coefficient function \(\beta_p(s,t)\) admits the principal component decomposition with the eigenfunctions \(\phi_k(t)\) and \(\psi_{pj}(s)\) as follows: \[\beta_p(s,t) = \sum_{k=1}^K \sum_{j=1}^{K_p} b_{pkj} \phi_k(t) \psi_{pj}(s),\] where \(b_{pkj} = \int_0^1 \int_0^1 \beta_p(s,t) \phi_k(t) \psi_{pj}(s) dt ds\). Then, the infinite-dimensional FFRM in \[eq:fof\] is approximated by the finite-dimensional regression model of principal component scores of the functional response on all the functional principal component scores as follows: \[\zeta_k = \sum_{p=1}^P \sum_{j=1}^{K_p} b_{pkj} \xi_{pj}.\] Finally, the following regression model is obtained for the functional response \[\mathcal{Y}(t) = \sum_{k=1}^K \left(\sum_{p=1}^P \sum_{j=1}^{K_p} b_{pkj} \xi_{pj} \right) \phi_k(t).\]

Main functions for the robust estimation of a FFRM and their arguments

The main function to estimate the FFRM robustly is called rob.ff.reg():

"robust"),
NULL,
NULL)

In the rob.ff.reg() function, the functional response is provided via the Y argument as a matrix. On the other hand, the functional predictors are provided in the argument X as a list object. Each element of X is a matrix containing the observations of \(p\)-th functional predictor. The model type to be fitted can be chosen with model argument. If model = "full", then all of the functional predictors are used in the model. On the other hand, if model = "selected", then only the significant functional predictor variables, determined by the forward variable selection procedure of (Beyaztas and Shang 2021), are used in the model. The functional principal component decomposition method is provided via the emodel argument. If emodel = "classical", then the classical functional principal component decomposition is performed to obtain principal components and the corresponding principal component scores and the coefficient vector of the regression problem of principal component scores of the functional response on the principal component scores are estimated via the least-squares method. If emodel = "robust", then the RFPCA of (Bali et al. 2011) is performed to obtain the principal components and the corresponding principal component scores. In this case, the method used to estimate the coefficient matrix of the regression problem constructed by the principal component scores is provided in the fmodel argument. Here, one of the methods, among MCD, MLTS, MM, S, and tau estimators, can be chosen to estimate the parameter matrix. The number of B-spline basis expansion functions used to approximate the functional principal components of response and predictor variables are provided in the nbasisY and nbasisX arguments, respectively. The argument nbasisY is a numeric value while the argument nbasisX is a vector with length \(P\). Suppose nbasisY = NULL and nbasisX = NULL, then \(\min(20, L_y)\) and \(\min(20, L_p)\) B-spline basis expansion functions are used to approximate the functional principal components of functional response and \(p\)-th the functional predictor, where \(L_y\) and \(L_p\) respectively denote the number of grid points for \(\mathcal{Y}(t)\) and \(\mathcal{X}_p(s)\). The grid points for the functional response and functional predictors are provided in the gpY and gpX arguments, respectively. The argument gpY is a vector consisting of the grid points of the functional response \(\mathcal{Y}(t)\). On the other hand, the argument gpX is a list object, and its \(p\)-th element is a vector containing the grid points of the \(p\)-th functional predictor \(\mathcal{X}_p(s)\). If gpY = NULL and If gpX = NULL, then equally spaced time points in the interval \([0,1]\) are used for all the functional variables. The number of functional predictors to be computed for the functional response and functional predictors are provided in the arguments ncompY and ncompX, respectively. The argument ncompY is a numeric value while the argument ncompX is a vector with length \(P\). If ncompY = NULL and ncompX = NULL, then the number whose usage results in at least 95% explained variation is used as the number of principal components for each functional variable.

The estimated bivariate regression coefficient functions from a fitted functional principal component regression model are obtained by the get.ff.coeffs() function:

get.ff.coeffs(object)

In this function, the argument object is the output object obtained using the function rob.ff.reg(). The function get.ff.coeffs() produces a list object whose \(p\)-th element is a matrix containing the \(p\)-th bivariate regression coefficient function \(\beta_p(s,t)\).

The image plots of the estimated bivariate regression coefficient functions can be obtained using the function plot_ff_coeffs():

b)

In Figure 5, the image plots of the regression coefficient functions obtained from simulated data using RFPCA and MM estimator are presented. In this function, the argument object is the output object obtained by the function get.ff.coeffs(). The argument b is an integer value indicating which regression parameter function is to be plotted.

image image image

Image plots of the estimated bivariate regression coefficient functions obtained from simulated data using RFPCA and MM estimator.

This figure and the results can be produced by the following code (refer to the supplement code file for all the reproducible code):

library(robflreg)
200
1]
model
set.seed(202)
101)
variable
sim.data\(Y # Predictors X <- sim.data\)X

Y
Xs

data
estimator:
"robust",
gpX)

functions
get.ff.coeffs(model.fit)
function
1)

Outlier detection in the functional response

Detection of outliers in functional data is an important problem (see, e.g., Sun and Genton 2011; Arribas-Gil and Romo 2014; Dai and Genton 2018). From a fitted functional principal component regression for scalar response and scalar predictors, the robflr package with the function rob.out.detect() allows to detection of outliers in the functional response. This is achieved by applying the function depth-based outlier detection method of (Febrero-Bande, Galeano, and Gonzalez-Mantelga 2008) together with the h-modal depth proposed by (Cuaves, Febrero, and Fraiman 2007) to the estimated residual functions obtained from rob.ff.reg() to determine the outliers in the response variable. In the outlier detection algorithm, the threshold value used to identify outliers is determined by the smoothed bootstrap procedure proposed by (Febrero-Bande, Galeano, and Gonzalez-Mantelga 2008). The rob.out.detect() is as follows:

FALSE)

Herein, the argument object is an output object obtained from rob.ff.reg(). alpha, whose default value is 0.01, denotes the percentile of the distribution of the functional depth. B denotes the number of bootstrap samples (the default value is B = 200). fplot is a logical argument, if fplot = TRUE, then the outlying points flagged by the method are plotted along with the values of functional response \(\mathcal{Y}(t)\).

To show how the function rob.ff.reg() works, we simulate an outlier-contaminated dataset for the FFRM. Then, we apply the outlier detection algorithm with the classical FPCA - least squares estimator and the RFPCA - MM estimator. The plots of the functional response and detected outlying observations are presented in Figure 6. The results show that the classical method fails to flag 13 outlying curves, while the robust procedure fails to flag only two outlying curves.

image image

Plots of the functional response and detected outliers: Classical method (left panel) vs. Robust method (right panel). The detected outlying curves are denoted by black, while the non-outlying curves are denoted by grey.

The following code can produce the results and Figure 6.

library(robflreg)
200
1]
model
set.seed(202)
0.1)
sim.data$out.indx
variable
sim.data$Y
Predictors
sim.data$X

Y
Xs

least-squares
"classical",
gpX)

MM-estimator
"MM",
gpX)

function
TRUE)
188
TRUE)
140
199

outliers
sort(out.indx)
199

Prediction

We review the prediction problem for a new set of functional predictors based on a fitted functional principal component regression model.

Prediction for the SFRM

When robustly predicting the unknown values of the scalar response variable for a given new set of functional predictors (\(\mathcal{X}^*(s)\)), the principal component scores of the new set of functional predictors (\(\xi^*\)) are obtained as follows: \[\xi_{pk}^* = \int_0^1 \mathcal{X}^*_{pk}(s) \widehat{\psi}_{pk}(s) ds,\] where \(\widehat{\psi}_{pk}(s)\) is the \(k\)-th eigenfunction of the \(p\)-the functional predictor obtained by the RFPCA. Then, the predictions corresponding to the new set of functional predictors are obtained as follows: \[\widehat{Y}^* = \sum_{p=1}^P \sum_{k=1}^{K_p} \widehat{b}_{pk} \xi_{pk}^*,\] where \(\widehat{b}_{pk}\) is the estimated parameter vector obtained from the fitted model rob.sf.reg().

Main function for the robust prediction of a SFRM and its arguments

The main function for the robust prediction of a SFRM is called predict_sf_regression():

NULL)

In the function predict_sf_regression(), the argument object is an output object obtained from rob.sf.reg. The new set of functional predictors is provided in the Xnew argument as a list object whose \(p\)-th element is a matrix denoting the new observations of \(\mathcal{X}_p(s)\). Xnew must have the same length and the same structure as the input X of rob.sf.reg. If scalar predictors are used in the SFRM, then, in the prediction process, the new set of scalar predictors is provided as a matrix in the Xnew.scl argument. The argument Xnew.scl must have the same length and the same structure as the input X.scl of rob.sf.reg.

To evaluate the prediction performance of classical and robust methods, we simulate a dataset with size \(n = 400\) for the SFRM. Then, the simulated data are divided into a training sample with a size of 280 and a test sample with a size of 120. Random outliers contaminate the training sample, and both the classical and robust methods with the tau estimator are applied to the training sample to predict the values of the response variable in the test sample. To compare both methods, we compute the mean squared prediction error (MSPE): \[\text{MSPE} = \frac{1}{200} \sum_{i=1}^{200} (Y_i^* - \widehat{Y}_i^*)^2,\] where \(Y_i^*\) and \(\widehat{Y}_i^*\) denote the observed and predicted values of the scalar response in the test sample. Our results indicate that the robust method considerably outperforms the classical method. The MSPE computed from the classical method is 20.9388, while the MSPE obtained from the robust method is 1.868. The reproducible code to obtain those results is as follows:

library(robflreg)
400
1]
model
set.seed(202)
0.1)
sim.data$out.indx
variable
sim.data$Y
Predictors
sim.data$X

samples.
120)
c(1:400)[-indx.test]

Y[indx.train,]
Y[indx.test,]
list()
1:5)
X[[i]][indx.train,]
X[[i]][indx.test,]

Xs

samples
gp)
model
tau-estimator
gp)
model.classical
X.test)
model.tau
X.test)
sample
method)
method)

Prediction for the FFRM

In the robust prediction of the FFRM for a given new set of functional predictors, as in the scalar-on-function regression case, the principal component scores of the new set of functional predictors are first obtained: \[\xi_{pk}^* = \int_0^1 \mathcal{X}^*_{pk}(s) \widehat{\psi}_{pk}(s) ds,\] where \(\widehat{\psi}_{pk}(s)\) is the \(k\)-th eigenfunction of the \(p\)-the functional predictor obtained by the RFPCA. Then, the predictions of functional response (\(\widehat{\mathcal{Y}}(t)\)) corresponding to the new set of functional predictors are obtained as follows: \[\widehat{\mathcal{Y}}^*(t) = \sum_{k=1}^K \left(\sum_{p=1}^P \sum_{j=1}^{K_p} \widehat{b}_{pkj} \xi_{pj}^* \right) \widehat{\phi}_k(t),\] where \(\widehat{\phi}_k(t)\) is the \(k\)-th eigenfunction of the functional response obtained by RFPCA and \(\widehat{b}_{pkj}\) is the estimated parameter matrix obtained from the fitted model rob.ff.reg().

Main function for the robust prediction of a FFRM and its arguments

The main function for the robust prediction of a FFRM is called predict_ff_regression():

Xnew)

Here, the argument object is an output object obtained from rob.ff.reg. The new set of functional predictors is provided in the Xnew argument as a list object whose \(p\)-th element is a matrix denoting the new observations of \(\mathcal{X}_p(s)\). Xnew must have the same length and the same structure as the input X of rob.ff.reg.

We simulate a dataset with size \(n = 200\) for the FFRM to investigate and compare the prediction performance of the classical and robust methods. The simulated data are divided into a training sample with a size of 140 and a test sample with a size of 60. Random outliers contaminate the training sample, and both the classical and robust methods with MM estimator are applied to the training sample to predict the values of the response variable in the test sample. To compare both methods, we compute the following MSPE: \[\text{MSPE} = \frac{1}{100} \sum_{i=1}^{200} \Vert \mathcal{Y}_i^*(t) - \widehat{\mathcal{Y}}_i^*(t)) \Vert^2_{\mathcal{L}_2},\] where \(\mathcal{Y}_i^*(t)\) and \(\widehat{\mathcal{Y}}_i^*(t))\) denote the observed and predicted values of the functional response in the test sample. Our results show that the robust method produces a significantly smaller MSPE value than the classical method. The MSPE computed from the classical method is 3.3213, while the MSPE obtained from the robust method is 0.5925. The reproducible code to obtain those results is as follows:

library(robflreg)
200
1]
model
set.seed(202)
0.1)
sim.data$out.indx
variable
sim.data$Y
variables
sim.data$X
samples.
60)
c(1:200)[-indx.test]
Y[indx.train,]
Y[indx.test,]
list()
1:5)
X[[i]][indx.train,]
X[[i]][indx.test,]

Y
Xs

samples
"full",
gpX)
regression
MM-estimator
"robust",
gpX)
model.classical
X.test)
model.MM
X.test)
sample
method)
method)

Data analysis

We consider the MaryRiverFlow dataset available in the robflreg package to present the superiority of the robust functional principal component regression models over the classical model when the dataset includes outliers. The MaryRiverFlow dataset consists of hourly river-flow measurements obtained from January 2009 to December 2014 (6 years in total) in the Mary River, Australia. River-flow time series varies throughout the years due to the variation of the seasons and the amount of rainfall received at particular catchments. This problem still needs to be solved in the hydrological domain to be addressed appropriately using theoretical models.

Our first aim with the MaryRiverFlow dataset is to assess how the previous river flow curve time series, \(\mathcal{X}_i(s)\), affects the current day’s maximum river flow, \(Y_i\). Here, the observations of predictor are assumed to be functions of hours, i.e., there are 2188 functional observations \(\mathcal{X}_i(t) \big(1\leq t\leq 24\),\(~ i=1, \ldots, 2188 \big)\) while the response is a scalar predictor. A graphical display of the response and predictor is presented in Figure 7. From this figure, both the scalar response and functional predictor include clear outliers, which motivates us to apply the robust functional principal component regression models to better model this dataset.

image image

Plot of the daily maximum river flow measurement (left panel) and the functional time series plot (right panel) of the flow level in the Mary River.

Figure 7 can be produced by the following code:

library(robflreg)
library(fda.usc)
data("MaryRiverFlow")
MaryRiverFlow[1:2188,]
max)
by="days")
"Response")
1:24)
"black",
0))

We assume the SFLRM, \(Y_i = \beta_0 + \int \mathcal{X}_i(s) \beta(s) ds\), where \(Y_i\) is the maximum river flow measurement for the current day and \(\mathcal{X}_i(s)\) is the true river flow measurements recorded in the previous day. An expanding-window approach is considered to compare the predictive performance of the robust methods with the classical one. While doing so, The entire data are divided into two parts: a training sample containing the days from 01/01/2009 to 20/12/2014 and a test sample including days from 21/12/2014 to 30/12/2014. First, the models are constructed using the entire training data to forecast the maximum river flow measurement on 21/12/2014. Then, the maximum river flow measurement is forecasted by increasing the training samples by one day. This procedure is repeated until the training samples cover the entire dataset. The MSPE is computed for each method, and the computed mean MSPEs (\(\times 10^{-3}\)) along with standard errors (\(\times 10^{-3}\) given in bracket) are 7.559 (5.674), 1.635 (1.303), 1.671 (1.376), 1.617 (1.274), and 1.617 (1.274) for the classical, LTS, MM, S, and tau-estimator based SFLRM, respectively. From the results, all the robust methods produce significantly smaller MSPE values than the classical method. These results can be obtained by the following code:

library(robflreg)
data("MaryRiverFlow")
as.matrix(MaryRiverFlow)
list(MaryRiverFlow[1:2188,])
max)
1)

10)
c("classical","LTS","MM","S","tau")
2178

1:10)

samples
Y[1:starting_value]
Y[(starting_value+1)]

list(X[[1]][1:starting_value,])
1))

models
gp)
gp)
gp)
gp)
gp)

day
X.test)
X.test)
X.test)
X.test)
X.test)

values
pred.classical)^2
pred.LST)^2
pred.MM)^2
pred.S)^2
pred.tau)^2
1

sd)

Our second aim with the MaryRiverFlow dataset is to assess how the previous river flow curve time series, \(\mathcal{X}_i(s)\) affects the current day’s river flow curve time series, \(\mathcal{Y}_i(t)\). In this case, the elements of both the response and predictor are assumed to be functions of hours. Here, we assume the FFLRM, \(\mathcal{Y}_i(t) = \beta_0(t) + \int \mathcal{X}_i(s) \beta(s,t) ds dt\). The similar expanding-window approach used in the SFLRM example is considered, i.e., the entire dataset is divided into two parts: a training sample containing the days from 01/01/2009 to 20/12/2014 and a test sample including days from 21/12/2014 to 30/12/2014. The functional principal component regression models are constructed using the entire training data to forecast the river flow curve time series on 21/12/2014. Then, the river flow curve time series is forecasted by increasing the training samples by one day. This procedure is repeated until the training samples cover the entire dataset. The MSPE is computed for each method, and the computed mean MSPEs (\(\times 10^{-3}\)) along with standard errors (\(\times 10^{-3}\) given in bracket) are 5.721 (3.884), 1.565 (0.965), 1.377 (0.847), 1.546 (0.949), 1.511 (0.920), and 1.377 (0.849) for the classical, MCD, MLTS, MM, S, and tau-estimator based FFLRMs, respectively. From the results, the robust methods produce improved MSPEs over the classical FFLRM, i.e., the robust method models the MaryRiverFlow data better than the classical functional principal component regression model. These results can be obtained by the following code:

library(robflreg)
data("MaryRiverFlow")
as.matrix(MaryRiverFlow)

predictor
function(data,order)
dim(data)[1]
data[((order+1):n),]
x=list()
a=1
b=order
1:order)
data[(a:(n-b)),]
a+1
b-1

return(list(x=x,y=y))

predictor
Y
Xs

10)
c("classical","MCD","MLTS","MM","S","tau")
2179
decomposed.
cases.
1:10)

try(

MaryRiverFlow[1:starting_value,]
predictor
order=1)
XY\(x Y = XY\)y

models
"classical",
gpX)
"robust",
gpX)
"robust",
gpX)
"robust",
gpX)
"robust",
gpX)
"robust",
gpX)
list(t(as.matrix(Y[dim(Y)[1],])))

day
Xnew)
Xnew)
Xnew)
Xnew)
Xnew)
Xnew)

values
MaryRiverFlow[starting_value+1,])^2)
MaryRiverFlow[starting_value+1,])^2)
MaryRiverFlow[starting_value+1,])^2)
MaryRiverFlow[starting_value+1,])^2)
MaryRiverFlow[starting_value+1,])^2)
MaryRiverFlow[starting_value+1,])^2)
+1

,silent=T)

na.rm=TRUE)

Conclusion

The R package robflreg provides an implementation of functional principal component regression model based on several robust procedures to fit and predict SFLRM and FFLRM. These methods are centered on the RFPCA of (Bali et al. 2011), a popular robust dimension reduction technique in functional data, and several robust regression parameter estimators. In addition, the package robflreg allows us to fit and predict SFLRM and FFLRM via the classical FPCA and least-squares estimator. Several simulation examples and empirical data analysis show that the robust procedures provide better inference for functional linear regression models when outliers are present in the response and predictor variables. The robflreg package code can be found at: https://github.com/cran/robflreg.

The aspects that the current version of the R package robflreg that may be improved upon in the future are listed below.

  1. In the current version, the scalar predictors are allowed in the SFLRM only. In the next versions of the package, it is planned to improve the functions rob.ff.reg and predict_ff_regression to include scalar predictors (whose effects are constant and/or vary along the continuum of the functional predictor).

  2. The current package version does not allow for modeling the function-on-scalar regression model, where the response consists of random functions, and the predictors include scalar observations. The robust procedures used in the FFLRM are planned to extend the function-on-scalar regression model in the subsequent versions of the package.

  3. In the current version, the observations of the functional variables are assumed to be densely observed. In the next versions of the package, the robust methods are planned to extend the models where the elements of functional variables can also be observed over irregular and curve-specific grids.

Acknowledgement

We thank three reviewers for their careful reading of our manuscript and valuable suggestions and comments, which have helped us produce an improved version of our manuscript and the R package. The first author was supported by The Scientific and Technological Research Council of Turkey (TUBITAK) (grant no: 120F270). The second author was partially supported by an Australian Research Council Discovery Project (grant no: DP230102250). This study is dedicated to the people who lost their lives in the earthquake that occurred in Turkey on February 6, 2023.

Ahn, M. K., J. D. Tucker, W. Wu, and A. Srivastava. 2020. “Regression Models Using Shapes of Functions as Predictors.” Computational Statistics and Data Analysis 151: 107017. https://doi.org/https://doi.org/10.1016/j.csda.2020.107017.
Arribas-Gil, A., and J. Romo. 2014. “Shape Outlier Detection and Visualization for Functional Data: The Outliergram.” Biostatistics 15 (4): 603–19. https://doi.org/https://doi.org/10.1093/biostatistics/kxu006.
Bali, J. L., G. Boente, D. E. Tyler, and J-L. Wang. 2008. “The Multivariate Least-Trimmed Squares Estimator.” Journal of Multivariate Analysis 99 (3): 311–38. https://doi.org/https://doi.org/10.1016/j.jmva.2006.06.005.
———. 2011. “Robust Functional Principal Components: A Projection-Pursuit Approach.” The Annals of Statistics 39 (6): 2852–82. https://doi.org/https://doi.org/10.1214/11-AOS923.
Ben, M. G., E. Martinez, and V. J. Yohai. 2006. “Robust Estimation for the Multivariate Linear Model Based on a \(\tau\) Scale.” Journal of Multivariate Analysis 97 (7): 1600–1622. https://doi.org/https://doi.org/10.1016/j.jmva.2005.08.007.
Beyaztas, U., and H. L. Shang. 2020. “On Function-on-Function Regression: Partial Least Squares Approach.” Environmental and Ecological Statistics 27 (1): 95–114. https://doi.org/https://doi.org/10.1007/s10651-019-00436-1.
———. 2021. “A partial least squares approach for function-on-function interaction regression.” Computational Statistics 36 (2): 911–39. https://doi.org/https://doi.org/10.1080/03610926.2022.2065018.
———. 2022. “A Robust Functional Partial Least Squares for Scalar-on-Multiple-Function Regression.” Journal of Chemometrics 36 (4): e3394. https://doi.org/https://doi.org/10.1002/cem.3394.
Bilodeau, M., and P. Duchesne. 2000. “Robust Estimation of the SUR Model.” The Canadian Journal of Statistics 28 (2): 277–88. https://doi.org/https://doi.org/10.2307/3315978.
Boente, G., M. Salibian-Barrera, and P. Vena. 2020. “Robust Estimation for Semi-Functional Linear Regression Models.” Computational Statistics & Data Analysis 152: 107041. https://doi.org/https://doi.org/10.1016/j.csda.2020.107041.
Cardot, H., F. Ferraty, and P. Sarda. 1991. “Functional Linear Model.” Statistics and Probability Letters 45 (1): 11–22. https://doi.org/https://doi.org/10.1016/S0167-7152(99)00036-X.
———. 2003. “Spline Estimators for the Functional Linearmodel.” Statistica Sinica 13 (3): 571–91.
Chen, D., P. Hall, and H-G. Müller. 2011. “Single and Multiple Index Functional Regression Models with Nonparametric Link.” The Annals of Statistics 39 (3): 1720–47. https://doi.org/https://doi.org/10.1214/11-AOS882.
Chiou, J. M., Y. F. Yang, and Y. T. Chen. 2016. “Multivariate Functional Linear Regression and Prediction.” Journal of Multivariate Analysis 146: 301–12. https://doi.org/https://doi.org/10.1016/j.jmva.2015.10.003.
Croux, C., and A. Ruiz-Gazen. 1996. “High Breakdown Estimators for Principal Components: The Projection-Pursuit Approach Revisited.” Journal of Multivariate Analysis 95 (1): 206–26. https://doi.org/https://doi.org/10.1016/j.jmva.2004.08.002.
Cuaves, A., M. Febrero, and R. Fraiman. 2007. “Robust Estimation and Classification for Functional Data via Projection-Based Depth Notions.” Computational Statistics 22 (3): 481–96. https://doi.org/https://doi.org/10.1007/s00180-007-0053-0.
Cuevas, A. 2014. “A Partial Overview of the Theory of Statistics with Functional Data.” Journal of Statistical Planning and Inference 147: 1–23. https://doi.org/https://doi.org/10.1016/j.jspi.2013.04.002.
Dai, W., and M. G. Genton. 2018. “Multivariate Functional Data Visualization and Outlier Detection.” Journal of Computational and Graphical Statistics 27 (4): 923–34. https://doi.org/https://doi.org/10.1080/10618600.2018.1473781.
Dou, W. W., D. Pollard, and H. H. Zhou. 2012. “Estimation in Functional Regression for General Exponential Families.” The Annals of Statistics 40 (5): 2421–51. https://doi.org/https://doi.org/10.1214/12-AOS1027.
Dryden, I. L., and K. V. Mardia. 2016. Statistical Shape Analysis, with Applications in R. New York: Wiley Series in Probability; Statistics.
Febrero-Bande, M., P. Galeano, and W. Gonzalez-Mantelga. 2008. “Outlier detection in functional data by depth measures, with application to identify abnormal NO\(_x\) levels.” Environmetrics 19 (4): 331–45. https://doi.org/https://doi.org/10.1002/env.878.
Ferraty, F., and P. Vieu. 2006. Nonparametric Functional Data Analysis. New York: Springer.
Goldsmith, J., J. Bobb, C. M. Crainiceanu, B. Caffo, and D. Reich. 2011. “Penalized Functional Regression.” Journal of Computational and Graphical Statistics 20 (4): 830–51. https://doi.org/https://doi.org/10.1198/jcgs.2010.10007.
Goldsmith, Jeff, Fabian Scheipl, Lei Huang, Julia Wrobel, Chongzhi Di, Jonathan Gellar, Jaroslaw Harezlak, et al. 2022. Refund: Regression with Functional Data. https://CRAN.R-project.org/package=refund.
Harezlak, J., B. A. Coull, N. M. Laird, S. R. Magari, and D. C. Christiani. 2007. “Penalized Solutions to Functional Regression Problems.” Computational Statistics and Data Analysiss 51 (10): 4911–25. https://doi.org/https://doi.org/10.1016/j.csda.2006.09.034.
Horváth, L., and P. Kokoszka. 2012. Inference for Functional Data with Applications. New York: Springer.
Hsing, T., and R. Eubank. 2015. Theoretical Foundations of Functional Data Analysis, with an Introduction to Linear Operators. Chennai, India: John Wiley & Sons.
Hullait, H., D. S. Leslie, N. G. Pavlidis, and S. King. 2021. “Robust Function-on-Function Regression.” Technometrics 63 (3): 396–409. https://doi.org/https://doi.org/10.1080/00401706.2020.1802350.
Ivanescu, A. E., A. M. Staicu, F. Scheipl, and S. Greven. 2015. “Penalized Function-on-Function Regression.” Computational Statistics 30 (2): 539–68. https://doi.org/https://doi.org/10.1007/s00180-014-0548-4.
James, G. M. 2002. “Generalized Linear Models with Functional Predictors.” Journal of Royal Statistical Society, Series B 64 (3): 411–32. https://doi.org/https://doi.org/10.1111/1467-9868.00342.
Jiang, C., and J. L. Wang. 2011. “Functional Single Index Models for Longitudinal Data.” The Annals of Statistics 39 (1): 362–88. https://doi.org/https://doi.org/10.1214/10-AOS845.
Kalogridis, I., and S. Van Aelst. 2019. “Robust Functional Regression Based on Principal Components.” Journal of Multivariate Analysis 173: 393–415. https://doi.org/https://doi.org/10.1016/j.jmva.2019.04.003.
Kokoszka, P., and M. Reimherr. 2017. Introduction to Functional Data Analysis. Boca Raton: CRC Press.
Koller, M., and W. A. Stahel. 2011. “Sharpening Wald-Type Inference in Robust Regression for Small Samples.” Computational Statistics & Data Analysis 55 (8): 2504–15. https://doi.org/https://doi.org/10.1016/j.csda.2011.02.014.
Kudraszow, N. L., and R. A. Moronna. 2011. “Estimates of MM Type for the Multivariate Linear Model.” Journal of Multivariate Analysis 102 (9): 1280–92. https://doi.org/https://doi.org/10.1016/j.jmva.2011.04.011.
Maronna, R. A., and V. J. Yohai. 2013. “Robust Functional Linear Regression Based on Splines.” Computational Statistics and Data Analysis 65: 46–55. https://doi.org/https://doi.org/10.1016/j.csda.2011.11.014.
Marron, J. S., J. O. Ramsay, L. M. Sangalli, and A. Srivastava. 2015. “Functional Data Analysis of Amplitude and Phase Variation.” Statistical Science 30 (4): 468–84. https://doi.org/https://doi.org/10.1214/15-STS524.
Matsui, H., S. Kawano, and S. Konishi. 2009. “Regularized Functional Regression Modeling for Functional Response and Predictors.” Journal of Math-for-Industry 1 (A3): 17–25.
Ocana, F. A., A. M. Aguilera, and M. Escabias. 2007. “Computational Considerations in Functional Principal Component Analysis.” Computational Statistics 22 (3): 449–65. https://doi.org/https://doi.org/10.1007/s00180-007-0051-2.
Ramsay, J. O., and C. J. Dalzell. 1991. “Some Tools for Functional Data Analysis.” Journal of the Royal Statistical Society, Series B 53 (3): 539–72. https://doi.org/https://doi.org/10.1111/j.2517-6161.1991.tb01844.x.
Ramsay, J. O., Spencer Graves, and Giles Hooker. 2022. Fda: Functional Data Analysis. https://CRAN.R-project.org/package=fda.
Ramsay, J. O., and B. W. Silverman. 2002. Applied Functional Data Analysis. New York: Springer.
———. 2006. Functional Data Analysis. New York: Springer.
Raña, P., G. Aneiros, and J. M. Vilar. 2015. “Detection of Outliers in Functional Time Series.” Environmetrics 26 (3): 178–91. https://doi.org/https://doi.org/10.1002/env.2327.
Reiss, P. T., and T. R. Ogden. 2007. “Functional Principal Component Regression and Functional Partial Least Squares.” Journal of the American Statistical Association: Theory and Methods 102 (479): 984–96. https://doi.org/https://doi.org/10.1198/016214507000000527.
Rousseeuw, P. J. 1984. “Least Median of Squares Regression.” Journal of the American Statistical Association: Theory and Methods 79 (388): 871–81. https://doi.org/https://doi.org/10.1080/01621459.1984.10477105.
Rousseeuw, P. J., K. V. Driessen, S. V. Aelst, and J. Agullo. 2004. “Robust Multivariate Regression.” Technometrics 46 (3): 293–305. https://doi.org/https://doi.org/10.1198/004017004000000329.
Salibian-Barrera, M., G. Willems, and R. Zamar. 2008. “The Fast-Tau Estimator for Regression.” Journal of Computational and Graphical Statistics 17 (3): 659–82. https://doi.org/https://doi.org/10.1198/106186008X343785.
Scheipl, F., A-M. Staicu, and S. Greven. 2015. “Functional Additive Mixed Models.” Journal of Computational and Graphical Statistics 24 (2): 477–501. https://doi.org/https://doi.org/10.1080/10618600.2014.901914.
Şentürk, D., and H-G. Müller. 2008. “Generalized Varying Coefficient Models for Longitudinal Data.” Biometrika 95 (3): 653–66. https://doi.org/https://doi.org/10.1093/biomet/asn006.
Shin, H., and S. Lee. 2016. “An RKHS Approach to Robust Functional Linear Regression.” Statistica Sinica 26: 255–72.
Srivastava, A., and E. P. Klassen. 2016. Functional and Shape Data Analysis. New York: Springer.
Sun, Y., and M. G. Genton. 2011. “Functional Boxplots.” Journal of Computational and Graphical Statistics 20 (2): 316–34. https://doi.org/https://doi.org/10.1198/jcgs.2011.09224.
Tucker, J. D., J. R. Lewis, and A. Srivastava. 2019. “Elastic Functional Principal Component Regression.” Statistical Analysis and Data Mining 12 (2): 101–15. https://doi.org/https://doi.org/10.1002/sam.11399.
Yao, F., H-G. Müller, and J-L. Wang. 2005. “Functional Linear Regression Analysis for Longitudinal Data.” The Annals of Statistics 33 (6): 2873–903. https://doi.org/https://doi.org/10.1214/009053605000000660.
Yohai, V. J. 1987. “High Breakdown-Point and High Efficiency Estimates for Regression.” The Annals of Statistics 15 (2): 642–65. https://doi.org/https://doi.org/10.1214/aos/1176350366.
Zhu, H., P. J. Brown, and J. S. Morris. 2011. “Robust, Adaptive Functional Regression in Functional Mixed Model Framework.” Journal of the American Statistical Association 106 (1): 1167–79. https://doi.org/https://doi.org/10.1198/jasa.2011.tm10370.