Jean Steve

Tamo Tchomgui

Julien Jacques
email: julien.jacques@univ-lyon2.fr

Vincent Barriac
email: vincent.barriac@orange.com

Guillaume Fraysse
email: guillaume.fraysse@orange.com

Stéphane Chrétien
email: stephane.chretien@univ-lyon2.fr

A Penalized Spline Estimator for Functional Linear Regression with Functional Response June 1, 2023

Keywords: Functional data analysis, function-on-function regression, penalized splines, Canadian weather data

Many scientific studies in recent years have been collecting data at a high frequency, which can be considered as functional data. When both the response variable to be modelled and the covariates are functions, we provide a novel and easy-to-implement method addressing function-on-function linear modelling and obtain interpretable parameters. Two main types of models are considered: the concurrent model which explains the response curve Y i (t) at time t from the values at same time t of the covariates X l i (t); the (feed-forward) integral model which explains Y i (t) based on the values of covariate curves X l i (s) observed at any times s ≤ t. A regularized inference approach is proposed, which accurately selects an appropriate set of basis functions that can be used for functional data reconstruction and at the same time provides smooth and interpretable functional parameters. Numerical studies on simulated data with different scenarios illustrate the good performance our the method to capture the relationship between covariates and response. The method is finally applied to the well-known data in order to compare it to some existed competitors. On Canadian weather data with the problem of predicting precipitations from temperature measurements and on Hawaii Ocean data for predicting ocean salinity from temperature, oxygen, chloropigments and density measurements, our method made significant improvements on prediction error.

Introduction

In the last few decades, data sets collected with new, fast and sometimes accurate sensors have become very common in various fields of applied sciences including economics, finance, geosciences, medicine, etc. As a result, new tools are needed to process and analyze this very fast growing resource of data. A quite natural idea that has emerged lately is to extend classical tools from data analysis to a new paradigm called Functional Data Analysis (FDA). This new paradigm has proved very successful at addressing the statistical analysis of data where at least one of the variables of interest need to be treated as a function. Extension of linear regression to the functional setting has therefore naturally become a major area of research in FDA. While the literature is too vast to cover here, the recommended references for this field are [START_REF] Ramsay | Functional Data Analysis[END_REF], [START_REF] Ramsay | Functional Data Analysis with R and MATLAB[END_REF], [START_REF] Horváth | Inference for Functional Data with Applications[END_REF], [START_REF] Kokoszka | Introduction to Functional Data Analysis[END_REF], which provide excellent introductions to FDA. Moreover [START_REF] Goldsmith | Penalized functional regression[END_REF] and [START_REF] Morris | Functional regression[END_REF] provide a broad overview of the methods of functional linear regression. In the functional setting, different types of functional linear regression have been considered, depending on the functional nature of the response and/or at least one of the covariates. Thus, using the convention that first term denotes response-type and second term denotes covariate-type, the following regression models are all the possible options to consider: function-on-scalar, scalar-on-function and function-onfunction. The scalar-on-function linear regression models is the most thoroughly studied model among the three models in the current literature. Some references include [START_REF] Cardot | Functional linear model[END_REF] and [START_REF] Hastie | Varying-coefficient models[END_REF].

Most of the inference approaches for these models rely on a basis expansion assumption. For instance [START_REF] Besse | Approximation spline de la prevision d'un processus fonctionnel autorégressif d'ordre 1[END_REF] and [START_REF] Ramsay | Functional Data Analysis[END_REF] proposed spline-type approximations of the functional covariates and then performed the estimation step by minimizing a least squares criterion. Among other useful references, Antoch et al. ( 2010) uses B-spline expansions for both the functional parameters and the functional covariates. The issue of possible non-identifiability was pursued in Scheipl and Greven (2016). In theses approaches, the functional regression models become equivalent to a multivariate model on the basis expansion coefficients.

An alternative way is to consider Functional Principal Components Analysis (FPCA, [START_REF] Ramsay | Functional Data Analysis[END_REF]), possibly using smoothness promoting penalization (Silverman, 1996;[START_REF] Besse | Simultaneous non-parametric regressions of unbalanced longitudinal data[END_REF]. Possible issues in determining the number of components to account for that seem to be still open. Indeed, it was shown in [START_REF] Crainiceanu | Generalized multilevel functional regression[END_REF] that the shape of the functional parameters can drastically change as one or two additional principal components are included, making the process quite unstable and relatively difficult to interpret.

In comparison with scalar-on-function problems, function-on-function models, that we address here, have been much less studied in the literature. For instance, [START_REF] Ivanescu | Penalized function-on-function regression[END_REF] proposes to estimate a function-on-function regression model using a penalized mixed model. In this setting as well, the main issue faced is not only the problem of accurately selecting the number of basis functions and the location of the knots [START_REF] Li | On the asymptotics of penalized splines[END_REF], but also the possible interpretability of the obtained estimators [START_REF] James | Functional linear regression that's interpretable[END_REF]. Signal compression approach (wSigcomp) designed by [START_REF] Luo | Functional wavelet regression for linear function-on-function models[END_REF] which is another way to address function-on-function models firstly apply wavelets transformation to covariates and with the functional response and the obtained multivariate covariates, proposed a method to estimate the functional bivariate parameter by characterize it as the solution of a generalized functional eigenvalue problem. The Optimal Penalized Functionon-Function Regression (OPFFR) proposed by (Sun et al., 2018), produce an estimator of the 2D functional parameter as optimizer of a form of penalized least squares where the penalty enforces a certain level of smoothness.

In mathematical terms, the problem considered in the present paper is the one of estimating a linear relationship between functional covariates and functional response based on the n-sample Y i (t), X i (t) = X 1 i (t), . . . , X p i (t) ⊤ , t ∈ [0, T] i = 1, . . . , n, where the output variable Y(t) and the p input variables (X l (t)) 1≤l≤p are assumed to belong to the separable Hilbert L 2 [0; T] . In the sequel, we focus in particular on the following two functional linear models:

Y i (t) = β 0 (t) + p l=1 β l (t) X l i (t) + ε i (t) = 1, X i (t) ⊤ β(t) + ε i (t) , (1)
Y i (t) = γ 0 (t) + p l=1 t 0 γ l (s, t) X l i (s) ds + ε i (t) = γ 0 (t) + t 0 X i (s) ⊤ γ(s, t) ds + ε i (t) (2)
where β(t) = β 0 (t), β 1 (t), . . . , β p (t) ⊤ , γ(s, t) = γ 0 (t), γ 1 (s, t), γ 2 (s, t), . . . , γ p (s, t)

⊤ are the unknown functional parameters and are assumed to be square integrable; ε i (t) is the model error and is a sample of centered random variables with variance σ 2 i , specific to the i th individual [START_REF] Ramsay | Functional Data Analysis[END_REF], Chapter 13); ε i (t) and X i (t) are assumed to be uncorrelated. The noise functions ε i (t) can be rigorously defined using white noise theory as presented in [START_REF] Hida | White Noise: An Infinite Dimensional Calculus[END_REF]. In our context, we will only use the fact that when sampled at various times from a finite set T , the vector (ε i (t)) t∈T can be expressed as a sum of a vector with i.i.d. components and a vector with prescribed covariance matrix, i.e. a vector with constant components in the simplest case.

