Gnanguenon Girault
email: girault.gnanguenon@gmail.com

Patrice Guesse

Bénedicte Loisel

Thierry Fontez

Nadine Simonneau

Hilgert

Girault Gnanguenon Guesse

Patrice Loisel
email: patrice.loisel@inrae.fr

Bénedicte Fontez
email: benedicte.fontez@supagro.fr

Thierry Simonneau
email: thierry.simonneau@inrae.fr

Nadine Hilgert
email: nadine.hilgert@inrae.fr

An exploratory penalized regression to identify

Keywords: joint distribution, penalized linear regression, information criteria, Generalized Fused Lasso, interpretable coefficient

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Nowadays, several fields of activity and in particular, agriculture, are being revolutionized by the emergence of sensor data. With regard to crops, the setting up of harvest can now be monitored with the aim of including/modeling the influence of multiple environmental conditions. Specifically, water scarcity and temperature increase are two major features which have long been analyzed as determining huge variation in crop yield. Their influences are increasing with climate change and are becoming a major concern for the sustainability of agriculture in many parts of the world. However, relationships between climatic conditions and quality of the harvest are still poorly understood and modeling approaches are still lacking. To better use newly available data from sensors, there is a need for methods able to explore which combination of climatic variables influences harvest quality and at which stage of plant development. Such data sets variate sample having as dimension the number of discretization points of the functional variables.

In this case, multivariate statistics meets limits, among which its failure to take into account the very strong colinearity existing between discretized variables. By contrast, colinearity can be considered in the Fused Lasso penalized regression which was therefore retained in our approach. We followed a "scalar on image" regression model, where the "image" is a contingency table of the joint class intervals, for a fixed candidate partition, to which we associated a graph of contiguity constraints defined in subsection §2. [START_REF] Arnold | genlasso: Path algorithm for generalized lasso problems[END_REF]. Identification was performed through a Generalized Fused Lasso (see §2.4) using each candidate contingency table as input variables. The selection of the best candidate and of its relative regression coefficients was achieved by minimizing an information criteria. To combine our exploratory objectives with technical constraints (such as potentially a small amount of data, for example), we have favored an iterative approach (see §2.5) which explores a larger space of solutions, while keeping aware of possible slight overestimation.

Transformation of both functional variables

Let us consider the samples (A i ) i=1,...,n and (B i ) i=1,...,n of two explanatory functional variables A and B, associated to a scalar response variable y with samples (y i ) i=1,...,n , where n is the number of statistical individuals. Both A and B are observed on the same set T containing equidistant observation times. It is assumed that they have no missing values. These practical conditions of use can be released with some pre-treatment of the data (like interpolation, smoothing and imputation), see Section 6 for more details.

The requirement for variable transformation is intrinsically linked to the goal of the approach: identify joint class intervals of the explanatory variables that influence the response. The transformation requires the definition of joint class intervals which can be used as linear regressors to predict the response. The steps to achieve this transformation are shown in the Figure 2.1.

Contingency table of the joint class intervals

Let's partition the first explanatory variable A in n A class intervals according to a linear scale. This partition generates n A + 1 breaks denoted L A (v), v = 1, . . . , n A + 1. We chose to have equidistant breaks, as defined in equation (2.1) : Let's define a partition vector u = (n A , n B ). For all i, it is then possible to obtain the frequency bivariate histogram of (A i , B i ) as a contingency table C u i , of dimension n A ×n B , whose components