Model (1), known as the "concurrent model" , assumes that the response function at time t, Y i (t), is explained by covariate functions X l i (t), at exactly the same time t, the functional parameters being allowed to vary with t as well. The second model (2), called the "integral model", represents Y i (t) using the values of the covariates curves X l i (s) for all the observed times s ≤ t. Clearly, Model (2) is more general and richer than Model (1). Exploring the "concurrent model" further at the first step is of great interest because, as mentioned in [START_REF] Hastie | Varying-coefficient models[END_REF], any functional linear model can be reduced to this form.

In the present paper, we develop an efficient approach for estimating the functional parameters β(t) of the concurrent model (1) and γ(s, t) of the integral model (2). For this purpose, we use cubic B-spline basis expansion for both functional covariates and functional parameters. We propose penalized estimator of the corresponding functional basis coefficients. As will be shown in the sequel, our approach allows to simply choose equispaced knots and a sufficient number of basis functions to capture the main features of the covariates. Overfitting will be naturally avoided by penalizing roughness via controlling the second derivatives of the functional parameters which are being maximized.

Plan of the paper. The paper is organized as follows: Section 2 presents the two types of models we focus on and Section 3 the estimation scheme. Our estimation scheme consists of two steps. The first one addresses recovering of the functional nature of the covariates, by approximating them into a functional basis. The second step consists of penalized estimation of the functional regression coefficients, which are themselves decomposed in another functional basis.

Section 4 contains a simulation based exploration of the method which confirms the efficiency of the proposed approach. Section 5 finally presents an illustration of the method on two real data sets. The first one is the well-known Canadian weather data set, in which the goal is to explain the precipitation as a function of the temperatures in different Canadian cities. The second one is the Hawaii ocean data set in which salinity is explained as a function of four functional covariates.

Finally, Section 6 concludes the paper.

Linear models for function-on-function regression

In this section, it is shown how the functional models (1) and ( 2) can be, under the basis expansion assumption of covariates and parameters, reduce to a linear mixed model onto the discrete observations of the functional response and functional covariates.

Functional concurrent model

Linear regression for a functional response involving one or more functional covariates in the concurrent model is a well-known problem. The main issue is to estimate an infinite dimensional parameter β(t) through a finite sample of observations. As shown in [START_REF] Hastie | Varying-coefficient models[END_REF], Model (1), also called the varying coefficient model, is interesting because any functional model can be reduced to this form. Chapter 14 in [START_REF] Ramsay | Functional Data Analysis[END_REF] describes how this model can be fitted by minimizing an unweighted least squares criterion. The method proposed in this paper addresses the estimation problem using a penalized function-on-function regression as proposed in [START_REF] Ivanescu | Penalized function-on-function regression[END_REF], where the problem is represented as a mixed model. Nevertheless, our work differs by the choice of the penalization criterion enforced on the functional parameter. The parameter β(t) is expanded in functional basis using q β basis functions to get back to a classical mixed model for which the estimations of the parameters are well known. Furthermore, we allow to choose the number of basis functions q β to be large enough to capture any desired variations of β(t), and we add a roughness penalty term to get a smooth solution for the parameter at the end.

As a first step of our modelling, we recover the underlying functional process, by using penalized cubic B-splines expansion for all the functional covariates.

Functional basis expansion of covariates and model parameters

In practice, we do not properly observe a continuous curve for each realization of both the response variable Y i (t) and the covariate variables X l i (t)

1≤l≤p

. In indeed, as opposed to the ideal observation setting, we only have access to a set of noisy observations at a finite number of points on a grid. As a result, the functional data can be presented as a numerical vector. In order to recover the continuous form, which generally belongs to an infinite dimensional space (e.g. Hilbert separable space L 2 ([0, T ])), one efficient way to proceed is by expanding the considered functions in a functional basis. The functional response, which is assumed in model (1) even in model (2) to be written as a linear combination of these predictors, is not necessary to be pre-processed. The advantage of this approach is the fact that by truncating the series at a given level q l , we obtain an approximation of the covariate function X l i (t) in a q l dimensional space. So for all the p covariates X l (t), we can therefore recover a representation in cubic B-splines functional basis. As indicated by [START_REF] Li | On the asymptotics of penalized splines[END_REF], the choice of the number of knots depends on the complexity of the variable and should be large enough to capture the patterns of the variable.

It is reasonable to suppose that this number and, thus, the number of basis functions depends on the covariate. So to distinguish the basis functions of each covariate, although they just differ by their number, we will adopt in the rest of this article the system B l 1 (t), B l 2 (t), . . . , B l q X l (t) as the basis function of X l (t). Then, any functional covariate can be written as:

X l i (t) = q x l j=1 x l ij B l j (t) = B l (t) ⊤ x l i with 1 ≤ l ≤ p. (3)
The basis functions B l j (t) being prescribed, the estimation of coefficients x l ij is done as a preliminary step [START_REF] Li | On the asymptotics of penalized splines[END_REF][START_REF] Ruppert | Selecting the number of knots for penalized splines[END_REF][START_REF] Ramsay | Functional Data Analysis[END_REF].

Similarly as for functional covariates, we expand all the functional parameters β l (t) l of the concurrent model in functional basis. The number of basis functions q β l must be chosen as sufficiently large to capture the patterns of any β l (t):

β l (t) = q β l j=1 b l j ϕ l j (t) = ϕ l (t) ⊤ b l with 0 ≤ l ≤ p. (4)
Using the expressions (3) and ( 4), the components in Model (1) become:

β(t) =         β 0 (t)
β 1 (t) . . .

β p (t)         =         ϕ 0 (t) ⊤ b 0 ϕ 1 (t) ⊤ b 1 . . . ϕ p (t) ⊤ b p         =         ϕ 0 (t) ⊤ 0 . . . 0 0 ϕ 1 (t) ⊤ . . . 0 . . . . . . . . . . . . 0 0 . . . ϕ p (t) ⊤         (p+1, l q β l ) -matrix         b 0 b 1 . . . b p         l q β l -vect. = Φ(t) b ,
and

X i (t) =         1 X 1 i (t) . . . X p i (t)         =         1 B 1 (t) ⊤ x 1 i . . . B p (t) ⊤ x p i         =         1 0 . . . 0 0 B 1 (t) ⊤ . . . 0 . . . . . . . . . . . . 0 0 . . . B p (t) ⊤         (p+1, l q X l ) -matrix         1 x 1 i . . . x p i         l q X l -vect. = B(t) x i .
By plugging-in these expressions into Model (1), we get:

Y i (t) = x ⊤ i B(t) ⊤ Φ(t) b + ε i (t) = R i (t) ⊤ b + ε i (t) (5) with R i (t) = Φ(t) ⊤ B(t)
x i which is used as design matrix and b the unknown parameters to be estimated.

Functional concurrent model on the observations

The concurrent model implicitly assumes that the functional covariates and the functional response are observed at the same timestamps. The observation grid will consist of m points {t 1 , . . . , t m }.

In mathematical terms we have:

Y i (t j ) = R i (t j ) ⊤ b + ε i (t j ) with 1 ≤ i ≤ n and 1 ≤ j ≤ m. (6)
One very specific issue to take care of is that the successive values of the observation noise ε i (t 1 ), . . . , ε i (t m ) can not be assumed independent.

One way to address the question of dependency is to use a linear mixed model (LMM, Wood (2006)). We thus assume that the model error can be decomposed as

ε i (t j ) = U i + η ij , with η ij
a Gaussian white noise and U i a random variable which takes into account the random effect in each individual i = 1, . . . , n. To summarize, our model consists of a LMM with fixed effects b and random effect U i . In matrix form we get:

Y = R ⊤ b + Z U + η, (7)
where

Y = Y 1 (t 1 ), . . . , Y 1 (t m ), Y 2 (t 1 ), . . . , Y n (t m ) ⊤ , R = R i (t j ) i,j
the design matrix of dimension q β ×nm with

q β = l q β l , U = U 1 , U 2 , • • • , U n ⊤ ∼ N (0, Γ), η = (η ij ) i,j ∼ N (0, σ 2 I nm ) and Z =         1 m×1 0 m×1 . . . 0 m×1 0 m×1 1 m×1 . . . 0 m×1 . . . . . . . . . . . . 0 m×1 0 m×1 . . . 1 m×1         (nm×n) -matrix
.

The specific notations we used are the matrices 0 k×l and 1 k×l of size k × l, which are composed of zeros and ones, respectively; The notations 0 refers to the corresponding null vector and Γ the unknown covariance matrix of the random effects.

The parameters are then the fixed effects vectors b and the variance components σ 2 and Γ. We describe how to perform the inference in Section 3.

Functional integral model

The integral Model (2) assumes cumulative effects of covariates. More clearly, the model we proposes use observations of covariates up until time t to predict the response at time t. It is important to note that in most models found in the literature [START_REF] Ramsay | Functional Data Analysis[END_REF][START_REF] Horváth | Inference for Functional Data with Applications[END_REF], when both covariates and response have the same domain, consider that the response at any time t depends on the influence of the covariates on the whole domain.

Such model implicitly assumes that the covariates at any time t + s can influence the response variable at time t. However, in the integral model, the functional parameters are bivariate functions γ l (s, t), except for the constant of the model, which remains univariate. In this section, we start by expanding the parameters in a finite-dimensional functional basis and then plug this expression into the model.

Functional basis expansion of covariates and model parameters

The functional parameters are therefore expanded in a bivariate basis which may or may not have the same number of basis functions on each of the two dimensions. Without loss of generality and for the sake of simplicity, we assume that the number of basis functions is the same in the two dimensions. This leads to the following expression:

γ l (t, s) = q γ l j,k=1 a l jk B l 1j (t) B l 2k (s) (8) where B l 1j (t) 1≤j≤q γ l and B l 2j (t) 1≤j≤q γ l
are the basis functions and (a l jk ) 1≤j,k≤q γ l the unknown basis coefficients to be estimated. We can rewrite this expression in matrix form by:

γ l (t, s) = a l ⊤ B l 1 (t) B l 2 (s) (9)
with

a l = a l 11 . . . a l 1q γ l a l 21 . . . . . . a l q γ l 1 . . . a l q γ l q γ l ⊤ , B l 1 (t) = diag B l 11 (t), . . . , B l 1q γ l (t), . . . , . . . , B l 11 (t), . . . , B l 1q γ l (t) , B l 2 (s) = B l 21 (s) . . . B l 21 (s) . . . . . . B l 2q γ l (s) . . . B l 2q γ l (s) ⊤ .
The functional constant being univariate, it can thus be written as in (4) in the form:

γ 0 (t) = q γ 0 j=1 a 0 j B 0 j (t) = B 0 (t) ⊤ a 0 .

Functional integral model on the observations

By plugging covariates and parameters functional basis expansion in the integral Model (2), we get:

Y i (t) = γ 0 (t) + p l=1 t 0 x l ⊤ i B l (s) B l 2 (s) ⊤ B l 1 (t) ⊤ a l ds + ε i (t) = γ 0 (t) + p l=1 x l ⊤ i t 0 B l (s) B l 2 (s) ⊤ ds B l 2 (t) B l 1 (t) ⊤ a l + ε i (t) = γ 0 (t) + p l=1 x l ⊤ i B l 2 (t) B l 1 (t) a l + ε i (t) = B 0 (t) ⊤ a 0 + p l=1 Q l i (t) ⊤ a l + ε i (t) , with Q l i (t) = B l 1 (t) ⊤ B l 2 (t) ⊤ x l i .
Finally we obtain:

Y i (t) = Q i (t) ⊤ a + ε i (t) (10)
with a = (a 0 , a 1 , a 2 , . . . , a p ) ⊤ and Q

i (t) = B 0 (t) ⊤ , Q 1 i (t) ⊤ , Q 2 i (t) ⊤ , . . . , Q p i (t) ⊤ ⊤ two vectors of length q γ = q γ 0 + p l=1
q 2 γ l . Once again, we are faced with the problem of lack of independence of the different measured values for the same individual. We will proceed exactly in the same way as with the concurrent model using a linear mixed model with fixed effects given by the vector a and random effects given by the random vector U = (U i ) i . The model will therefore be written as a LMM given by:

Y = Q ⊤ a + Z U + η, (11)
with Z, U and η define similarly to (7

). Q = Q i (t j ) i,j
the design matrix of dimension q γ × nm.

As in the concurrent model, the parameters we need to estimate are the fixed effects vectors a and the variance components σ 2 and Γ. The inference scheme is described in Section 3.

B-spline-based penalized estimator

In both the concurrent and the integral models presented in Section 2.1 and Section 2.2 respectively, we have used the decomposition of the infinite-dimensional functional covariates and parameters into a truncated functional basis depending on the chosen number of basis functions.

These values naturally needed to be correctly selected in order to avoid over-or under-fitting.

Nevertheless, precise adjustment of these values often induces a high computational effort. In the case of the B-spline basis, even more parameters have to be properly tuned such as the choice of the spline order and the location of the knots. In order to reduce the expected cost of such a computationally demanding procedure, we made the choice of choosing a sufficiently large a priori value for q β (or q γ ) and then apply a roughness penalty. This approach brings the benefit of reducing the overall computational cost, and of possibly improving the interpretability of the estimated functional coefficients. This last point is very interesting in the case of the linear model because as we already know, the interpretation of the predictors-response relationship becomes more difficult as the shape of the functional parameter β (or γ) does not have any simple structure.

Various approaches to regularize the parameter shape have been proposed in the literature.

In our setting of interest, the main idea is oftentimes to enhance the model performance and interpretability by adding a roughness penalty. [START_REF] Leurgans | Canonical correlation analysis when the data are curves[END_REF] is among the first to explore the functional penalization and show that the obtained estimator β(t) (resp. γ(s, t))

becomes less sensitive to the rather subjective choice of the number of basis functions q β (resp.

q γ ). More recently, [START_REF] James | Functional linear regression that's interpretable[END_REF] proposed a method called Functional Linear Regression

That is Interpretable (FLiRTI) which addresses the issue of choosing relevant penalties. Based on variable selection ideas such as the Lasso penalty, they produce accurate, flexible and highly interpretable estimates of the functional parameters. The main idea in James et al. ( 2009) is, instead of enforcing sparsity on the function themselves, to enforce sparsity of the derivatives instead. Using the notation β (l) (t) for the l th derivative of β(t), we may deduce that β (0) (t) = 0 guarantees X(t) has no effect on Y(t) at t; β (1) (t) = 0 implies that β(t) is constant at t; β (2) (t) = 0 means that β(t) is linear at t and so on. The FLiRTI approach also combine sparsity enforcing penalties for more than one derivative at a time, which can be useful for smooth parameters that may even vanish on some intervals.