L A (v) = A + v -1 n A Ā -A , v = 1, . . . , n A + 1 (2.
C u i,(v,w) are obtained through (2.
2). The new variables are the joint class intervals, noted

I u A (v) × I u B (w), v = 1, . . . , n A , w = 1, . . . , n B .
They are called joint modalities.

C u i,(v,w) = t∈T 1 A i (t)∈I u A (v), B i (t)∈I u B (w) = Card {t ∈ T |A i (t) ∈ I u A (v), B i (t) ∈ I u B (w)}, (2.2)
for all v = 1, . . . , n A , w = 1, . . . , n B and each u = (n A , n B ), with:

n A v=1 n B w=1 C u i,(v,w) = Card(T ).
C u i,(v,w) is the number of times that the observations of A i and B i are at the same time in

I u A (v) × I u B (w). C u i,(v,
(v) × I u B (w)) v,w is a collection of 2D
intervals in which the pairs (A i , B i ) will be projected. This is detailed in the following subsection.

Functional model

The functional model used in the present work can be formulated as follows:

y = F u,β (A(t), B(t))dt + ε (2.3)
where F u,β = I u .β (scalar product of two vectors), I u is a vector of components functions I u v,w defined by

I u v,w (a, b) = 1 a∈I u A (v), b∈I u B (w) for a ∈ [A, Ā] and b ∈ [B, B]
, and β = (β (v,w) ) v,w is the vector of unknown coefficients to estimate. ε is an i.i.d. Gaussian error. Through (2.3), we note that the relationship between the variable to be explained y and the predictors A(.) and B(.) is stable over time: the function F u,β (and thus β) is not time-dependent. It is hence assumed that certain ranges of cross-values have an influence and that only the time spent in these ranges is important.

A discrete version of the model (2.3) can be written as follows (up to a constant), for each i :

y i = t∈T F u,β (A i (t), B i (t)) + ε i = t∈T I u (A i (t), B i (t)).β + ε i , which, combined with (2.2), yields C u i,(v,w) = t∈T I u v,w (A i (t), B i (t)
) and:

y i = n A v=1 n B w=1 C u i,(v,w) β u (v,w) + ε i (2.4)
This model is a linear multiple regression model where the regressors are the frequencies associated to the joint class intervals

(I u A (v) × I u B (w)) v,w . The indicator functions (I u v,w ) v,w define a functional basis.
From the contingency tables C u i , we construct the design matrix X u associated to model (2.4) as follows. After vectorization (stacking column by column) and transposition of the contingency table C u i (see part 2 of Figure 2.1), we obtain, for a fixed partition vector u = (n A , n B ), a row vector of length n A .n B :

X u i = t V ect(C u i ) , X u i ∈ R n A n B (2.5)
which represents the number of time observations t during which an individual i has been observed in each of the n A × n B levels described by the joint class intervals. The n stacked row vectors form the matrix

X u = t (X u 1 , X u 2 , . . . , X u n ) ∈ R n×n A n B .

Creation of a graph of contiguity constraints

To each matrix X u corresponds a graph G u (V u , E u ), which contains the contiguity constraints between modalities of the contingency table. V u represents the columns (new variables) of the candidate matrix X u and E u all the edges connecting two close joint modalities. We used the Rook's case contiguity rule [START_REF] Plant | Spatial Data Analysis in Ecology and Agriculture Using R[END_REF] where two joint modalities are said to be close if the bins following the variable A (indexed by v) or (exclusive) the bins following the variable B (indexed by w) are consecutive.

Selection of class intervals and related regression coefficients

The Fused Lasso is a variant of the Lasso introduced in 2005 by [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF], in order to take into account the existence of a structure in the variables. In its original form, the Fused Lasso aims not only at parsimony of coefficients but also at parsimony of differences in consecutive coefficients.

This version of the Fused Lasso can be interpreted as a one-dimensional Fused Lasso (1D-Fused Lasso). The Generalized Fused Lasso (GFL) [START_REF] Xin | Efficient generalized fused lasso and its application to the diagnosis of alzheimer's disease[END_REF] aims to promote smoothness over neighboring variables on a general graph G = (V, E) made of V knots and E edges. Each explanatory variable corresponds to a node on the graph and an edge symbolizes the link between a pair of separate nodes in G.

For a fixed partition vector u = (n A , n B ), the GFL criterion to minimize is written:

1 2 n i=1 (y i -X u i β) 2 + λ p j∈V u |β j | + λ f (j,j ) ∈ E u |β j -β j | (2.6)
with respect to β, where:

• β = t (β (1,1) , β (2,1) , . . . , β (1,2) , . . . , β (n A ,n B ) ) ∈ R n A n B the unknown coefficients,
• λ p ≥ 0 and λ f > 0 the regularization parameters (of parsimony and fusion) to be optimized,

• for j = (v, w) fixed, the couples (j, j ) relative to j and contained in E u are (j, j ) 1 = ((v, w), (v + 1, w)) and (j, j ) 2 = ((v, w), (v, w + 1)). In the following and depending on the context, the index j will refer either to the pair (v, w), or to the j th element of the vector obtained from the matrix stored by columns.

The argmin solution of (2.6), denoted β u , is computed as a function of the regularization parameters λ p and λ f , for a fixed value of u.

2.4.1. The Generalized Fused Lasso in the Generalized Lasso framework

If differences of contiguous coefficients were not penalized in (2.6) (i.e. if λ f was zero), then the criterion would reduce to the Lasso criterion presented by [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. The Lasso minimizes the residual sum of squares subject to the the constraint that the sum of the absolute value of the coefficients is less than a constant. In this case, there is only one regularization parameter to estimate. [START_REF] Efron | Least angle regression[END_REF] proposed the Least-Angle Regression (LARS) algorithm able to solve the problem for all λ ∈ [0, ∞[, producing a full piece-wise linear solution path. The result of the path algorithm is the finite set of increasing λ values, where each λ delimits a model reduction dimension (the number of non-zero β components). Criterion (2.6) has two regularization parameters λ p and λ f to be optimized. The path algorithm is no longer suitable to identify them.

Our proposal is to parameterize (2.6) with the ratio:

γ = λ p λ f (2.7)
(choice made in the R package genlasso [START_REF] Arnold | genlasso: Path algorithm for generalized lasso problems[END_REF] used for implementing SPICEFP). This ratio represents a balance between parsimony and fusion. Then, for a fixed value of γ, criterion (2.6) can be equivalently rewritten as:

1 2 ||y -X u β|| 2 2 + λ||D u,γ β|| 1 , (2.8)
where y = t (y 1 , y 2 , . . . , y n ) ∈ R n is the response vector and D u,γ is a specified penalty matrix (see below). This model corresponds to the Generalized Lasso model, introduced by [START_REF] Tibshirani | The solution path of the generalized lasso[END_REF] as an encapsulation of statistical models using the L 1 norm to impose additional constraints. Through this new parametrization, the value of λ = λ f can be optimized with the path algorithm and a corresponding β u,γ (λ) can be estimated. Finally, there are as many pairs (λ, β u,γ (λ)) solutions as there are u and γ parameters set. So, several models are available and we have to select one of them in order to deduce the optimal pair. Selection of the best model is done with an information criteria which requires an estimation of the degree of freedom for each model. Parameter λ p will be deduced from γ through λ p = γλ f .

In our context, the penalty matrix D u,γ is a row-binding of two sub-matrices, namely:

• D u,γ,p ∈ R n A n B ×n A n B : the penalty sub-matrix associated to the regularization of parsimony

λ p j∈V u |β j | • D u,f ∈ R 2n A n B ×n
A n B : the sub-matrix associated to the regularization of the fusion according to the two dimensions λ f

(j,j ) ∈ E u |β j -β j | . Affecting two dimensions, D u,f can be subdivided into D u,f 1 and D u,f 2 . Hence D u,γ =   D u,f 1 D u,f 2 D u,γ,p   ∈ R 3n A n B ×n A n B with: D u,f 1 (v,w)(v ,w ) =            1 if (v , w ) = (v + 1, w) -1 if (v , w ) = (v, w) and v < n A 0 if not D u,f 2 (v,w)(v ,w ) =            1 if (v , w ) = (v, w + 1) -1 if (v , w ) = (v, w) and w < n B 0 if not D u,γ,p = γ.I n A .n B (2.9)
where γ ≥ 0 and I n A n B is an identity matrix.

The Generalized Fused Lasso presented in criterion (2.6) is a 2 Dimensional -Sparse Fused Lasso (2d-SFL).

Degrees of freedom of the Generalized Fused Lasso fit

Our approach follows that of [START_REF] Tibshirani | Degrees of freedom in lasso problems[END_REF], who established the calculation of the degree of freedom for any Lasso-type regression written as a generalized Lasso problem as presented in (2.8).

We first introduce some notations. For any penalty matrix D ∈ R m×p involved in a generalized Lasso problem of type (2.8), let S be the active set corresponding to a particular solution β, it is defined as:

S = {r ∈ {1, . . . , m} : (D β) r = 0} = support(D β).
Let D -S be the matrix D from which were removed the rows indexed by S. And let null(D -S ) be the null space or kernel of D -S .

The calculation of the degrees of freedom of the generalized Lasso fit is stated in [START_REF] Tibshirani | Degrees of freedom in lasso problems[END_REF]Theorem 3] and reminded here :

Theorem (Generalized lasso degrees of freedom -Tibshirani and Taylor (2012))

Assume that y ∈ R n follows a normal distribution (y ∼ N (µ, σ 2 I) with given (unknown) mean vector µ ∈ R n and marginal variance σ 2 ). For any fixed and nonrandom predictor matrix X ∈ R n×p , penalty matrix D ∈ R m×p and λ ≥ 0, the degree of freedom of the generalized Lasso fit can be expressed as

df (X β) = E[dim(X(null(D -S )))] (2.10)
with S = S(y) the active set corresponding to any generalized Lasso solution β(y) at y.

The notation X(V ) represents the image space of a subspace V by X. It is the space generated by the columns of the X matrix projected on V . The assumptions required for this Theorem are those usually made in regression estimations (Gaussian i.i.d. errors). No assumptions are made on matrix D nor on X. This result can thus be applied to the 2d-Sparse Fused Lasso constraint.

To compute df (X u βu,γ ), we need β u,γ , X u , D u,γ , λ defined in (2.8) and γ in (2.7). Note that, in our context, m = 3n A n B the number of rows in D u,γ and p = n A n B the number of columns in X u .

With fused lasso constraints, we are interested in sets of parameters that share the same value.

By this, we refer to connected components, which are illustrated on Figure 2.4.2 and defined as follows:

Definition 2.1. A connected component cc is a set of indexes of non zero coefficients β u,γ v,w that are linked together via the D u,γ matrix (2.9) and that all share the same real value. Let us consider that the estimated coefficient β u,γ contains Q u,γ connected components denoted cc q , q = 1, 2, . . . , Q u,γ . In order to identify the respective coefficients involved in each non-zero connected component, we introduce the matrix Θ through Definition 2.2.

Definition 2.2. Let (cc q ) q=1,...,Q be a set of connected components. The connected component membership matrix is a binary matrix Θ ∈ R n A n B ×Q whose q-th column Θ (q) indicates the membership or not of the βu,γ components to the connected component cc q , as follows:

∀j = 1, . . . , n A n B , Θ (q) j = 0 if j / ∈ cc q 1 if j ∈ cc q , q = 1, . . . , Q. (2.11)
Let's now compute the null space of D -S , denoted null(D -S ) :

null(D -S ) = {Θ ∈ R n A n B ×Q |D -S Θ (q) = 0 (m-s) , q = 1, . . . , Q} (2.12)
where D -S ∈ R (m-s)×n A n B , D -S Θ (q) ∈ R m-s and s = Card(S).

As the D u,γ matrix has a simple structure, adapting Theorem 3 of [START_REF] Tibshirani | Degrees of freedom in lasso problems[END_REF] in the context of 2d-Sparse Fused Lasso is equivalent to looking for the components of β u,γ which have different values. This is the subject of the next corollary.

Corollary 2.3. The degree of freedom df (X β u,γ ) associated to the criterion (2.8) in the context of 2d-Sparse Fused Lasso is equal to the number of connected components Q u,γ .

Proof:

df (X β u,γ ) = dim X u (null(D u,γ -S )) = dim (X u (Θ)) = dim V ect X u Θ (q) , q = 1, . . . , Q u,γ = rank([X u Θ (1) , . . . , X u Θ (Q u,γ ) ]) = Q u,γ
where we omitted Θ dependencies on u and γ to lighten the notations in the corollary.

Choice of the best candidate matrix and selection of its variables

SPICEFP requires the construction of different candidate explanatory matrices X u from both functional variables and partition's vector u. Constructing a GFL for a matrix of predictors associated to a fixed u requires identifying the optimal values of the penalty parameters: λ and γ in (2.8). In penalized regressions, cross-validation is often used to optimize regularization parameters, but it is time consuming. We suggest using an information criterion to achieve the same purpose [START_REF] Garcia | Real-time inflation forecasting with high-dimensional models: The case of brazil[END_REF]. To that aim, we computed an adapted information criterion for each model indexed by u, λ and γ. The best model is obtained by minimizing the information criterion chosen, which yields the best partition u and allows to select the best variables from X u . The variable selection is done through the selection of γ and λ and the associated non-zero β components are deduced.

There exist various information criteria including Akaike Information Criterion (AIC) [START_REF] Akaike | Information Theory and an Extension of the Maximum Likelihood Principle[END_REF] and Bayesian Information Criterion (BIC) [START_REF] Schwarz | Estimating the dimension of a model[END_REF]. These criteria penalize log-likelihood by the number of model parameters. The BIC also penalizes log-likelihood by the sample size. A well defined penalization is essential to compare regression models involving different explanatory matrices. These information criteria require computing the degree of freedom Q u,γ of the GFL model, obtained in Corollary 2.3.

For each u and γ defined on respective grids (given by the users) and λ taken from the set (λ e ) e=1,...,N λ delivered by the path algorithm, we considered the following information criteria:

• AIC u,γ e = -2 log(L u,γ e ) + 2Q u,γ e • BIC u,γ e = -2 log(L u,γ e ) + log(n)Q u,γ e with L u,γ
e the likelihood function of the following model: y = X u β u,γ e + ε with ε ∼ N (0, σ 2 I), associated with criterion (2.8). We have:

-2 log(L u,γ e ) = 2n log(σ) + n log(2π) + 1 σ 2 ||y -X u β u,γ e || 2 2
where the variance of the residuals σ 2 is unknown. As mentioned by [START_REF] Hirose | Tuning parameter selection in sparse regression modeling[END_REF], the same variance estimator must be used in the calculation of the criteria for all the constructed models, and at each iteration (see bellow). We thus decided to estimate σ 2 by the variance of the response variable:

σ 2 = 1 n-1 ||y -y|| 2 2 .
It's a biased estimator of σ 2 , but this bias remains fixed for all models compared [START_REF] Hirose | Tuning parameter selection in sparse regression modeling[END_REF]. Such an estimator leads to an overestimation of the variance, which penalizes the introduction of new coefficients in the model. This bias can be partly offset by an iterative approach. That's why we have implemented an iterative algorithm where residuals are used as new response variable and so on, until shutdown conditions are verified.

SPICEFP: an iterative approach

To capture all potential non-zero coefficients, one should take a fine partition with values of n A and n B large enough. This could imply a very low or even zero number of points in the joint class intervals (making the method ineffective) and prohibitively long computation times [START_REF] Mairal | Complexity analysis of the lasso regularization path[END_REF]. As a trade off between thinness and work-ability we chose to develop an iterative approach to explore a large space of solutions (that allows addition of different thinness of partition).

Data required for the use of SPICEFP are the functional explanatory variables A i and B i discretized on a grid T and a response variable y i with i = 1, . . . , n. Other elements are also required at the input of the algorithm: Γ, a set of positive reals representing γ ratios of regularization parameters, U A and U B the sets of numbers of class intervals n A and n B , n λ the selected number of pairs (among N λ ) (λ, β u,γ (λ)), the information criterion to be used, and K the maximum number of iterations to explore. The n λ values of λ are chosen equally spaced on the log scale (see Genlasso package [START_REF] Arnold | genlasso: Path algorithm for generalized lasso problems[END_REF]).

SPICEFP constructs for each couple (u, γ), a matrix of explanatory variables X u and a penalty matrix D u,γ . This first step (line 1 to 13, algorithm 1) is performed only once and the set of candidate explanatory matrices remains unchanged. The second step (line 14 to 34, algorithm 1) of the approach is iterative. For each candidate matrix a solution path (λ, β u,γ (λ)) is obtained.

Once the criterion associated with each model is computed, the optimal triplet ( u, γ, e) is the argument that minimizes Crit u,γ e , where e is the index of the solution path-coefficients (λ e , βu,γ e ).

X u , D u, γ and β u, γ e respectively represent the optimal matrices of the explanatory, penalty and coefficient variables. At each iteration k, we denote u k = u and β k = β u, γ e . SPICEFP then checks if the selected coefficients according to the criterion Crit u,γ e correspond to a zero vector or if the maximum number of iterations K is reached. The algorithm is stopped when at least one of these conditions is verified. When none of these conditions are verified, the residuals of the optimal model are computed and used as a response variable at the next iteration of the algorithm.

The third step (line 35, algorithm 1) gives the final prediction as the sum of all the predictions obtained at each iteration:

k * k=1 X u k β k .
The function corresponding to the final prediction can be written as follows :

F * = k * k=1 F u k ,β k = k * k=1 I u k .β k (2.13)
hence the predicted value is given by: E(

y i ) = t∈T F * (A i (t), B i (t)
). We remark that each vector β k at each iteration may have different length dim(u k ).

Adaptation of SPICEFP to the partitioning of functional variables according to a non-linear scale

A linear partitioning of type (2.1) of the functional variables is not always suitable. Other types of partitioning may have to be chosen. Assuming for example that functional variable B requires a partitioning according to a logarithmic scale, the following breaks can be used instead of Equation (2.1):

L B (w) = B + 1 α B 1 + α B ( B -B) w-1 n B -1 , w = 1, . . . , n B + 1, (2.14)
where parameters α B > 0 and n B have to be set to determine L B (w). For a fixed n B , high α B value is related to high proportion of breaks close to B and vice versa. So, two partitioning parameters

Input : K; U A ; U B ; Γ; n λ ; Crit ∈ {AIC, BIC} Output: k * , (u k ) k≤k * ; (β k ) k≤k * Data : {y i , i = 1, . . . , n}; {A i (t); B i (t), i = 1, . . . , n, t ∈ T } foreach u ∈ U A . U B do for i ← 1 to n do for v ← 1 to n A do 4 for w ← 1 to n B do 5 C u i,(v,w) = Card {t ∈ T |A i (t) ∈ I u A (v), B i (t) ∈ I u B (w)} 6 end end X u [i, ] = t V ect(C u i ) end foreach γ ∈ Γ do
Construct D u,γ as shown in the equation (2.9)

end end for k ← 1 to K do foreach u ∈ U A . U B do
Center y and each joint modality in X u foreach γ ∈ Γ do 18

Find the solution path β u,γ (λ) = argmin Various types of breaks (chosen between Equations (2.1), (2.14) or user-defined) can be used for partitioning explanatory functional variables and the SPICEFP algorithm can be easily adapted. Experts in viticulture assume that the accumulation of chemical compounds affecting the quality of the grape berry is jointly influenced by these initial explanatory variables. This assumption is reinforced by results of [START_REF] Tarara | Berry temperature and solar radiation alter acylation, proportion, and concentration of anthocyanin in merlot grapes[END_REF], which underlined that the anthocyanin composition of Merlot grapes was influenced by a complex combined effect of berry temperature and solar irradiation. According to [START_REF] Varlet-Grancher | Mise au point : rayonnement solaire absorbé ou intercepté par un couvert végétal[END_REF], solar radiation can be characterized by three different quantifiers, including Photosynthetic Photon Flux Density (PPFD), measured in 10 -6 mol m -2 s -1 . It corresponds to the number of incident photons useful for photosynthesis, received per unit of time on a horizontal surface unit. Rows were roughly oriented south-north, and irradiance was separately measured on bunches located on the east and west side of the row. This measurement system therefore enabled the observation of different modalities of the couple (temperature, irradiance) affecting the grape berries. Temperature and Irradiance were recorded every twelve minutes throughout the maturation period when anthocyanins are known to accumulate.

β∈R n A n B 1 2 ||y -X u β|| 2 2 + λ||D u,γ
u k = u, β k = β u, γ e if β k = 0 n A n B then k * = k if k < K then 29 y ← y -X u k β k
So the experimental design contained: rows (1, 2, 3) × 2 sun exposure orientations (East, West) × 16 vine stocks = 96 "statistical individuals". Anthocyanin contents were measured weekly via the Ferari Index F I i for each individual i. This is a non-destructive measure of anthocyanin content of the bunches [START_REF] Agati | Assessment of anthocyanins in grape (vitis vinifera l.) berries using a noninvasive chlorophyll fluorescence method[END_REF]. The objective of our study is to understand how the couple (Temperature, Irradiance) acts on a weekly variation of the Ferari Index ∆F I.

Temperature and Irradiance are variables of different natures. Temperature is a variable whose variations are regular enough to be partitioned according to Equation (2.1). The Irradiance variable is partitioned according to (2.14), as explained in the following subsection.

Partitioning of the Irradiance variable

Observed on a one-day scale, PPFD increases exponentially from sunrise to a daily peak (observation time t max ), decreases until sunset, and remains almost constant until the next sunrise.

Irradiance primarily influences plant photosynthesis in a nonlinear way with a maximal reached at high irradiance. Therefore, the Irradiance variable was not partitioned according to a linear scale as proposed by Equation (2.1), but rather according to a logarithmic scale as in Equation (2.14).

The use of logarithmic transformation has consistently been used in the development of models involving solar radiation [START_REF] Salminen | A measuring system for estimating the frequency distribution of irradiance within plant canopies[END_REF], [START_REF] O'connor | Acacia karroo invasion of grassland: Environmental and biotic effects influencing seedling emergence and establishment[END_REF], [START_REF] Bergqvist | Sunlight exposure and temperature effects on berry growth and composition of cabernet sauvignon and grenache in the central san joaquin valley of california[END_REF].

Objective

The objective of our study was to identify the ranges of temperature and irradiance that jointly influence or not the accumulation of anthocyanins between sunrise and noon. This study is an appropriate application framework to test the SPICEFP method.

Our study is composed of two parts. First a simulation part is presented in section 4. Temperature and Irradiance, measured during the week of July 17 to 24, 2014, are the input variables of a model which simulates output variables to be predicted, using two different known β: one with two distinct patches of coefficients and another with a concentric gradient of coefficients (Table 2, column 1). This simulation study made it possible to evaluate the functioning and performance of the SPICEFP algorithm. Second, in section 5, the algorithm was tested on a complete dataset, obtained from July 24 to August 1, 2014, to understand the effect of Temperature and Irradiance interaction on ∆F I. Dataset and script are available online at https://forgemia.inra.fr/exploratory-penalizedregression/paper-script-and-data.git.

Simulation study

We present in this section simulations that help to better understand the SPICEFP characteristics. Remember that the approach must be able to identify an optimal partition or joint class intervals used as linear regressors.

Simulation design and SPICEFP setting

In order to carry out the simulations, a few steps were required:

• We considered the observations of temperature (A) and irradiance (B) in the Vine dataset obtained between sunrise and noon during the week of July 17 to 24, 2014.

• We then arbitrarily set a partition vector u 0 = (n A = 17, n B = 20, α B = 0.05).

• From A, B, and u 0 we constructed X u 0 using equations (2.2) and (2.5). X u 0 was computed based on observed data set in order to respect realistic frequencies for the joint class intervals.

• Based on dimensions of X u 0 , we drew coefficients β u 0 from a random distribution (values are available on git@forgemia.inra.fr) and computed the response variable of the simulation:

Y = X u 0 β u 0 + ε where ε ∼ N (0, σ 2 ε I). (4.1)
After that, we used the SPICEFP approach to make the estimation.

In this simulation study, we chose to simulate two different coefficient vectors. They are as followed:

• Coefficients β u 0 for simulation 1 is made of two distinct patches (column 1, row 1 of the Table 2). Two noise levels were used and presented in table 1.

• Coefficients β u 0 for simulation 2 is made of a concentric gradient of coefficients (column 1, row 2 of the Table 2). Two noise levels were also used and presented in Table 1.

= 1.50 (σ 2 ε /σ 2 Y = 0.03) σ ε = 2.5 (σ 2 ε /σ 2 Y = 0.08) Simulation 2 σ ε = 0.25 (σ 2 ε /σ 2 Y = 0.01) σ ε = 1 (σ 2 ε /σ 2 Y = 0.14)
We thus generated four simulation datasets (1 matrix X u 0 × 2 simulated coefficient vectors × 2 noise levels).

The inputs required by the algorithm were as follows:

• U A = {15, 16, . . . , 20}

• U B = {18, 19, . . . , 22}

• V B = {0.015, 0.05, 0.135, 0.405, 1.215}

• K = 2
• Γ = {0.0001, 0.05, 0.15, 0.45, 2, 8}

• n λ = 100

Simulation results

The results of the simulations are presented in Tables 2 to 5 (Tables 4 and5 are in the appendix).

Table 2 is related to the estimated coefficient vectors and Tables 4 and5 present respectively histograms of the residuals and scatter plots showing the quality of the estimates. For each of the four response variables (one per row), three estimations are provided (columns 3 to 5, Table 2), computed as follows:

• column 3: at the first iteration, an exploratory matrix X u 1 indexed by a partition vector u 1 is identified. The estimated coefficient at iteration 1 of SPICEFP is noted β 1 (see algorithm 1, line 25). The estimated response Y 1 = X u 1 β 1 and the residuals

ε 1 = Y -Y 1 are computed.
Figures in column 3 provide the visualization of β 1 . For a suitable visualization, the vector of coefficients is transformed into a matrix of dimension n A × n B , where n A et n B are the numbers of class intervals associated to u 1 .

• column 4: at the second iteration, an exploratory matrix X u 2 indexed by a partition vector u 2 is identified. The estimated coefficient at iteration 2 of SPICEFP is noted β 2 . The response variable used at the second iteration in the model is ε 1 , the residuals obtained at iteration

1. The estimated response Y 2 = X u 1 β 1 + X u 2 β 2 and the residuals ε 2 = ε 1 -X u 2 β 2 are computed. Figures in column 4 give the visualization of F (a, b) = F 1 (a, b) + F 2 (a, b)
where F is defined as in Equation (2.3). In the case of the visualizations in Table 4, 500 equally spaced values of a ∈ [A, Ā] and 1000 equally spaced values b ∈ [B, B] were used.

• column 5: the results presented in this column are slightly out of the scope of the SPICEFP approach and concerns only the outputs of iteration 1. The idea here is to consider an average of the 1% best models. The best models are defined in the sense of the information criterion.

Best coefficients relative to the different partitions available in these 1% best models are identified. Let's assume that there are n m . They are then indexed by u 1,(1) , u 1,(2) , . . . , u 1,(nm) .

The coefficients of interest can be noted β 1,(1) , β 1,(2) , . . . , β 1,(nm) . The estimated response is Y 1,(1:nm) = 1 nm nm m=1

X u 1,(m) β 1,(m) . All the selected models have the same weight in the computation of this average. The residuals can be obtained by ε 1,(1:nm) = Y -Y 1,(1:nm) .

Figures in column 5 Black color refers to never-observed joint modalities whereas white color refers to joint modalities with null estimated coefficients. The histograms of residuals related to each estimation are presented in Table 4. On top of these histograms, the ratio between the variance of residuals and the variance of response variable are provided. Table 5 shows the goodness of fit of each estimation, by presenting the slope of the regression "predicted versus simulated variables". The closer to 1 the slope is and the closer to 0 the ratio of variances is, the better the estimate is.

From the estimates provided by the algorithm (Table 2), we notice that SPICEFP effectively identifies the simulated zones of influence and assigns the right color to the coefficients: the graphics show two distinct areas, one with positive coefficients (red area), the other with negative coefficients (blue area). The approach tends to assign the same value (same colour) to groups of estimated coefficients, although with a gradient within the group. This behavior can be explained on the one hand by the Generalized Fused Lasso which penalizes the difference between two related coefficients and on the other hand by the variance estimate used for computing the information criterion. This variance was overestimated, which penalized the introduction of new coefficients into the model.

With respect to the noise level contained in the response variable, we note that the more noisy Y is, the more false positives are observed. Looking at the differences in the outputs of iterations 1 and 2, we observe that coefficients are added to those previously obtained. These new coefficients lead to improved estimation when the noise is low. When the noise is high, the β estimate of the second iteration has non-zero elements in some areas that are not relevant: for example for high temperature and irradiance values. The average of the 1% best models underestimates, in both simulations, the amplitude of the β values but, in simulation 2, it restores the β support much better than the estimators obtained in iterations 1 or 2.

5.

Modeling the evolution of a grape berry quality index

Methodology used for data analysis

We focus in this section on the modeling of the Ferari Index variation ∆F I from July 24 to August 1, 2014. We selected the n 1 = 32 individuals which have the highest contribution to final Ferari Index, with an initial index around 0.2 at the beginning of the week.

Most of the photosynthetic functioning of the grapevine happens in the morning. The time period between sunrise and noon is denoted T 1 below. The following variables and parameters are the input objects of SPICEFP.

• y i = ∆F I i , A i (t); B i (t); i = 1, . . . , n 1 ; t ∈ T 1 ;

• U A = U B = {10, 11, 12, . . . , 29, 30}

• 10 values were chosen between 0.0025 and 1 for α based on an exponential function α ∈ {0.0025, 0.0048, 0.0094, 0.0183, 0.0357, 0.0695, 0.1353, 0.2636, 0.5134, 1.0000}

• K = 3 • n λ = 20

Results

The results are presented in Table 3. The first column shows the estimated coefficients, the second column the histograms of residuals and the last column the scatter plots related to the goodness of fit. For this data set, SPICEFP stopped at the third iteration with a shutdown rule of zero coefficients.

The first row of Table 3 contains the results observed at the first iteration of SPICEFP. In terms of model quality, the slope is 0.558 and the residuals follow a normal distribution centered in 0. The visualization of the coefficients indicates conditions (irradiance < 100 µmol m -2 s -1 , temperature from 15°C to 33°C) that affect negatively the Ferari Index. The second row of Table 3 is relative to the results of the second iteration. Due to the fact that the coefficients retained in iteration 3 are all zero, the result of the approach is therefore that of row 2. The slope of model goodness of fit is 0.678 and the residuals follow a normal distribution centered in 0. The visualization of the coefficients indicates, in addition to what is observed at iteration 1, a small positive influence (compared to the amplitude of the negative effect) of the temperature below 30°C for an irradiance higher than 100 µmol m -2 s -1 .

When focusing on the average of the 1% best models (presented in the third row), we remark that the quality of this model (as indicated by the slope : 0.472) is not better than that of the models chosen by SPICEFP at iteration 1 or 2. However, it should be noted that the model obtained has more fused coefficients, which allows the identification of a border zone between the positive and negative zones of influence.

As an interpretation, we can note that, in the morning (sunrise to noon), for low irradiance values (< 100 µmol m -2 s -1 ), there is a range of temperature values that are not suitable for an increase of the Ferari index. On the contrary, a combination of irradiance values above 150 µmol m -2 s -1 and temperature below 30°C is suitable for increasing the Ferari index. The average of the coefficients shows, within each of the non-zero zones of influence, a relative variability of the coefficients, suggesting the importance of some temperature and irradiance amplitudes.

The slope coefficients of the goodness of fit of the models are away from 1. It should be reminded that the response variable is studied with respect to the variations of temperature and irradiance only between sunrise and twelve. There is also a lot of information hidden in the residuals.

Discussion

The SPICEFP approach is a scalar-on-function approach. The response variable is real and predictors are functional variables. The approach is based on a transformation of functional variables, which yields a contingency table. To construct this table, it is assumed that no data are missing and that the times of observation are identical for both functional variables. The method takes into account the multicolinearity resulting from the auto-correlations existing in the processes. Moreover, the constructed candidate explanatory matrices are not nested but they covered the same domain. That's why model and variable selection methods (Fused Lasso, information criterion selection) had to be adapted and generalized to the framework of SPICEFP. In the implementation of the approach, all the candidate explanatory matrices are evaluated at each iteration. This is not due to a re-sampling as is the case for bootstrap or bagging methods, but to the evaluation of different partitions.

Table 3: Visualization of the combined effects of irradiance and temperature on the Ferari index (From sunrise to twelve, I: joint modalities that have never been observed). Rows 1 and 2 present the results of iterations 1 and 2, respectively. The third row presents the average of the 1% best models obtained at iteration 1.

Estimation of coefficients Histogram of residuals Goodness of fit Iteration 1 Iteration 2

Mean of top 1% (iteration 1)

SPICEFP: a functional approach

In recent years, several studies in frequentist and Bayesian statistics, parametric as well as non-parametric, have focused on functional data analysis. The solutions developed can be used to achieve a wide range of goals (dimension reduction, regression, clustering, classification, etc.). They provide models that are predictive but often difficult to interpret. This lack of interpretability is partly due to the pre-processing step required on the functional data, which must be projected into bases of functions (splines, kernels...). The transformation of functional variables is a fundamental step, no matter if the variable is a response variable or a predictor.

The SPICEFP approach is primarily explanatory and not necessarily predictive. The first step is a transformation of the predictors into categorical variables. The choice of this transformation was motivated by the potential to interpret the results while considering the hypothesis of a joint influence of the predictors on the response variable. Each partition is a collection of 2D intervals, see §2.1. We constructed linear regressors on the basis of indicator functions associated to these 2D intervals, see Equation (2.4). This basis of indicator functions is not common compared to polynomial basis, Fourier basis, wavelet basis, etc., [START_REF] Wang | Functional data analysis[END_REF]. In a previous work [START_REF] Marx | Multidimensional penalized signal regression[END_REF], the authors have focused on multidimensional penalized signal regression, a single surface was estimated with smooth regression coefficient using B-spline tensor products. In our case, the indicator functions facilitate the interpretation of the results, provided that both functional predictors (T in grape berry dataset) are discretized over the same set of equidistant observation times. This constraint can be released with usual pretreatment such as:

• imputation of missing data [START_REF] Stekhoven | MissForest-non-parametric missing value imputation for mixed-type data[END_REF][START_REF] Josse | missmda: A package for handling missing values in multivariate data analysis[END_REF] • interpolation, smoothing [START_REF] Ramsay | Functional Data Analysis with R and MAT-LAB[END_REF] or restriction of the functional variable to an identical set of observation times. Indeed, functions have uncountable supports.

The use of 2D class intervals implies also the assumption that the structure of the underlying process does not change over the observation variable (time in the grape berry dataset). This is what we can call hypothesis of stationarity. For example, in our use case on grapevine, this hypothesis of stationarity requires to work at the scale of a week but also to split the day (working with observations obtained in the morning (sunrise to noon)). When analyzing the data in section 5, we assumed that the underlying process is time-invariant in the mornings of the week under consideration.

The method is design-dependent: it will not be able to properly estimate the β coefficient in an area with little or no data. It is therefore necessary to have data of (A, B) in areas where there is potentially something happening. There is also a limitation due to the curse of dimensionality [START_REF] Giraud | Introduction to High-Dimensional Statistics[END_REF]: "If the number n of observations remains fixed while the dimension p of the observations increases, the observations get rapidly very isolated and local methods cannot work.". The more fluctuations in many directions, the more data will be needed. Our approach will work best if the shapes are simple.

SPICEFP: an iterative approach

The SPICEFP approach can be summarized as follows: first, the functional predictors are transformed into a set of candidate matrices X u (and D u,γ ) for each (u, γ). Then we iterate the following steps: estimation of the coefficients for each of the candidate matrices, evaluation and selection based on the chosen validation information criteria with respect to (u, γ) and then checking of the shutdown conditions. The number of iterations performed is controlled by the shutdown conditions. The residuals obtained at iteration k are used as response variable at iteration k + 1. This iterative procedure has been implemented in order to:

• extract knowledge still available in the residuals with the risk of overestimating the number of non-zero coefficients: This risk is linked to several factors such as:

1. the noise level in the response variable: the higher the noise level in the response variable, the higher the risk of model misspecification in the SPICEFP approach [START_REF] Gustafson | On the simultaneous effects of model misspecification and errors in variables[END_REF]. This situation is well illustrated by the results of the simulations at iteration 2 (row 1 vs row 2 and row 3 vs row 4, column 4, Table 2).

2. the shape of the area of influence to be identified: the number of connected components to cover the shape increases when the shape is not a combination of rectangles. For a shape composed of rectangles of different sizes, the iterative approach allows changing the scale of the grid to get a refined estimation. This tolerance to the overestimation of the coefficients is made in return for the identification of easily interpretable areas of influence. This contributes to the explanatory character of the approach. The number of non-zero coefficients is also related to the value of the γ parameter presented in the Equation (2.7) which controls the ratio between parsimony and fusion of the coefficients. The lower the parameter γ (conditionally to the data), the larger the estimated area of influence. This increases the risk of getting unnecessary /untrue non-zero coefficients.

1 )

 1 with A ∈ R and Ā ∈ R the minimum and maximum scalar values in A. The bins used for partitioning all(A i ) i=1...n are I A (v) = [L A (v), L A (v + 1)[, v = 1, . . . , n A .The partition is the same for all i, i = 1, . . . , n. Using the same approach for partitioning the second explanatory variable B, we obtain n B + 1 breaks L B (w). B ∈ R and B ∈ R are the minimum and maximum scalar values in B. The bins used for partitioning all (B i ) i=1,...,n are I B (w) = [L B (w), L B (w + 1)[, w = 1, . . . , n B . The numbers of class intervals n A and n B have to be set to compute the breaks L A (v) and L B (w), v = 1, . . . , n A , w = 1, . . . , n B .

  w) can also be interpreted as a discrete approximation of the density of the time spent by the individual i with variable A around L A (v) and variable B around L B (w). Part 1 of Figure 2.1 shows the transformation of the functional explanatory variables A and B for the fixed u = (4, 3). Note that, for a fixed u, the (I u A

Figure 2 . 1 :

 21 Figure 2.1: Transformation of both functional explanatory variables for the SPICEFP approach

Figure 2 . 2 :

 22 Figure 2.2: Example of β u,γ coefficient values with 4 connected components (here u = (9, 8)). Two coefficients β u,γ v,w and β u,γ v ,w are linked via the D u,γ matrix if and only if |v -v| = 1 or (exclusive) |w -w| = 1 (i.e., if and only if |v -v| + |w -w| = 1). Let us consider that the

β|| 1 19 20 AICσ 2

 1202 Select n λ equally spaced couples (λ e , β u,γ e ) on the log scale with respect to λ over the solution path and compute Crit u,γ e for e = 1 . . . n λ : ||y -X u β u,γ e || 2 2 + log(n)Q u,γ e end end ( u, γ, e) ← argmin u∈U A .U B , γ∈Γ, 1≤e≤n λ Crit u,γ e

X u k β k Algorithm 1 :

 1 SPICEFP algorithm (n B and α B ) should be optimized for partitioning B. In this case, the partition vector is written u = (n A , n B , α B ). Let V B be the set containing the possible values of α B . V B is an additional input to the algorithm. The only change in the core of the algorithm is to replace u ∈ U A . U B by u ∈ U A . U B . V B (lines 1 to 15).

3 .

 3 Use Case: Grapevine dataset 3.1. Data presentation Data were collected during an experiment conducted in a vineyard of the Institut Agro campus at Montpellier in 2014 (Syrah vines). The aim was to study the influence of the micro-climate (temperature, solar irradiation) at the grape level on the anthocyanin content of the berries.

  The experimental plot was made of three rows of vines within the vineyard, each with eight vines equipped with open-top chambers to warm the basis of the plant, and eight under control conditions (without open-top chambers). The chambers were made up of 2 translucent polycarbonate panels placed on the ground at about 10cm below the bunches on each side of the vines and inclined to form a mini two-pitched greenhouse roof open at its ridge. The greenhouse effect created during the day in the chambers generated a flow of warm air that escaped through the open top, raising the temperature of the bunches by 2 to 3 °C, mimicking global warming. The microclimate at bunch level was recorded through the measurement of temperature and irradiance.

show the visualization of 1 nm nm m=1 F 1 ,

 m=11 Figures in column 5 show the visualization of 1 nm

Table 1 :

 1 Noise simulation design

	Low noise	High noise
	Simulation 1 σ ε

Table 2 :

 2 Simulation results: estimation with the SPICEFP algorithm of two coefficient functions β with two types of noise (high and low). I: Joint modalities that have never been observed.

		Explanatory variables
		Simulated coefficients
	Simulation 1	Simulation 2
		Responses
	Low noise	High noise	Low noise	High noise
		Iteration 1
		Iteration 2
		Mean of top 1% (iteration 1)

Table 4 :

 4 Histogram of residuals Joint modalities that have never been observed (no t counted for these joint modalities for all individuals)

	Simulations

I:

Table 5 :

 5 Quality of the estimate Joint modalities that have never been observed (no t counted for these joint modalities for all individuals)

	Simulations

I:

• Explore the space of observations by using different partitions/pavements (fine scale, large scale): from the second iteration onwards, the candidate matrix selected covers the same space as the previous one, with the difference that it can cover it with different widths. Thus new areas of influence can be identified, or existing ones reinforced, etc. There is therefore no independence between the different candidate matrices because they are all constructed from the observed predictors. Even if the response of the next iteration corresponds to the residuals of the previous iteration, this approach is not orthogonal.

Eventually, the shutdown condition of SPICEFP is rather natural. The approach stops when the vector of coefficients at one iteration is null: in this case, the residuals of the model are equivalent to the response variable. This criterion is internal to the approach.

Acknowledgments

The authors dedicate this work to their colleague Eric Lebon, who passed away. He participated in the very fruitful discussions that initiated the present work, and provided us with the data. Data