Instead of the Lasso penalty applied in the FLiRTI method, where choosing the derivatives remains a difficult computational issue, our approach uses a Ridge penalty on the second derivative of the functional parameters. The choice of penalizing the second derivative is mainly motivated by the desire to obtain a possibly locally linear relationship if needed. Moreover, the use the Ridge penalty is motivated by the lack of exact sparsity observed in real problems and the clear benefits of getting a closed form formula for the estimators.

Penalized estimator for the concurrent model

Let us first consider the concurrent model in the classical mixed model form as in (7). In order to obtain an interpretable estimator, we will use a roughness penalty in the form of a Ridge-type penalty on the second derivatives, as advocated for in the previous paragraph. The objective function is the penalized log-likelihood function given by

L pen (b, Γ) = -2 L(b, Γ | Y) + p l=0 Pen(β l ) ; (12) with, L(b, Γ | Y) = nm log(2π) + log |V| + (Y -R ⊤ b) ⊤ V -1 (Y -R ⊤ b) (13)
using that V = Var(ZU + η) and the penalty

Pen(β l ) = λ l β ′′ l (t) 2 dt = λ l q β l j=1 b l j ϕ l " j (t) 2 dt = λ l q β l s,k=1 b l s b l k Φ l sk , with Φ l sk = ϕ l " s (t) ϕ l " k (t) dt.
When λ l is too large, the estimation of β l (.) will be too smooth, and we will not be able to account for the possible variations of the regression coefficients. When, instead, λ l is too small, the estimators might become too rough and overfitting might occur.

For a given value λ l , the estimation of β l (t) 0≤l≤p is obtained by solving:

min b, Γ L pen (b, Γ) = min b, Γ -2 L(b, Γ | Y) + p l=0 λ l q β l s,k=1 b l s b l k Φ l sk = min b, Γ -2 L(b, Γ | Y) + b ⊤ (λP) b , (14)
where λP ∈ R q β ×q β is given by:

λP =         λ 0 Ψ 0 0 q β 0 ×q β 1 . . . 0 q β 0 ×q β p 0 q β 1 ×q β 0 λ 1 Ψ 1 . . . 0 q β 1 ×q β p . . . . . . . . . . . . 0 q β p ×q β 0 0 q β p ×q β 1 . . . λ p Ψ p         with Ψ l =         Φ l 11 Φ l 12 . . . Φ l 1q β l Φ l 21 Φ l 22 . . . Φ l 2q β l . . . . . . . . . . . . Φ l q β l 1 Φ l q β l 2 . . . Φ l q β l q β l        
.

Here, 0 q 1 ×q 2 is the standard notation for the null matrix of size q 1 × q 2 . As Ψ l is a symmetric positive-definite matrix for any 0 ≤ l ≤ p, we can easily find its Cholesky decomposition, which can be efficiently leveraged in the implementation.

We first rewrite Model (7) in the form :

Y = R ⊤ b + ε * , with ε * = Z U + η and V = Var(ε * ) = Z Γ Z ⊤ + σ 2 I.
By setting the partial derivatives with respect to b and V to 0 and then solving the resulting linear system, we get:

b λ = R ⊤ V-1 R + λP -1 R ⊤ V-1 Y. ( 15
)
(see Appendix 7.2 for more details).

Let us now address the problem of choosing the smoothing parameters λ = (λ l ) 0≤l≤p . The correct choice will make great use of the observed accuracy of the prediction. For this purpose, for a fixed value of λ, we resort to a leave-one-out cross-validation type approach and compute b

(-i) λ
based on the sample except for the i th observation. We then compute the prediction Y (-i) λ at observation i. Finally, we can compute the prediction error or cross-validation score associated with the parameter λ as

V n (λ) = 1 n n i=1 Y i -Y (-i) λ 2 .
The value of λ that achieves the lowest estimated risk will be selected.

Penalized estimator for the integral model

For the integral model, we can also optimize the penalized log-likelihood function as in ( 12).

However, the main difference lies in the specific form of the penalty. The log-likelihood of the model will thus have the expression:

L(a, Γ | Y) = nm log(2π) + log |V| + (Y -Q ⊤ a) ⊤ V -1 (Y -Q ⊤ a) (16)
using V = Var(ZU + η), and the penalized log-likelihood:

L pen (a, Γ | Y) = -2 L(a, Γ | Y) + Pen(γ 0 ) + p l=1
Pen(γ l ), ;

For this model, the parameters γ l will be bivariate functions except for γ 0 , which is univariate.

The penalties for bivariate parameters will take the following expression:

Pen(γ l ) = λ l H γ l (t, s) 2 ds dt = λ l   ∂ 2 γ l (t,s) ∂t 2 ∂ 2 γ l (t,s) ∂t ∂s ∂ 2 γ l (t,s) ∂s ∂t ∂ 2 γ l (t,s) ∂s 2   2 ds dt.
Here H f (t, s) denotes the Hessian matrix of the bivariate function f and ∥ • ∥ is as is standard the Frobenius norm. To simplify expressions we will use the notation: ∂ 2 γ l (t, s) ∂t ∂s ≡ γ ts l (t, s), and then we have

Pen(γ l ) = λ l γ tt l (t, s) 2 + 2 γ ts l (t, s) 2 + γ ss l (t, s) 2 ds dt. ( 18
)
We know from (9) that

γ l (t, s) 2 = a l ⊤ B 1 (t)B 2 (s) 2 = q γ l i,j=1 a l ij B l 1i (t)B l 2j (s) 2 = q γ l i,j,k,m=1 a l ij a l km B l 1i (t)B l 1k (t) B l 2j (s)B l 2m (s),
so we then have the following expressions for the partial derivatives:

                                                         γ tt l (t, s) 2 ds dt = q γ l i,j,k,m=1 a l ij a l km B l ′′ 1i (t)B l ′′ 1k (t) dt B l 2j (s)B l 2m (s) ds = q γ l i,j,k,m=1 a l ij a l km Φ l ′′ 1,ik Φ l 2,jm ; γ ts l (t, s) 2 ds dt = q γ l i,j,k,m=1 a l ij a l km B l ′ 1i (t)B l ′ 1k (t) dt B l ′ 2j (s)B l ′ 2m (s) ds = q γ l i,j,k,m=1 a l ij a l km Φ l ′ 1,ik Φ l ′ 2,jm ; γ ss l (t, s) 2 ds dt = q γ l i,j,k,m=1 a l ij a l km B l 1i (t)B l 1k (t) dt B l ′′ 2j (s)B l ′′ 2m (s) ds = q γ l i,j,k,m=1 a l ij a l km Φ l 1,ik Φ l ′′ 2,jm . with the notation Φ l ′ u,sk = B l ′ us (t)B l ′ uk (t) dt.
Back to the problem of minimizing the penalized log-likelihood (17), for fixed values (λ l ) 0≤l≤p , the estimation of γ 0 (t), γ 1 (t, s), . . . , γ p (t, s) is obtained by solving the problem:

min a, Γ L pen (a, Γ) = min a, Γ -2 L(a, Γ | Y) + λ 0 q γ 0 s,k=1 a 0 s a 0 k Φ 0 sk + p l=1 λ l q γ l i,j,k,m=1 a l ij a l km Φ l ′′ 1,ik Φ l 2,jm + 2 Φ l ′ 1,ik Φ l ′ 2,jm + Φ l 1,ik Φ l ′′ 2,jm .
The penalty term for any bivariate parameter γ l (t, s) can be seen as the tensor product between the three following terms: the (a ij ) 1≤i,j≤q γ l matrix, the 4 th order square tensor of dimension q γ l and the matrix (a km ) 1≤k,m≤q γ l . We can rearrange this tensor product as a matrix product by flattening the matrix to a vector and the 4 th order tensor to a matrix. So we get a matrix product between the row vector of length q 2 γ l , the square matrix of dimension q 2 γ l × q 2 γ l and the column matrix of length q 2 γ l . The minimization problem can be written in matrix form:

min a, Γ L pen (a, Γ) = min a, Γ -2 L(a, Γ | Y) + a ⊤ (λP) a , (19)
where λP the matrix of dimension q γ × q γ with q γ = q γ 0 + p l=1 q 2 γ l defined as in ( 14). The main difference lies in the expression of the block matrix Ψ l for l > 0 given by:

Ψ l = Φ l ′′ 1,ik Φ l 2,jm + 2 Φ l ′ 1,ik Φ l ′ 2,jm + Φ l 1,ik Φ l ′′ 2,jm 1≤i,j,k,m≤q γ l
.

With this expression, we proceed in the same way as in the concurrent model to obtain the penalized estimator of a and then γl (t, s) l .

Prediction interval

We know, based on earlier works of [START_REF] Ruppert | Semiparametric Regression[END_REF] and Wood (2006), that the variance of our estimators and joint and point wise confidence intervals in the mixed effects model can easily be obtained. The penalized estimator of fixed effects and variance components is given by:

     b λ = (R ⊤ V-1 R + λP) -1 R ⊤ V-1 Y, U = σ 2 Z ⊤ V-1 (Y -R ⊤ b λ ) ;
where λ is the smoothing parameter. Since the above expressions depend on the variance components, they can only be calculated if they are known.

With a model of the form

Y = R ⊤ b + ε * with ε * = ZU + η, we have ε * ∼ N (0, V) with V = ZΓZ ⊤ + σ 2 I nm and then Y ∼ N (R ⊤ b, V). Thus: Cov(Y, U) = Cov(R ⊤ b + ZU + η, U) = Cov(R ⊤ b, U) + Z Var(U) + Cov(η, U) = Z Γ We can thus deduce that E(U | Y) = ΓZ ⊤ V -1 (Y -R ⊤ b), and since Y ∼ N (R ⊤ b, V) holds, Var( b λ ) = Var R ⊤ V-1 R + λP -1 R ⊤ V-1 Y = R ⊤ V-1 R + λP -1 R ⊤ V-1 Var(Y) V-1 R R ⊤ V-1 R + λP -1 = R ⊤ V-1 R + λP -1 R ⊤ V-1 V V-1 R R ⊤ V-1 R + λP -1 = R ⊤ V-1 R + λP -1 R ⊤ V-1 R R ⊤ V-1 R + λP -1
where one assumes that V is fixed and does not depend on Y. Therefore, the diagonal elements of this matrix are considered as estimates of Var( bλ,j ) j even though it is known to often underestimate it target. With these ideas in hand, we get bλ,j ± z 1-α/2 Var( b λ ) jj gives an approximate 100(1 -α)% confidence interval of bλ,j .

Based on this confidence interval of b we can easily build a pointwise confidence interval of the prediction Y i (.) at any desired time t.

Simulation study of functional models

The aim of this section is to illustrate and validate the estimation procedure described in Section 3 in the framework of "perfectly controlled" data, i.e. in the set-up where the assumptions about the distribution are the ones underlying our theory. The main properties of interest are accuracy, interpretability and smoothness of the estimated parameters, as well as the prediction quality.

Hence, we will first conduct a simulation study without penalization of the functional parameters in order to assess the relevance of the method based on prediction accuracy. Secondly, another simulation experiment compares unpenalized parameter estimation with its penalized counterpart.

Only the concurrent model is considered in this section, but similar results are expected for the integral model.

Data simulation process

The framework of the linear model we are considering in the present work can be helpful to explain the variations of a functional response variable through a set of controlled factors, which are the covariates or explanatory variables. In order to illustrate the relevance of our model, we are going to simulate an artificial data set and evaluate if the parameters are correctly recovered. For this purpose, we will first simulate the covariates and then use them as input to our regression model in order to simulate the corresponding response.

The p = 5 functional covariates are simulated at m = 50 equidistant viewpoints (t j ) j over the domain T = [0, 1] according the following procedure :

U l i (t j ) = ξ l i,1 + log(10 + t j ) ξ l i,2 + ξ l i,3 sin 2πt j ξ l i,4 (20)
where ξ l i,r is drawn from U [-1, 1] (1 ≤ r ≤ 4). This data, as we can see in Figures (1a)-(1e) inside Figure 1 for one randomly chosen individual, is generated at discrete timestamps over T = [0, 1] (blue dots). Before estimating the model, we first compute the underlying expansion into a basis of B-Splines.

We obtain the red curves in Figure 1 using K = 25 basis functions and equidistant distributed nodes. In other words, we observe U l i (t) with 1 ≤ l ≤ p and 1 ≤ i ≤ n, and the functional covariates X l i (t) are obtained as:

U l i (t) = X l i (t) + δ l i (t) with X l i (t) = K j=1
x l ij b j (t) and δ i (t) ∼ N (0, u 2 i ).

The functional parameters β(t) = β 0 (t), β 1 (t), . . . , β p (t) are chosen as follows: β 0 (t j ) = log(10 + t j ) ρ 0 and β l (t j ) = ρ l 1 sin 2πt j ρ l 2 with ρ 0 , ρ l 1 , ρ l 2 some constants given in Appendix 7.1. Figure 2 shows the corresponding representations of the functional parameters.

(a) β 0 (t) (b) β 1 (t) (c) β 2 (t) (d) β 3 (t) (e) β 4 (t) (f) β 5 (t)
Figure 2: Functional parameters in the concurrent model.

Given the proposed functional covariates and functional parameters, we can now compute the functional response using the concurrent Model (1), for two different sampling rates n ∈ {200; 500}.

In this experiment, ε i (t) is a Gaussian noise with mean 0 and two levels of variance σ 2 ∈ {1; 4}.

For each configuration, we run N = 50 Monte Carlo simulations.

Assessment criteria

We assess the performance of our estimation procedure. Two criteria are considered: prediction accuracy and estimation error for the model parameters. We extend to the functional framework the well-known Mean Relative Prediction Error (MRPE), which is used to quantify the distance between the actual and the predicted value of the functional response:

MRPE = 1 m m j=1 n i=1 Y i (t j ) -Y i (t j ) 2 n i=1 Y i (t j ) 2 . ( 21
)
We also define one extension of the determination coefficient, which consists of a simple arithmetic average of the classical determination coefficient along the time observation of the functional response. This determination coefficient noted R 2 is defined as follows:

R 2 = 1 m m j=1 1 - n i=1 Y i (t j ) -Y i (t j ) 2 n i=1 Y i (t j ) -Y i (t j ) 2 , (22)
where Y i (t j ) is the predicted output of the sample i at time t j and Y i (t j ) the mean function of the output sample at time t j .

To evaluate the performance of the estimation parameters, we compare the actual functional parameters with those provided by our models using the Mean Square Error (MSE) given by:

MSE β l (.) = p l=0 1 m m j=1 β l (t j ) -β l (t j ) 2 1/2 . ( 23
)

Simulation results

Figure 3 compares the boxplots of the Mean Square Error (23) of the estimated parameters.

As expected, the MSE decreases as the number of observations increases. This is the case in the scenarios when we have either a small or a high variance of the model error. Additional information about the results, the estimated functional coefficients versus the actual ones, are available in Appendix 7.3 in Figure 13. The functional determination coefficient computed over all the scenarios is presented in Table 1.

0.0 0.1 0.2 β 0 (t) β 1 (t) β 2 (t) β 3 (t) β 4 (t) β 5 (t) MSE Scenarios : n = 200, σ 2 = 4 n = 500, σ 2 = 4 n = 200, σ 2 = 1 n = 500, σ 2 = 1
From this table, we observe that when the additive noise increases, the coefficient of determination gets smaller, and increasing the sample size improves the coefficient of determination. Table 1: Functional determination coefficient R 2 over all the repetitions of any scenarios of simulation.

Scenarios

Let us now turn to our main objective, which is performance in prediction. We generate a test sample with n = 2000 observations, and we compare the difference between the actual values of the functional response and the prediction given for each model in a Monte Carlo simulation when n = 200 and n = 500. Accuracy is measured by the Mean Relative Prediction Error (MRPE). The boxplots for our four simlutation setups are given in Figure 4 while Figure 5 gives the actual values and the prediction over time. The simulations corroborate our expectations that our prediction scheme is able to cope with large variations in the functional response.

Effect of regularization

A close observation of the results shows that the basis expansion of β 0 (t) requires a large number of zero coefficients, which makes the estimation problem difficult to address in the small sample setting. As a notable consequence, estimation for this parameter may become unreliable for insufficiently large samples. Since the choice of the number of basis functions L β strongly affects the estimation of the functional parameters, our subsequent strategy will rely on introducing a regularization term. As explained in Section 3, regularization is a flexible and often robust way to adjust the number of basis functions. In order to illustrate the potential positive impact of using regularization, we propose a simulation study with p = 3 functional predictors whose parameters are constant (or linear) in some region and present high variability in other regions. We compare the unpenalized setting with the penalized setting when L β is chosen arbitrarily large. More precisely, we investigate the three following scenarios: L β = 50 without regularization, L β = 50 with regularization and L β = 5 without regularization.

The following functional parameters are considered:

β 0 (t) = 8t, β 2 (t) = β 1 (1 -t) and β 1 (t) =      0 if t ≤ 0.4 , 1.44 sin 2πt 0.38
otherwise.

;

β 3 (t) =      8t if t ≤ 0.4 , 3.2 + 1.44 sin 2πt 0.38 otherwise.
For the proposed model, L β = 50 basis functions is most certainly too large a number. Therefore, in addition to the prediction accuracy for these three models, we focus on the interpretability of the obtained functional parameters. The objective is to have a parameter that fits very well both in the regions where the function is constant as well as where the function undergoes high variability. L β = 50 with penalization in green which is completely hidden by the actual parameter in black.

Figure 6 plots one random instance, among the N = 50 simulations, of the functional parameters β 0 (t), β 1 (t), β 2 (t) and β 3 (t) in our three scenarios. Table 2 shows the average MSE (and standard deviation) between actual parameters and the estimated ones in the three considered scenarios. For the β 0 (t) parameter, which has a linear shape, the best fit comes from Scenario 1

where we have a small number of basis functions. The penalization process (Scenario 3) does not have the best performance on this parameter, but it has an acceptable shape as compared with the non-penalized process for the same number of basis functions (Scenario 2). The unpenalized estimator does not perform well especially at the start and at the end of the domain, while Scenario 1 does provide a sufficient number of basis functions to cope with the more complex behaviour in the very non-linear areas. Our main observation in this setup is that Scenario 3 is globally the best approach among the three since it is the only one that correctly adjusts its complexity to the oscillations of the function to estimate.

Mean Square Error β 0 (t) Table 2: Average (and standard errors) obtained over the N = 50 repetitions we performed of MSE between estimated parameters and actual ones in the three different scenarios.

β 1 (t) β 2 (t) β 3 (t) Scenario 1: L β =
An additional important observation is that choosing equispaced knots seems a sufficient strategy for most of the estimation problems we encountered. Sticking to this strategy allows avoiding the cumbersome task of selecting the locations of the knots using cross-validation.

Application to real data

In this section, we apply the proposed methodology for function-on-function regression (subsequently denoted by PenFFR which stands for "penalized function-on-function regression" and by FFR for "(unpenalized) function-on-function regression" ) for concurrent and integral models. These models are applied to two well-known data sets in FDA: Canadian Weather (CW) data available in the R package fda and Hawai Ocean (HO) data available in the R package FRegSigComp. We compare the prediction accuracy obtained using our method with the accuracy obtained with other existing methods: integral and concurrent Penalized Function-on-Function Regression (PFFR, [START_REF] Ivanescu | Penalized function-on-function regression[END_REF]) implemented in the R package refund; the signal compression approach (wSigcomp) designed by [START_REF] Luo | Functional wavelet regression for linear function-on-function models[END_REF] for the integral model and implemented in the R packages FRegSigComp; the Optimal Penalized Function-on-Function Regression (OPFFR) for the integral model (Sun et al., 2018), the Functional Principal Component Analysis (FPCA) and Functional Data Analysis method (FDA) [START_REF] Ramsay | Functional Data Analysis[END_REF]. Due to the unavailability of code for the OPFFR approach, we simply use the published results as presented in their paper (Sun et al., 2018).

Hyper-parameter tuning For our methods (FFR and PenFFR), we consider cubic B-splines basis functions for both functional predictors and regression coefficients. On the CW data set, we use 100 basis functions to address the functional and complex nature of the predictors and on HO data set, we use 40 basis functions. This choice is motivated by the fact that on the raw data, predictors on CW data set has 365 measurements while predictors in the HO data set have 200 measurements. The number of basis functions of parameters is set to 15 on CW data, both for integral and the concurrent models. For the HO data, based on the fact that we have 4 functional predictors and we know that the number of features of design matrix depends on the squared of the number of basis functions in the integral model. So for this complexity, we choose 40 basis functions for the concurrent model and only 6 for the integral model. The penalty parameters λ l of any predictor is selected using cross-validation on a predefined grid of values (10 equispaced values between 0.1 and 2.0).

For the PFFR method we used the default settings prescribed in the software and only set the number of basis functions for both the functional parameters and predictors. To correctly compare to our proposed method, we also used a cubic splines basis for both the functional predictors and parameters for the two (CW and HO) data sets. We use as our method the same number of basis functions to recover the functional nature of the predictors and on parameters.

For the wSigcomp method designed for the integral model, the default settings of the software are also used. For the HO data set which is tested by authors in their package description, the number of basis functions is set to 40 for the functional parameters and 20 for predictors. For the CW data, we slightly change but in the same proportion these value and set the number of basis functions involved for the functional parameters to 80 and the predictors to 40. We have detailed the choices of the hyperparameters but it should be noted that the performance of all these methods remains slightly sensitive to a reasonable variation of these values.

Canadian Weather Data

Hawaii ocean data Methods Type of basis

X ℓ i (t) β ℓ (t) Type of basis X ℓ i (t) β ℓ (t

Canadian weather data

The (2018) in terms of prediction accuracy over the 365 days of the year through the leave-one-out cross-validation integrated square error (ISE) given by:

ISE i = 365 0 Y i (t) -β (-i) X i (t) 2 dt
where the predictor X i (.) derives from the noisy daily temperature measurements; the functional response Y i (.) is the log daily precipitation and β (-i) is the functional parameter estimated in the data set of all the observations except for the i th observation.

For sake of reducing the computational burden, instead of the ISE, the L 2 -norm between the actual and prediction values on a grid of values t is used as a surrogate. It is given by:

ISE i = 365 j=1 Y i (j) -β (-i) X i (j) 2 . ( 24
)
The average ISE i values for the different models are given in Table 4. They show the numerical advantage of our proposed PenFFR method over the other methods. We also note that the variance observed in our predictions remains quite high for the different models. This is due to the quality of the input data. Also shown in Figure 8 is the prediction obtained using the different methods. We restrict our attention to the integral model since it appeared to be the best model for this data set, independent of the estimation method (PenFFR, PFFR and wSigcomp). The prediction is given for two randomly chosen weather stations (Iqaluit and Arvida) and are compared with the actual precipitation. Similar results are illustrated by Figure 14 in the appendix for the concurrent model.

- actual data, the red two-dashed line represents the prediction given by our integral PenFFR, the magenta dashed line is the prediction given by the integral PFFR method, and the long-dashed green line is the prediction given by the wsigcomp method.

Hawaii ocean data

This data set is one of those used by [START_REF] Luo | Functional wavelet regression for linear function-on-function models[END_REF] to apply their wSigcomp approach. The as functional predictors X(t) = X 1 (t), X 2 (t), X 3 (t), X 4 (t) . We split the full data set into two train/test sub-data sets where the training data consists of the 50 first days (observations) only.

First of all, we expand all the functions considered into a cubic B-spline basis with 40 basis functions. Figure 9 displays the sample curves for these variables. Our PenFFR method is compared with PFFR and wSigcomp in the setting of integral models.

We also consider PenFFR and PFFR for the concurrent model. Figure 10 and 11 show the estimated parameters γ0 (t) and γj (t, s), 1 ≤ j ≤ 4 obtained for the three methods in the case of the integral model. We first notice that the shape of the estimated parameters is smooth for our method (third column). In addition, Figure 15 in the appendix shows the estimates βj (t), 0 ≤ j ≤ 4 of the concurrent model with the PenFFR and PFFR methods. Prediction accuracy using ISE on a test set of size 66 is shown in Table 5. Since the number of individuals for this data ( 116) is larger than the size of the previous data set, we evaluate the performance on a single test set rather than using cross-validation in order to circumvent the potentially heavy computational burden. Our method is seen once again to outperform all other methods as illustrated in Figure 12 which shows predictions on two randomly chosen individuals.

Methods

ISE (×10

Conclusion

In this article, we have presented a new estimation process for the linear regression model with functional responses and functional covariates. We approach the problem via expanding the functions onto a common B-spline basis, hence allowing the reduction of the functional model to a linear mixed model. Adaptation to unknown smoothness is performed by adding a roughness penalty on second derivatives. Unlike any estimator based on basis functions, our estimates have a smooth shape and sufficient flexibility to capture the encountered variability in various experiments with real-world data sets. We then illustrate the performance of our proposed estimation process in terms of prediction accuracy and parameter interpretability on simulated and real data sets.

Perspectives for future work on this model are manifold. First, prediction confidence bounds can be obtained using various methods such as conformal prediction [START_REF] Angelopoulos | A gentle introduction to conformal prediction and distribution-free uncertainty quantification[END_REF], which can handle black box models and could be adapted to our setting as well. Another avenue for future investigations is to explore mixture function-on-function models. This type of mixture model can be safely expected to be extremely relevant when heterogeneous clusters are present in the population [START_REF] Desarbo | A maximum likelihood methodology for clusterwise linear regression[END_REF]). Mixture of experts can also be explored as an additional extension which could prove very efficient in predictive modelling; see [START_REF] Chamroukhi | Functional mixturesof-experts[END_REF], where a new family of FME is proposed, albeit restricted to scalar responses.

Appendix

Simulation parameters

The values of the chosen constants is drawn from uniform law between -5 and 5. The values is given by: ρ 0 = 0.439, ρ 1 1 = -3.562, ρ 2 1 = -1.058, ρ 3 1 = -2.955, ρ 4 1 = -0.585, ρ 5 1 = -0.298, ρ 1 2 = 0.228, ρ 2 2 = 2.641, ρ 3 2 = 4.462, ρ 4 2 = 2.757 and ρ 5 2 = 2.283.

Mixed model estimator

We first rewrite the model in the form :

Y = R ⊤ b + ε * , (25)
with ε * = Z U + η, from which we get V = Var(ε * ) = Z Γ Z ⊤ + σ 2 I. We aim to estimate the fixed effects b and the error variance V from the observed data. The most popular estimation methods for the parameters in Model ( 7) are maximum likelihood (ML) and restricted maximum likelihood (ReML) as described in [START_REF] Lindstrom | Newton-raphson and em algorithms for linear mixedeffects models for repeated-measures data[END_REF]. The log-likelihood of the model is written

∂L pen ∂b = ∂ ∂b Y ⊤ -(R ⊤ b) ⊤ V -1 Y -(R ⊤ b) + b ⊤ (λP) b = ∂ ∂b Y ⊤ V -1 Y -Y ⊤ V -1 R ⊤ b -(R ⊤ b) ⊤ V -1 Y + (R ⊤ b) ⊤ V -1 R ⊤ b + b ⊤ (λP) b = -(Y ⊤ V -1 R ⊤ ) ⊤ -R V -1 Y + 2 R V -1 R ⊤ b + 2 (λP) b = -2 R V -1 Y + 2 (R V -1 R ⊤ + λP) b.
and by equalizing to 0, i.e. ∂L pen ∂b = 0, we get:

b(V) = R V -1 R ⊤ + λP -1 R V -1 Y. ( 27
)
By replacing b by its estimator in the likelihood expression, we get the profiled log-likelihood given by:

L p (V) = - 1 2 N log(2π) + log |V| + Y -R ⊤ R V -1 R ⊤ -1 R V -1 Y ⊤ V -1 Y -R ⊤ R V -1 R ⊤ -1 R V -1 Y = - 1 2 N log(2π) + log |V| + Y ⊤ V -1 -Y ⊤ V -1 R ⊤ R V -1 R ⊤ -1 R V -1 Y -R ⊤ R V -1 R ⊤ -1 R V -1 Y = - 1 2 N log(2π) + log |V| + Y ⊤ V -1 Y -Y ⊤ V -1 R ⊤ R V -1 R ⊤ -1 R V -1 Y - Y ⊤ V -1 R ⊤ R V -1 R ⊤ -1 R V -1 Y + Y ⊤ V -1 R ⊤ R V -1 R ⊤ -1 R V -1 R ⊤ R V -1 R ⊤ -1 R V -1 Y = - 1 2 N log(2π) + log |V| + Y ⊤ V -1 Y -Y ⊤ V -1 R ⊤ R V -1 R ⊤ -1 R V -1 Y L p (V) = - 1 2 N log(2π) + log |V| + Y ⊤ V -1 I -R ⊤ R V -1 R ⊤ -1 R V -1 Y .
On the other hand, there holds V = Var(ε * ) = σ 2 U Z Z ⊤ + σ 2 I , and thus, L p (V) = L p (σ 2 U , σ 2 ). It is obviously not easy to derive this likelihood which no longer depends on b. Moreover, maximizing this last function gives the MLE which is nevertheless biased. For these reasons, and in order to account for the degrees of freedom of the fixed effects in the model, we propose to use the Restricted Maximum Likelihood (ReML) which reads:

L R (V) = L p (V) - 1 2 log |R V -1 R ⊤ | (28)
From a numerical viewpoint, we obtain the estimator V of the variance V by maximizing this last likelihood from which we finally deduce the value of Û given by:

     b = (R ⊤ V-1 R + λP) -1 R ⊤ V-1 Y, U = σ 2 Z ⊤ V-1 (Y -R ⊤ b). (29

Figure 1 :

 1 Figure 1: Simulated functional predictors and functional (concurrent) response for a randomly chosen individual.

Figure 3 :

 3 Figure 3: Boxplots of Mean Square Error of estimated parameters over the N = 50 Monte Carlo simulation when n = 200 (red) and n = 500 (blue).

R 2 n

 2 = 200, σ 2 = 4 0.868 (0.0065) n = 500, σ 2 = 4 0.878 (0.0040) n = 200, σ 2 = 1 0.946 (0.0022) n = 500, σ 2 = 1 0.947 (0.0013)

nMRPEFigure 4 :Figure 5 :

 45 Figure 4: MRPE on a test sample of length n = 2000 in all the scenarios of simulation for Monte Carlo simulation

Figure 6 :

 6 Figure 6: Estimated parameters vs actual ones in the 3 scenarios for any parameters in one of the N = 50 simulations. L β = 5 without penalization in cyan ; L β = 50 without penalization in red ;

  data set consists of m = 365 daily temperature measurements (average over the years 1961 to 1994) at n = 35 weather stations in Canada and their corresponding daily precipitation (in log scale). The weather stations are located in K = 4 climate zones: Atlantic, Pacific, Continental and Arctic and the aim is to use the daily temperature to predict the precipitation at each station.

Figure 7 Figure 7 :

 77 Figure7gives the daily average over the years 1961 to 1994 (temperature on the left, precipitation on the right). Note that the stations in the Pacific zone have the highest precipitation values, and stations from this zone also have the highest temperatures in the winter. The same can be said about the stations in the Arctic zone for low temperatures and precipitation. A positive relationship between temperature and precipitation can therefore be suspected.

Figure 8 :

 8 Figure 8: Prediction on two randomly chosen stations. For each figure, the black points are the

  data set includes physical and biochemical oceanographic observational data from the Hawaii Ocean Time-series (HOT) Program, including thermosalinograph, Conductivity, Temperature and Depth (CTD), bottle and biochemical data. The HOT program makes repeated observations of the physics, biology and chemistry at a site approximately 100 km north of Oahu, Hawaii. In the data set, five variables: Salinity, Potential Density, Temperature, Oxygen and Chloropigment, are observed every two meters between 0 and 200 meters below the sea surface on 116 different days. This data set is available from the R package "FRegSigComp", under the name Ocean data. It consists of 5 functional variables with 116 individuals, each having 101 measurement points. Here, we consider the function-on-function regression model with the salinity curves as the response variable Y(t) and (Potential Density, Temperature, Oxygen, Chloropigment) curves

Figure 9 :

 9 Figure 9: Original sample curves of predictors expanded by cubic B-splines basis with 40 basis functions.

Figure 10 :

 10 Figure 10: Estimates γj (s, t), 1 ≤ j ≤ 4 for the three methods: wSigcomp (left column), integral PFFR (middle column) and integral PenFFR (right column).

Figure 11 :

 11 Figure 11: Estimates γ0 (t), 1 ≤ j ≤ 4 for the three methods: wSigcomp (left column), integral PFFR (middle column) and integral PenFFR (right column).

Figure 12 :

 12 Figure 12: Prediction given by the three methods for integral model on two randomly chosen observations in the test sample: PenFFR in blue, PFFR in green and wSigcomp in red. Black dots are the true values.

L

  pen (b, V) = nm log(2π) + log |V| + (Y -R ⊤ b) ⊤ V -1 (Y -R ⊤ b) + b ⊤ (λP) b (26)First order condition: ∂ ∂b L pen (b, V) = 0.

Figure 15 :

 15 Figure 15: Estimates βj (t), 0 ≤ j ≤ 4 for the two methods (pffr and PenFFR) on concurrent model. The first row shows the estimation provided by our PenFFR method. The second row shows the estimation provided by the pffr method.

  5 without penalization 0.20 (0.04) 0.77 (0.001) 0.77 (0.000) 0.82 (0.001)

	Scenario 2: L β = 50 without penalization	166.58 (9.1)	4.09 (3.0)	4.46 (3.4)	3.61 (2.77)
	Scenario 3: L β = 50 with penalization	0.38 (0.002) 0.06 (0.002) 0.06 (0.002) 0.06 (0.002)

Table 3 :

 3 Number of basis functions for the regression coefficients β ℓ (t) and the covariates X ℓ i (t)

	)

Table 4 :

 4 The average (and standard deviation) of ISE for the Canadian Weather data set. The best result is in boldface.

	For recall that we are trying to predict precipitation from temperature on a
	dataset of 35 very different weather stations

Table 5 :

 5 The average (and standard deviation) of ISE for the Hawaii ocean data set. The best result is in boldface.

		35.3
							35.3
		35.2					35.2
	Salinity	35.1				Salinity	35.1
		35.0					35.0
		34.9					34.9
		0	50	100	150	200	0	50	100	150	200
				meter below the sea surface					meter below the sea surface

models. Electronic Journal of Statistics, 10(1):495 -526. Silverman, B. W. (1996). Smoothed functional principal components analysis by choice of norm.

The Annals of Statistics, 24(1):1 -24. Sun, X., Du, P., Wang, X., and Ma, P. (2018). Optimal penalized function-on-function regression under a reproducing kernel hilbert space framework. Journal of the American Statistical Association, 113(524):1601-1611. Wood, S. N. (2006). On confidence intervals for generalized additive models based on penalized regression splines. Australian & New Zealand Journal of Statistics, 48(4):445-464.

Parameter representation on simulated data

In each scenario, we estimate the functional parameters with cubic B-splines basis, regular knots over the grid and L β l = 50 basis functions. The parameters we obtain with our model are close to the true parameters. However, we note that estimation of β 0 (t) is noised by the two large number of basis functions considered. This confirms the previously mentioned concerns about interpretability (smoothness) of the estimated parameters without regularization. Figure 13 also confirms that estimation accuracy increases with the number of observations.