Baptiste Gregorutti
email: baptiste.gregorutti@safety-line.fr

Bertrand Michel
email: bertrand.michel@upmc.fr

Philippe Saint-Pierre
email: philippe.saintpierre@upmc.fr

Grouped variable importance with random forests and application to multiple functional data analysis

Keywords: Random forests, functional data analysis, group permutation importance measure, group variable selection. 2010 MSC

The selection of grouped variables using the random forest algorithm is considered. First a new importance measure adapted for groups of variables is proposed. Theoretical insights into this criterion are given for additive regression models. Second, an original method for selecting functional variables based on the grouped variable importance measure is developed. Using a wavelet basis, it is proposed to regroup all of the wavelet coefficients for a given functional variable and use a wrapper selection algorithm with these groups. Various other groupings which take advantage of the frequency and time localization of the wavelet basis are proposed. An extensive simulation study is performed to illustrate the use of the grouped importance measure in this context. The method is applied to a real life problem coming from aviation safety.

Introduction

In the high dimensional setting, identification of the most relevant variables has been the subject of much research during the last two decades [START_REF] Guyon | An introduction to variable and feature selection[END_REF]. For linear regression, the lasso method [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] is widely used. Many variable selection procedures have also been proposed for non linear methods. In the context of random forests [START_REF] Breiman | Random forests[END_REF], it has been shown that the permutation importance measure is an efficient tool for selecting variables (Díaz-Uriarte and Alvarez de [START_REF] Díaz-Uriarte | Gene selection and classification of microarray data using random forest[END_REF][START_REF] Genuer | Variable selection using random forests[END_REF][START_REF] Gregorutti | Correlation and variable importance in random forests[END_REF].

In many situations such as medical studies and genetics, groups of variables can be clearly identified and it is of interest to select groups of variables rather than to select them individually [START_REF] He | Stable feature selection for biomarker discovery[END_REF]. Indeed, interpretation of the model may be improved along with the prediction accuracy by grouping the variables according to a priori knowledge about the data. Furthermore, grouping variables can be seen as a solution to stabilize variable selection methods. In the linear setting, and more particularly for linear regression, the group lasso has been developed to deal with groups of variables, see for instance Yuan and Lin (2006a). Group variable selection has also been proposed for kernel methods [START_REF] Zhang | Variable selection for the multicategory svm via adaptive sup-norm regularization[END_REF] and neural networks [START_REF] Chakraborty | Selecting useful groups of features in a connectionist framework[END_REF]. As far as we know, this problem has not been studied for the random forest algorithm introduced by [START_REF] Breiman | Random forests[END_REF]. In this paper, we adapt the permutation importance measure for groups of variables in order to select groups of variables in the context of random forests.

The first contribution of this paper is a theoretical analysis of the grouped variable importance measure. Generally speaking, the grouped variable importance does not reduce to the sum of the individual importances and may even be quite unrelated to it. However, in more specific models such as additive regression ones, we derive exact decompositions of the grouped variable importance measure.

The second contribution of this work is an original method for selecting functional variables based on the grouped variable importance measure. Functional Data Analysis (FDA) is a field in statistics that analyzes data indexed by a continuum. In our case, we consider data providing information about curves varying over time [START_REF] Ramsay | Functional Data Analysis[END_REF][START_REF] Ferraty | Nonparametric Functional Data Analysis: Theory and Practice[END_REF]Ferraty, 2011). One standard approach in FDA consists in projecting the functional variables onto a finite dimensional space spanned by a functional basis. Classical bases in this context are splines, Fourier, wavelets or Karhunen-Loève expansions, for instance. Most of the papers about regression and classification methods for functional data consider only one functional predictor; references include [START_REF] Cardot | Functional linear model[END_REF][START_REF] Cardot | Spline estimators for the functional linear model[END_REF]; Rossi et al. (2006); [START_REF] Cai | Prediction in functional linear regression[END_REF] for linear regression methods, [START_REF] Amato | Dimension reduction in functional regression with applications[END_REF]; [START_REF] Araki | Functional logistic discrimination via regularized basis expansions[END_REF] for logistic regression methods, [START_REF] Górecki | A comparison of tests for the one-way anova problem for functional data[END_REF] for ANOVA problem, [START_REF] Biau | Functional classification in hilbert spaces[END_REF]; [START_REF] Fromont | Functional classification with margin conditions[END_REF] for k-NN algorithms and Rossi andVilla (2006, 2008) for SVM classification. The multiple FDA problem, where p functional variables are observed, has been less studied. Recently, [START_REF] Matsui | Variable selection for functional regression models via the regularization[END_REF] and [START_REF] Fan | Functional additive regression[END_REF] have proposed solutions to the linear regression problem with lasso-like penalties. The logistic regression case has been studied by [START_REF] Matsui | Variable and boundary selection for functional data via multiclass logistic regression modeling[END_REF]. Classification based on several functional variables has also been considered using the CART algorithm [START_REF] Poggi | Classification supervisée en grande dimension. application à l'agrément de conduite automobile[END_REF] and SVM [START_REF] Yang | A supervised feature subset selection technique for multivariate time series[END_REF][START_REF] Yoon | Feature subset selection on multivariate time series with extremely large spatial features[END_REF].

We propose a new approach for multiple FDA using random forests and the grouped variable importance measure. Indeed, various groups of basis coefficients can be proposed for a given functional decomposition. For instance, one can choose to regroup all coefficients of a given functional variable. In this case, the selection of a group of coefficients corresponds to the selection of a functional variable. Various other groupings are proposed for wavelet decompositions. For a given family of groups, we adapt the recursive feature elimination algorithm [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF] which is particularly efficient when predictors are strongly correlated [START_REF] Gregorutti | Correlation and variable importance in random forests[END_REF]. In the context of random forests, this backward-like selection algorithm is guided by the grouped variable importance. Note that by regrouping the coefficients, the computational cost of the algorithm is drastically reduced compared to a backward strategy that would eliminate only one coefficient at each step.

An extensive simulation study illustrates the application of the grouped importance measure for FDA. The method is then applied to a real life problem coming from aviation safety. The aim of this study is to explain and predict landing distances. We select the most relevant flight parameters regarding the risk of long landings, which is a major issue for airlines.

The group permutation importance measure is introduced in Section 2. Section 3 deals with multiple FDA using random forests and the grouped variable importance measure. The application to flight data analysis is presented in Section 4. Note that additional experiments about the grouped variable importance are given in Appendix B. In order to speed up the algorithm, the dimension of the data can be reduced in a preprocessing step. In Appendix C, we propose a modified version of a well-known shrinkage method [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] that simultaneously shrinks to zero the coefficients of the observed curves of a functional variable.

The grouped variable importance measure

Let Y be a random variable in R and X = (X 1 , . . . , X p ) a random vector in R p . We denote by f (x) = E[Y |X = x] the regression function. Let Var(X) and Cov(X) denote the variance and variancecovariance matrices of X.

The permutation importance introduced by [START_REF] Breiman | Random forests[END_REF] measures the accuracy of each variable X j for predicting Y . It is based on the elementary property that the quadratic risk E (Yf (X)) 2 is the minimum error for predicting Y knowing X. The formal definition of the variable importance measure of X j is:

I(X j ) := E Y -f (X (j) ) 2 -E (Y -f (X)) 2 , (1)
where X (j) = (X 1 , . . . , X j , . . . , X p ) is a random vector such that X j is an independent replicate of X j which is also independent of Y and of all other predictors. This criterion evaluates the increase of the prediction error after breaking the link between the variable X j and the outcome Y (see [START_REF] Zhu | Reinforcement Learning Trees[END_REF] for instance).

In this paper, we extend the permutation importance to groups of variables. Let J = (j 1 , . . . , j k ) be a k-tuple of increasing indices in {1, . . . , p}, with k ≤ p. We define the permutation importance of the sub-vector X J = (X j1 , X j2 , . . . , X j k ) of predictors by

I(X J ) := E Y -f (X (J) ) 2 -E (Y -f (X)) 2 ,
where X (J) = (X 1 , . . . , X j1 , X j1+1 , . . . , X j2 , X j2+1 , . . . , X j , X j +1 , . . . , X p ) is a random vector such that X J = (X j1 , X j2 , . . . , X j k ) is an independent replicate of X J , which is also independent of Y and of all other predictors. We call this quantity the grouped variable importance since it only depends on which variables appear in X J . By abuse of notation and ignoring the ranking, we may also refer to X J as a group of variables.

Decomposition of the grouped variable importance

Let X J be a subgroup of variables from the random vector X. Let X J denote the group of variables that does not appear in X J . Assume that we observe Y and X in the following additive regression model:

Y = f (X) + ε = f J (X J ) + f J (X J ) + ε, (2)
where the f J and f J are two measurable functions and ε is a random variable such that [START_REF] Gregorutti | Correlation and variable importance in random forests[END_REF] on the permutation importance of individual variables can be extended to the case of a group of variables.

E[ε|X] = 0 and E[ε 2 |X] is finite. Results in
Proposition 1. Under model (2), the importance of group J satisfies

I(X J ) = 2 Var [f J (X J )] .
The proof is given in Appendix A. The next result gives the grouped variable importance for more specific models. It can be easily deduced from Proposition 1.

Corollary 1. Assume that we observe Y and X in model (2) with

f J (x J ) = j∈J f j (x j ), (3)
where the f j are measurable functions and x J = (x j ) j∈J .

1. If the random variables (X j ) j∈J are independent, then

I(X J ) = 2 j∈J Var (f j (X j )) = j∈J I (X j ) .
2. If for any j ∈ J, f j is a linear function such that f j (x j ) = α j x j , then

I(X J ) = 2α J Cov(X J )α J , (4)
where α J = (α j ) j∈J .

If f is additive and if the variables of the group are independent, the grouped variable importance is nothing more than the sum of the individual importances. As affirmed in the second point of Corollary 1, this property is lost as soon as the variables in the group are correlated. The experiments presented in Appendix B allow us to compare the grouped variable importance with the individual importances in various models. In essence, these experiments suggest that in general, the grouped variable importance is not comparable with the sum of the individual importances. This is not surprising since the grouped variable importance is a more accurate measure of the importance of a group of variables than a simple sum of the individual importances.

Next, let us introduce a rescaled version of the grouped variable importance:

I nor (X J ) := 1 |J| I(X J ).
To see why this quantity makes sense, let us consider the situation where two groups of variables have almost the same group permutation importances. Then, one would prefer to select first the smaller group in order to obtain a sparser set of predictors. We thus propose to normalise I(X J ) by an increasing function of group size |J| so as to compare permutation importances of groups of variables. Moreover, Corollary 1 tells us that normalisation by |J| is reasonable since it corresponds to comparing the means of the individual permutation importances over each group in the case where the variables are independent.

In Section 3, we propose a variable selection algorithm based on the grouped variable importance. This algorithm uses the rescaled version to take into account group sizes in the selection process. More generally, we choose to work with the rescaled version when comparing groups of variables of different sizes.

Grouped variable importance and random forests

Classification and regression trees are well-performing techniques for estimating f . A popular method in this context is the CART algorithm of [START_REF] Breiman | Classification and regression trees[END_REF]. Though efficient, tree methods are also known to be unstable insofar as a small perturbation of the training sample may change radically the predictions. To counter this, [START_REF] Breiman | Random forests[END_REF] introduced random forests as a substantial improvement to decision trees. The permutation importance measure was also introduced in this seminal paper. We now recall how individual permutation importances can be estimated with random forests before presenting a natural extension to the estimation of grouped variable importances.

Assume that we observe n i.i.d. replicates D n = {(X 1 , Y 1 ), . . . , (X n , Y n )} of (X, Y ). The random forest algorithm consists in aggregating a collection of random trees as in the bagging method, also proposed by [START_REF] Breiman | Bagging predictors[END_REF]; trees are built over M bootstrap samples D 1 n , . . . , D M n of the training data D n . Instead of the CART algorithm, a subset of variables is randomly chosen for the splitting rule at each node. Each tree is then fully grown or grown until each node is pure, and not subsequently pruned. The resulting learning rule is the aggregation of all of the tree-based estimators, denoted by f1 , . . . , fM . In the regression setting, aggregation is based on the average of the predictions.

For any m ∈ {1, . . . , M }, let Dm n := D n \ D m n be the corresponding out-of-bag sample. The risk of fm is estimated on the out-of-bag sample as follows:

R( fm , Dm n ) = 1 | Dm n | i:(Xi,Yi)∈ Dm n (Y i -fm (X i )) 2 .
Let Dmj n be the permuted version of Dm n obtained by randomly permuting the variable X j in each out-of-bag sample Dm

n . An estimation of the permutation importance measure of variable X j is then given by

Î(X j ) = 1 M M m=1 R( fm , Dmj n ) -R( fm , Dm n ) . (5)
This random permutation mimics the replacement of X j by X j in (1) and breaks the link between X j and Y and the other predictors. We now extend the method for estimating the permutation importance of a group of variables X J . For any m ∈ {1, . . . , M }, let DmJ n be the permuted version of Dm n obtained by randomly permuting the group X J in each out-of-bag sample Dm n . Note that the same random permutation is used for each variable X j of the group. In this way the (empirical) joint distribution of X J is left unchanged by the permutation whereas the link between X J and Y and the other predictors is broken. The importance of

X J can be estimated by Î(X J ) = 1 M M m=1 R( fm , DmJ n ) -R( fm , Dm n ) . (6)
We define Înor (X J ) as the rescaled version of this estimate. In the following section, we use the grouped variable importance as a criterion for selecting features in the context of multiple functional data.

Multiple functional data analysis using grouped variable importance

In this section, we consider an application of grouped variable selection for multiple functional regression with scalar response Y . Each covariate X 1 , . . . , X p takes its values in the Hilbert space L 2 ([0, 1]) equipped with the inner product

f, g L 2 = f (t)g(t)dt, for f, g ∈ L 2 ([0, 1]).
One common approach of functional data analysis is to project the variables onto a finite dimensional subspace of L 2 ([0, 1]) and to use the basis coefficients in a learning algorithm [START_REF] Ramsay | Functional Data Analysis[END_REF]. For instance, the wavelet transform is widely used in signal processing and for nonparametric function estimation (see for instance [START_REF] Antoniadis | Wavelet estimators in nonparametric regression: A comparative simulation study[END_REF]. Unlike Fourier bases and splines, wavelets are localized both in frequency and time.

For j ≥ 0 and k = 0, . . . , 2 j -1, define a sequence of functions φ jk (resp. ψ jk ), obtained by translations and dilatations of a compactly supported function φ (resp. ψ), called a scaling function (resp. wavelet function). For any j 0 ≥ 0, the collection

B = {φ j0k , k = 0, . . . , 2 j0 -1} ∪ {ψ jk , j ≥ j 0 , k = 0, . . . , 2 j -1}
forms an orthonormal basis of L 2 ([0, 1]) (see for instance [START_REF] Percival | Wavelet Methods for Time Series Analysis[END_REF]). Then, a function s ∈ L 2 ([0, 1]) can be decomposed as

s(t) = 2 j 0 -1 k=0 s, φ j0k L 2 φ j0k (t) + j≥j0 2 j -1 k=0 s, ψ jk L 2 ψ jk (t). ( 7
)
The first term in Equation ( 7) is the smooth approximation of s at level j 0 while the second term is the detail part of the wavelet representation. We assume that each covariate X is observed on a fine sampling grid t 1 , . . . , t N with t = N . Note that a wavelet decomposition of X can also be given in a similar form as in (7). For j 0 = 0, we have

X(t ) = ζφ 00 (t ) + J-1 j=0 2 j -1 k=0 ξ jk ψ jk (t ), (8)
where J := log 2 (N ) is the maximal number of wavelet levels and ζ and ξ jk are respectively the scale and wavelet coefficients of the discretized curve X at position k for resolution level j. These empirical coefficients can be efficiently computed using the discrete wavelet transform described in Chapter 4 of [START_REF] Percival | Wavelet Methods for Time Series Analysis[END_REF].

For a given wavelet basis, we introduce the wavelet support at time t as the set of all indices of wavelet functions that are non null at t: interval T , we define the wavelet support of T by

S(t) = {(j, k) : ψ jk (t) = 0}.
S(T ) = {(j, k) : ψ jk (t) = 0, ∀t ∈ T } = t∈T S(t).
This set corresponds to wavelet functions localized on the interval T .

Grouped variable importance for functional variables

In this section, we show how the grouped variable importance can be fruitfully used for comparing the importance of different wavelet coefficients in the context of functional prediction. Remember that p functional covariates X 1 , . . . , X p are observed, along with a scalar response Y . For the sake of simplicity, covariates are decomposed on the same wavelet basis B but the methodology presented above could be also adapted with a specific basis for each covariate. For any u ∈ {1, . . . , p}, let W u = (ζ u , ξ u jk ) jk be the random vector composed of the wavelet coefficients of the functional variable X u .

Groups of wavelet coefficients

Wavelet coefficients are characterised by their frequency, time location and the functional variables they describe. Consequently, they can be grouped in many ways. We give below a nonexhaustive list of groups for which we are interested in computing the importance:

• A group related to a variable. The vector W u defines the group G(u).

• A group related to a frequency level of a variable. For a fixed variable X u , the group is composed of the wavelet coefficients of frequency level j: G(j, u) := {ξ u j,1 , . . . , ξ u j,2 j -1 }.

• A group related to a frequency level. The group is composed of all the wavelet coefficients of frequency level j for all the variables:

G(j) := u=1,...,p

G(j, u).

• A group related to a given time. Define the group of "active" wavelet coefficients associated to a given time t by

G(t) := u=1,...,p {ζ u } ∪ u=1,...,p, (j,k)∈S(t) {ξ u jk }.
Depending on the size of the support of φ and ψ, the group G(t) may be very large. For instance with a Daubechies wavelet basis with two vanishing moments, in Figure 1 the group G(t) is composed of the darkest points of the row corresponding to time t.

• A group related to a time interval. Let [a, b] be a time interval. The group of "active" wavelet coefficients associated with

[a, b] is G([a, b]) := t∈[a,b] G(t).
Many other groupings could be proposed. For instance, one could regroup pairs of correlated variables, or consider a group composed of the wavelet coefficients taken in an interval of frequencies, a group related to a given time and a fixed variable, etc. By computing the importance of such groups, one directly obtains a rough detection of the most important groups of coefficients for predicting Y . When grouping by frequency levels or time locations, all groups do not have equal sizes. As explained in Section 2, it is preferable to use the rescaled version of the grouped variable importance in order to compensate the effect of group size in the grouped variable importance measure.

Grouped variable selection

We now propose a more elaborate method for selecting groups of coefficients. The selection procedure is based on the Recursive Feature Elimination (RFE) algorithm proposed by [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF] in the context of support vector machines. In this paper, we propose a random forests version of the RFE algorithm which is guided by grouped variable importance. The procedure is summarized in Algorithm 1. This backward grouped elimination approach produces a collection of nested subsets of groups. The selected groups are obtained by minimizing the validation error computed in step 2.

Algorithm 1 Grouped Variable Selection

1: Train a random forest model 2: Compute the error using a validation sample 3: Compute the grouped variable importance measure 4: Eliminate the least important group of variables 5: Repeat 1-4 until no further groups remain This algorithm is motivated by the results of our previous work [START_REF] Gregorutti | Correlation and variable importance in random forests[END_REF] about variable selection using the permutation importance measure for random forests. Strong correlation between predictors has a strong impact on the permutation importance measure. It was also shown in this previous paper that when predictors are strongly correlated, the RFE algorithm provides a better performance than the "non-recursive" strategy (NRFE) that computes the grouped variable importance just once and does not recompute the importance at each step of the algorithm. In the present paper, we continue this study by adapting the RFE algorithm for the grouped variable importance measure. We give below two applications of this algorithm.

• Selection of functional variables. Each vector W u defines a group G(u) and the goal is to perform a grouped variable selection over the groups G(1), . . . , G(p). The selection allows us to identify the most relevant functional variables.

• Selection of wavelet levels for a given functional variable. For a fixed u, we make a selection over the groups G(j, u) to identify the frequency levels which yield predictive information.

Remark. Algorithm 1 for grouped variable selection is appropriate for groups defining a partition over the wavelet coefficients. This is not the case for groups related to time locations. The algorithm can not easily be adapted to these groups because most wavelet coefficients belong to several groups and the elimination of a whole group may not be an efficient strategy. For instance the coefficient ζ u , which approximates the smooth part of the curves and which is usually a good predictor, is common to all times t.

Experiments

In this section, we present various numerical experiments for illustrating the interest of grouped variable importance for analyzing functional data. We first describe the simulation set-up.

Presentation of the general simulation design

The experiments presented below consider one or several functional covariates for predicting an outcome Y . Except for the second simulation in Experiment 1 which is presented in detail in the next section, the functional covariates are defined as functions of Y .

First, an n-sample of the outcome variable Y is simulated from a given distribution specified for each experiment. The realization of a functional covariate X (denoted by X u when there are several) is an nsample of independent discrete time random processes X i = (X i (t )) =1,...,N , for all i ∈ {1, . . . , n} under a model of the form

X i (t ) = s(t , Z i ) + σε i, , = 1, . . . , N, (9)
where the ε i, are i.i.d standard Gaussian random variables, and t = N . The random variable Z i is correlated with Y i and will be specified for each experiment; it is equal to the outcome variable Y i for most of the experiments. The functional covariates are actually simulated in the wavelet domain from the following model: for i = 1, . . . , n, j = 0, . . . , J -1 and k = 0, . . . , 2 j -1,

ζ i = ω 0 + h ζ (Z i ) + ση iζ , (10)
and

ξ ijk = ω jk + h jk (Z i ) + ση ijk if j ≤ j , k ∈ {0, . . . , 2 j -1}, 0 if j < j ≤ J -1, ( 11
)
where j is the highest wavelet level of the signal. The random variables η ijk and η iζ are i.i.d. standard Gaussian. The "signal" part of Equation ( 11) is the sum of a random coefficient ω jk whose realisation is the same for all i, and a link function h jk . The coefficients ω 0 and ω jk in ( 10) and ( 11) are simulated as follows:

ω jk ∼ N 1 (0, τ 2 j ), if j ≤ j , k ∈ {0, . . . , 2 j -1}, ω 0 ∼ N 1 (3, 1),
where τ j = e -(j-1) . Note that the standard deviation τ j decreases with j and thus less noise is added to the first wavelet levels. The link function h jk describes the link between the wavelet coefficients ξ ijk and the variable Z (or the outcome variable Y ). Two different link functions are considered in the experiments:

• a linear link: h jk (z) = θ jk z, • a logistic link: h jk (z) = θ jk 1 + e -z ,
where the coefficients θ jk parametrize the strength of the relationship between Z and the wavelet coefficients. The n discrete processes X 1 , . . . , X n are simulated according to (10) and ( 11) before applying the inverse wavelet transform. We choose a Daubechies wavelet filter with four vanishing moments to simulate the observations. We use the same basis for the projection of the functional observations.

Experiment 1: Detection of important time intervals

In this first experiment, we illustrate the use of the grouped variable importance for detection of the most relevant time intervals. We simply estimate the importance of time intervals without applying Algorithm 1 (see Remark 3.1). We only consider one functional covariate X since it will be sufficient to illustrate the method. Let T = [t 50 , t 55 ], we propose two simulation designs for which the outcome Y is correlated to the signal X on the interval T .

• Simulation 1. For this first simulation, we follow the general simulation design presented before by considering linear link functions h for all wavelet coefficients belonging to the wavelet support S(T ).

The outcome variable Y is simulated from a Gaussian distribution N 1 (0, 3). We simulate the wavelet coefficients as in ( 11) and the scaling coefficients as in ( 10) with Z = Y . We take linear link functions and set n = 1000, σ = 0.01, N = 2 8 and thus J = 8. We take j = 7 which means that even the wavelet coefficients of highest level j = 7 are not Dirac distributions at zero. The wavelet coefficients and the scaling coefficient are generated as follows: for any j ∈ {0, . . . , J -1} and any k ∈ {0, . . . , 2 j -1},

ξ ijk = ω jk + Y i + ση ijk , if (j, k) ∈ S(T ) ω jk + ση ijk , otherwise, (12)
and

ζ i = ω 0 + Y i + ση iζ , (13)
for i = 1, . . . , n.

• Simulation 2. In contrast to the previous simulation, we first simulate the functional variable X and then simulate the outcome variable Y as a function of X. The functional variable is simulated in the wavelet domain according to ( 10) and ( 11) with h jk = h ζ = 0 for all j and k. We also take n = 1000, σ = 0.01, N = 2 8 and j = 7. By applying the wavelet inverse transform for any i, we obtain an n-sample of discrete time random processes X i = (X i (t )) =1,...,N . Figure 2 displays a set of ten of these processes. The outcome variable Y is then obtained via

Y i = 1000 |T | t ∈T |X i (t ) -X i (t -1 )|.
Thus, Y i is a measure of the oscillations of the curve X i over the interval T .

The aim is to detect T using the grouped variable importance. In both cases, the grouped variable importance I(G(t)) is evaluated at 50 equally spaced time points. Figure 3 displays the importance of the time points, averaged over 100 iterations. The first and third quartiles are also represented in order to highlight estimation variability. In the two cases, the importance estimation makes it possible to detect T . Note that the detection problem is tricky in the second case because the link between Y and the wavelet coefficients is complex here. Consequently the estimated importances are low and the important intervals are difficult to detect. Experiment 2: Selection of wavelet levels This simulation deals with the selection of wavelet levels for one functional variable. We follow the general simulation design presented above. The outcome variable Y is simulated from a Gaussian distribution N 1 (0, 3). We simulate the wavelet coefficients as in (11) and the scaling coefficients as in ( 10) with Z = Y . We perform two types of simulation:

• In the first case we use linear link functions:

h ζ (y) = 0.1y and h jk (y) = θ j y if j ≤ 3, k ∈ {0, . . . , 2 j -1}, 0 otherwise,
where the θ j decrease linearly from 0.1 to 0.01.

• For the second simulation, we use logistic link functions:

h ζ (y) = 0.1 1 + e -y and h jk (y) =    θ jk 1 + e -y if j ≤ 3, k ∈ {0, . . . , 2 j -1}, 0 otherwise,
where the θ j decrease linearly from 0.1 to 0.01.

We fix n = 1000, σ = 0.05, N = 2 8 (thus J = 8) and j = 7 in both cases. The aim is to identify the most relevant wavelet levels for the prediction of Y using the grouped importance. We regroup the wavelet coefficients by wavelet levels: for j ∈ {0, . . . , J -1},

G(j) = ξ jk , k ∈ {1, . . . , 2 j -1} and G ζ = {ζ}.
We apply Algorithm 1 with these groups. As the group sizes are different, the rescaled grouped importance criterion given in Section 2.2 is used.

The trials are both repeated 100 times. Figures 4 and5 respectively give the results for the linear and logistic links. Let us first look at the trials with a linear link. The boxplots of the grouped permutation importances at the first step of the algorithm over the 100 trials are given in Figure 4a. The fifth group G(3) being not strongly correlated with Y , its importance is close to zero. It is selected 40 times out of the 100 simulations (Fig. 4c) whereas G ζ , G(0), G(1) and G(2) are almost always selected. The other groups are not correlated with Y and are almost never selected. For each trial, the mean squared error (MSE) is computed as a function of the number of variables in the model. Figure 4b shows the average of the errors over the 100 simulations. On average, the model selected by minimizing the MSE includes four groups but the model with five groups also has an error close to the minimum. The simulation with the logistic link gives similar results. However the fifth group G(3) is more frequently selected (Fig. 5c). The minimization of the MSE leads to the selection of five groups as shown in Figure 5b. Note that this random forest and grouped variable importance-based approach performs well even with a nonlinear link.

In both experiments, the grouped variable importances obtained at the first step of the algorithm are ranked in the same order as the θ j . Indeed the impact of the correlation between predictors is almost null in both cases because of the orthogonality of the wavelet bases. In this context, the backward Algorithm 1 does not provide additional information compared to the "non-recursive" strategy (see the discussion following Algorithm 1). 4) G( 5) G( 6) G( 7 Experiment 3: Selection of functional variables in the presence of strong correlation This simulation illustrates the interest of Algorithm 1 when selecting functional variables in the presence of correlation. First, we simulate n = 1000 i.i.d. realizations of p = 10 functional variables X 1 , . . . , X p according to the general simulation design detailed before. For all i ∈ {1, . . . , n}, let Z 1 i , . . . , Z p i be latent variables drawn from a standard Gaussian distribution. The outcome variable Y is defined as:

Gζ G(0) G(1) G(2) G(3) G(
Y i = 3.5 Z 1 i + 3 Z 2 i + 2.5 Z 3 i + 2.5 Z 4 i .
Then for u ∈ {1, . . . , p}, the wavelet coefficients are simulated according to ( 11) and ( 10) with a linear link:

ξ u ijk = ω u jk + Z u i + ση u ijk if j ≤ j , k ∈ {0, . . . , 2 j -1} and ζ u i = ω u 0 + Z u i + ση u iζ
, with σ = 0.1, N = 2 9 and thus J = 9. We take j = 3 in order to make the functional variables smooth enough. Among the 10 variables X u , the first four variables have decreasing predictive power whereas the others are independent of Y . Next, we add q = 10 i.i.d. variables X 1,1 , . . . , X 1,q which are strongly correlated with X 1 : for any v ∈ {1, . . . , q}, any i ∈ {1, . . . , n},

ζ 1,v i = ζ 1 i + ση 1,v iζ and ξ 1,v ijk = ξ 1 ijk + ση 1,v ijk if j ≤ j , k ∈ {0, . . . , 2 j -1} 0 otherwise,
where σ = 0.05. The η 1,v ijk and the η 1,v iζ are i.i.d. standard Gaussian random variables. The discrete processes X 1,v i are obtained using the inverse wavelet transform. In the same way, we add q = 10 i.i.d. variables X 2,1 , . . . , X 2,q which are strongly correlated with X 2 . To sum up, the vector of predictors is composed of X 1 , X 1,1 , . . . , X 1,q , X 2 , X 2,1 , . . . , X 2,q , X 3 , X 4 , . . . , X p .

The aim is to identify the most relevant functional variables for the prediction of Y . For each functional variable, we regroup all wavelet coefficients and apply Algorithm 1. This process is repeated 100 times. Boxplots of the group permutation importances at the first step of the algorithm, over the 100 trials and for each functional variable, are given in Figure 6a. We see that the importances of the variables X 1 and X 2 and their noisy replicates are much lower than those of variables X 3 and X 4 . This is due to the strong correlations between the two first variables and their noisy replicates. Indeed, the importance measure decreases when the correlation or the number of correlated variables increase. Note that the importances of X 1 and X 2 are slightly lower than those of their noisy replicates. This can be explained by the fact that the correlation between X 1 and its noisy replicates is higher than, for instance, the correlation of X 1,1 with X 1 , X 1,2 , . . . , X 1,q .

Figure 6b gives a comparison of the performance of Algorithm 1 with the "non-recursive" strategy (NRFE). Algorithm 1 clearly shows better prediction performances. In particular, it reaches a minimum error faster than the NRFE: only five variables for Algorithm 1 whereas NRFE needs about twelve. The RFE procedure is more efficient than the NRFE when predictors are highly correlated.

Additional information is displayed in Figure 6c. The selection frequencies using Algorithm 1 show that the variables X 3 and X 4 are always selected. Indeed, these two variables have predictive power and they are not correlated with the other predictors. Note that X 1 and X 2 are selected less often than their replicates, even though they are more correlated with Y . This also comes from the fact that the correlation between X 1 and its replicates is higher than the correlation of X 1,1 with X 1 , X 1,2 , . . . , X 1,q . We observe that X 1 and X 2 are eliminated in the first steps of the backward procedure, but this has no consequence on the prediction performance of Algorithm 1. These results motivate the use of this algorithm in practice: it reduces the effect of the correlation between predictors on the selection process and provides better prediction performances. including aircraft speed, accelerations, the heading, position, and warning signals. Each flight provides a multivariate time series corresponding to this family of functional variables.

We focus here on the risk of long landing. A data sequence for N = 512 seconds before touchdown is observed for predicting the landing distance. The evaluation of the risk of long landings is crucial for safety managers to avoid runway excursions and more generally to keep a high level of safety. One answer to this problem is to select the flight parameters that best explain the risk of long landings. Therefore, we attempt to find a sparse model showing good predictive performances. In the future, flight data analysis could be used for pilot training or development of new flight procedures during approach.

Following aviation experts, 23 variables are preselected, see Table 1. A sample of 1868 flights from the same airport and the same company is considered. The functional variables are projected on a Daubechies wavelet basis with four vanishing moments using the discrete wavelet transform algorithm as in Section 3. The choice of the wavelet basis is determined by the nature of the flight data, which here contains information on both the time and frequency scales.

Abbreviation

Flight

Preliminary dimension reduction

The design matrix formed by the wavelet coefficients for all of the chosen flight parameters has dimension 23 × 512 = 11 776. Selecting the variables directly from the whole set of coefficients is computationally prohibitive, so we first need to significantly reduce the dimension. The naive method that shrinks the n curves independently according to [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] and then brings the non-zero coefficients together in a second step would lead us to consider a large block of coefficients with many zero values. This is not relevant in our context, so we propose an alternative method which consists of shrinking the wavelet coefficients of the n curves simultaneously. More precisely, this method is adapted from [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] for the particular context of n independent (but not necessary identically distributed) discrete random processes. Shrinkage is performed on the norm of the n-dimensional vector containing the wavelet coefficients. The complete method is described in Appendix C.

Selection of flight parameters

We obtain a selection of the functional parameters by grouping together the wavelet coefficients of each flight parameter and applying Algorithm 1 with these groups. At each iteration, we randomly split the dataset into a training set containing 90 % of the data and a validation set containing the remaining 10 %. In the backward algorithm, the grouped variable importance is computed on the training set and the validation set is only used to compute the MSE errors. The selection procedure is repeated 100 times to reduce the variability of the selection. The final model is chosen by minimizing the averaged prediction error. Figure 7a represents the boxplots of the grouped variable importance values computed on the 100 runs of the selection algorithm. According to this ranking, five variables are found to be significantly relevant. Looking at the averaged MSE estimate on Figure 7b, we see that the average number of selected variables is ten, but taking only five is sufficient to get a risk close to the minimum.

Figure 7c gives additional information, displaying the proportion of times each flight parameter is selected. First, it confirms the previous remarks: five variables are always selected by the algorithm and the first ten are selected more than 60 times over the 100 runs. Second, it shows that flight parameters related to aircraft trajectory during approach are among the most relevant for predicting long landings. Indeed, the elevators (ELEVL, ELEVR) are used by the pilots to control the pitch of the aircraft. These have an effect on the angle of attack (AOA) and consequently on the landing. The variable GLIDE.DEVC is the glide slope deviation, i.e., the deviation between the aircraft trajectory and a glide path of approximatively three degrees above horizontal. Other significant variables are the gross weight (GW.KG) which has an effect on the deceleration efficiency, the airspeed (CASC) and the engine rating (N11C, N12C).

It should be noted that the ranking based on selection frequency is close to the direct ranking given by the importance measures when all variables are included in the model. This suggests that for this data, correlation between predictors is not strong enough to influence their permutation importance measures. Moreover, if we regroup several variables such as for example flight parameters N11C and N12C, N21C and N22C, and ELEVR and ELEVL into three new variables N1, N2 and ELEV, Figure 8 shows that the ranking remains unchanged.

Selection of wavelet levels

We now determine, for a given flight parameter, which wavelet levels are the most able to predict the risk of long landing. We perform selection of wavelet levels independently for the gross weight (GW.KG) and angle of attack (AOA), which are among the most commonly selected flight parameters. Figures 9a and9c show that the average number of selected levels for GW.KG is less than for AOA. Indeed, the selection frequencies in Figure 9b indicate that for GW.KG, the first approximation levels are selected at each run (groups ζ, G(0) and G(1)) whereas the last levels are selected in less than 40 of the 100 trials. The situation is quite different for AOA: all levels are selected more than 50 times. The predictive power of this functional variable is contained in both the high levels of approximation and in the details of the wavelet decomposition (Fig. 9d).

Detection of important time intervals

The importance of time interval is now computed for two of the most important flight parameters, the altitude (ALT.STDC) and angle of attack (AOA). Figure 10 displays the averaged grouped importance G(t) evaluated on 50 equally-spaced time points. This number is arbitrarily chosen in order to find a tradeoff between computational time and accuracy. The time t = 0 refers to touchdown of the aircraft.

Two intervals are detected with high predictive power for the altitude. These results are consistent with the view of aviation safety experts. Indeed, during the interval [-460, -400], the aircraft has to level off for stabilizing before the final approach. An overly high altitude at this moment can induce a long landing. During the interval [-60, 0], the final seconds before touchdown, an overly high altitude can also induce a long landing.

The interval detected for the angle of attack is [-200, -100]. This make sense because according to flight procedure pilots have to reduce the vertical speed a few seconds before touchdown (the flare). Consequently, they must increase the angle of attack in order to keep sufficient lift.

G W .K G A O A E L E V L E L E V R G L ID E .D E V C A L T .S T D C C A S C .G S C C A S C N 1 2 C N 1 1 C S A T P IT C H N 2 1 C V A P P IV V N 2 2 C D R IF T H E A D .M A G V R T G R U D D R O L L F
E L E V L G W .K G A O A E L E V R G L ID E .D E V C A L T .S T D C N 1 1 C N 1 2 C C A S C .G S C C A S C S A T N 2 1 C V A P P P IT C H N 2 2 C D R IF T IV V H E A D .M A G V R T G R O L L R U D D F L A P C P IT C H .R A T E 0.0 0.2 0.4 0.6 0.8 1.0 (c)
G W .K G A O A E L E V L E L E V R G L ID E .D E V C A L T .S T D C C A S C .G S C C A S C N 1 2 C N 1 1 C S A T P IT C H N 2 1 C V A P P IV V N 2 2 C D R IF T H E A D .M A G V R T G R U D D R O L
G W .K G A O A E L E V G L ID E .D E V C A L T .S T D C C A S C .G S C C A S C N 1 S A T P IT C H V A P P N 2 IV V D R IF T H E A D .M A G V R T G R U D D R

Conclusion

We have considered the selection of grouped variables using random forests and proposed a new permutation-based importance measure for groups of variables. Our theoretical analysis provided exact decompositions of the grouped importance measure into a sum of the individual importances for specific models such as additive regression models. A simulation study highlighted the fact that in general the importance of a group does not reduce to the sum of the individual importances. Since the idea of variable importance is not restricted to the random forest algorithm, it would be of interest to extend our method to other learning algorithms such as different bagging methods, neural networks and SVMs.

We have also introduced a rescaled version of the grouped variable importance with a heuristic choice for the normalization factor. The question of choosing a good value for the normalization factor is an open question that could perhaps be examined with further mathematical work.

The second contribution of the article was an original method for selecting functional variables using the grouped variable importance measure and a projection-based approach. The various functional variables were projected onto a wavelet basis and the corresponding coefficients grouped in different ways. For instance, one can choose to regroup all coefficients of a given functional variable. The selection of these groups involved selection of the functional variables. Various other groupings were proposed for wavelet decompositions. A backward algorithm adapted from recursive feature elimination in the context of random forests was used as a selection algorithm. The selection process was guided by the grouped importance measure in order to select the most relevant group of variables for predicting the variable of interest.

Lastly, an extensive simulation study was performed to illustrate application of the grouped importance measure for functional data analysis, and encouraging results were then obtained for a real life problem in aviation safety. Of course, the grouped variable importance measure may be used for any other application as soon as the features can be partitioned into known disjoint pieces. Potential applications can be found in the extensive literature about the group lasso method (Yuan and Lin, 2006b), see for instance [START_REF] Meier | The group lasso for logistic regression[END_REF]; [START_REF] Ma | Supervised group lasso with applications to microarray data analysis[END_REF] for applications in Genetics and [START_REF] Chatterjee | Sparse group lasso: Consistency and climate applications[END_REF] for applications in Climatology. Let 0 p and I p denote the null vector and identity matrix of R p . Let 1 p be the vector in R p with all coordinates equal to one and let 0 p,q denote the null matrix of dimension p × q.

Experiment 1: linear link function.

We simulate X and Y from a multivariate Gaussian distribution. More precisely, we simulate samples from the joint distribution

X Y =   W Z Y   ∼ N 2p+1   0 2p+1 ,   C w 0 p,p τ 0 p,p I p 0 p τ 0 p 1     ,
where τ is the vector of the covariances between W and Y . In this context, the conditional distribution of Y over X is normal and the conditional mean f is a linear function: f (x) = p+q j=1 α j x j with α = (α 1 , . . . , α p , 0, . . . , 0) a sequence of deterministic coefficients (see for instance [START_REF] Rao | Linear statistical inference and its applications[END_REF], p. 522.

• Experiment 1a: independent predictors. We take τ = 0.9 1 p and C w = I p . All variables of W are independent and correlated with Y .

• Experiment 1b: correlated predictors. We take τ = 0.9 1 p and C w = (1 -0.9)I p + 0.91 p 1 p . The variables of W are correlated, and also correlated with Y .

• Experiment 1c: independent predictors, sparse case. We take τ = (0.9, 0, . . . , 0) and C w = I p .

Only the first variable in the group W is correlated with Y .

• Experiment 1d: correlated predictors, sparse case. We take τ = (0.9, 0, . . . , 0) and C w =

(1 -0.9)I p + 0.91 p 1 p . The variables of W are correlated. Only the first variable in the group W is correlated with Y .

Experiment 2: additive link function.

We simulate X from a multivariate Gaussian distribution:

W = W Z ∼ N 2p 0 2p , C w 0 p,p 0 p,p I p ,
and the conditional distribution of Y is

(Y |X) ∼ N   p j=1 f j (X j ), 1   ,
where f j (x) = sin(2x) + j for j < p/2 and f j (x) = cos(2x) + j for j ≥ p/2.

• Experiment 2a: independent predictors. We take C w = I p . All variables of W are independent and correlated with Y .

• Experiment 2b: correlated predictors. We take C w = (1 -0.9)I p + 0.91 p 1 p . The variables of W are correlated and also correlated with Y .

Experiment 3: link function with interactions.

We simulate X from a multivariate Gaussian distribution: and the conditional distribution of Y is

W = W Z ∼ N 2p (0 2p I 2p ) , 20
(Y |X) ∼ N   p j=1 X j + X p X 1 + p-1 j=1 X j X j+1 , 1   .

Results

Experiments 1a-b and 2a-b illustrate the results of Corollary 1 (see Figures B.11 and B.12). Indeed, the regression function in both cases satisfies the additive property (3) for these experiments. In Experiments 1a and 2a, the variables of the group W are independent and the grouped variable importance is nothing more than the sum of the individual importances in this case. In Experiments 1b and 2b, the variables of group W are positively correlated. In these situations, the grouped variable importance is larger than the sum of the individual importances, which agrees with Equation 4. Note that the grouped variable importance increases with p, which is natural because the amount of information for predicting Y increases with the group size in these models. On the other hand, it was shown in [START_REF] Gregorutti | Correlation and variable importance in random forests[END_REF] that individual importances decrease with correlation between the predictors. Indeed we observe that the sum of the individual importances decreases with p in the correlated cases.

The regression function in Experiment 3 does not satisfy the additive form (3). Although the variables in the group are independent, the grouped variable importance is not equal to the sum of the individual importances (Figure B.13). In a general setting therefore, it appears that these two quantities differ.

We now comment the results of the sparse Experiments 1c and 1d (Figure B.14). It is clear that I(W) = I(X 1 ) = j=1...p I(X j ) for these two experiments (see Proposition 1 for instance). Regarding Experiment 1c, the boxplots of the estimated values Î(W) and j=1...p Î(X j ) agree with this equality. On the other hand, in Experiment 1d, we observe that j=1...p Î(X j ) is significantly higher than Î(W). Indeed, it has been noticed by [START_REF] Nicodemus | The behaviour of random forest permutationbased variable importance measures under predictor correlation[END_REF] that the individual importances of predictors that are not associated with the outcome tend to be overestimated by random forests when there correlation between predictors. In contrast, the estimator Î(W) seems to correctly estimate I(W) even for large p. Indeed, for both experiments 1c and 1d, the importance Î(W) is unchanged when the size of the group p varies from 2 to 16. Note that for variable selection, we may prefer to consider the rescaled importance Înor so as to select with priority small groups of variables.

Appendix C. Curve dimension reduction with wavelets

The analysis of flight data in Section 4 required, for computational reasons, to preliminarily reduce the dimension of the wavelet decomposition of the flight parameters. We require to adapt the well-known wavelet shrinkage method to the context of independent random processes. Using the notation of Section 3, we first recall the hard-thresholding estimator introduced by [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] in the case of one random signal. This approach is then extended to deal with n independent random signals.

Signal denoising via wavelet shrinkage

The problem of signal denoising can be summarized as follows. Suppose that we observe N noisy samples X(t 1 ), . . . , X(t N ) of a deterministic function s (the signal): X(t ) = s(t ) + σε , = 1, . . . , N, (C.1)

where the ε are independent standard Gaussian random variables. We assume that s belongs to L 2 ([0, 1]).

The goal is to recover the underlying function s from the noisy data {X(t ), ∈ {1, . . . , N }} with small error.

Using the discrete wavelet transform, this model can be rewritten in the wavelet domain as ξ jk = ω jk + ση jk , ∀j ∈ {0, . . . , J -1}, ∀k ∈ {0, . . . , 2 j -1}, and the scaling domain as ζ = ω 0 + ση 0 , where ξ jk and ζ are the empirical wavelet and scaling coefficients of X(t ) as in Equation ( 8). The random variables η jk and η 0 are i.i.d. random variables from the distribution N 1 (0, 1).

A natural approach for estimating ω jk is to shrink the coefficients ξ jk to zero. An estimator of s in this context has the form ŝ(t ) = ω0 φ(t ) + This method is referred to as the hard-thresholding estimator in the literature. [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] propose the universal threshold δ N = σ 2 log(N ). In addition, the standard deviation σ can be estimated by the median absolute deviation (MAD) estimate of the wavelet coefficients at the finest levels, i.e., σ = Med(|ξ jk -Med(ξ jk )| : j = J -1, k = 0, . . . , 2 J-1 -1) 0.6745 , where the normalisation factor 0.6745 comes from the normality assumption in (C.1). This estimator is known to be a robust and consistent estimator of σ. The underlying idea is that the variability in the wavelet coefficients is essentially concentrated at the finest level.

Consistent wavelet thresholding for independent random signals • Identically distributed case. We start by assuming that the observations come from the same distribution: for any i ∈ {1, . . . , n}, X i (t ) = s(t ) + σε i, , = 1, . . . , N, (C.3)

where the ε i, are independent standard Gaussian random variables. The wavelet coefficients of s can be easily deduced from the mean signal X := 1 • Non identically distributed case. In many real life situations, assuming that the n signals are identically distributed is not a realistic assumption. For the study presented in Section 4 for instance, the flight parameters have no reason to follow the same distribution in safe and unsafe conditions. We propose a generalization of the model (C.3) by introducing a latent random variable Z taking its values in a set Z.

Roughly speaking, the variable Z represents all the phenomena that have an effect on the mean signal.

Conditionally on Z i = z i , the distribution of the process X i is now defined, for any i ∈ {1, . . . , n}, by X i (t ) = s(t , z i ) + σε i, , = 1, . . . , N, (C.4)

where the ε i, are independent standard Gaussian random variables. This regression model allows us to consider various situations of interest arising in functional data analysis. In supervised settings where a variable Y has to be predicted using X, one reasonable model is to take Z = Y . We now propose a hard-thresholding method which simultaneously shrinks the wavelet decomposition of the n signals.

Let • n denote the 2 -norm in R n : u n := n i=1 u 2 i for any u ∈ R n . Let ξ jk be the vector (ξ 1jk , . . . , ξ ijk , . . . , ξ njk ) where ξ ijk is the coefficient of level (j, k) in the wavelet decomposition of

Figure 1 Figure 1 :

 11 Figure1displays the matrix giving the correspondence between a time location and the associated wavelet functions for a Daubechies wavelet basis with two vanishing moments. In a similar way but for a time

Figure 2 :

 2 Figure 2: Experiment 1 -Example of ten processes drawn from the protocol used for Simulation 2.

Figure 3 :

 3 Figure 3: Experiment 1 -Averaged time importances and first and third quartiles over 100 iterations. The time interval T is located between the two vertical lines.

Figure 4 :

 4 Figure 4: Experiment 2, linear links -selection of the wavelet levels. From left to right: (a) Boxplots of the grouped variable importances, (b) MSE error versus the number of groups and (c) Selection frequencies.

Figure 5 :

 5 Figure 5: Experiment 2, logistic links -selection of the wavelet levels. From left to right: (a) Boxplots of the grouped variable importances, (b) MSE error versus the number of groups and (c) Selection frequencies.

Figure 7 :

 7 Figure 7: Application to long landings -from left to right: (a) Boxplots of the grouped variable importance, (b) MSE error versus the number of groups and (c) selection frequencies.

Figure 8 :

 8 Figure 8: Application to long landings -grouped variable importance measured before and after grouping the correlated flight parameters N11C and N12C, N21C and N22C, ELEVR and ELEVL, as N1, N2 and ELEV.

Figure 9 :

 9 Figure 9: Application to long landings -selection of wavelet levels for the gross weight and angle of attack.

  Figure B.11: Boxplots of the importance measures for Experiments 1a and 1b. The number of variables in W varies from 1 to 16. For Experiment 1a, the sum of the individual importance@s and Î(W) overlap.

  Figure B.12: Boxplots of the importance measures for Experiments 2a and 2b. The number of variables in W varies from 1 to 16. For Experiment 2a, the sum of the individual importances and Î(W) overlap.

  Figure B.13: Boxplots of the importance measures for Experiment 3. The number of variables in W varies from 1 to 16.

  jk if |ξ jk | > δ N 0 otherwise.

  . . . , N, where the ε are independent standard Gaussian random variables. By applying the hard-thresholding rule to this signal, we obtain the following estimation of the wavelet parameters of s: ω0 = ζ and ωjk = ξjk if | ξjk | > δN 0 otherwise, where ξjk is the wavelet coefficient of level (j, k) of X. Here the threshold is δN = σ √ n 2 log(N ).

Table 1 :

 1 List of preselected flight parameters.

	parameter

Grouped importanceX 1 X 1 ,1 X 1 ,2 X 1 ,3 X 1 ,4 X 1 ,5 X 1 ,6 X 1 ,7 X 1 ,8 X 1 ,9 X 1 ,1 0 X 2 X 2 ,1 X 2 ,2 X 2 ,3 X 2 ,4 X 2 ,5 X 2 ,6 X 2 ,7 X 2 ,8 X 2 ,9 X 2 ,1 0 X 3 X 4 X 5 X 6 X 7 X 8 X 9 X 1 0

Selection frequenciesX 4 X 3 X 1 ,8 X 1 ,5 X 1 ,9 X 1 ,1 X 1 ,3 X 1 ,1 0 X 1 ,7 X 1 ,2 X 2 ,1 X 1 ,6 X 1 ,4 X 2 ,5 X 2 ,3 X 2 ,1 0 X 2 ,8 X 2 ,6 X 2 ,4 X 2 ,9 X 2 ,7 X 2 ,2 X 1

A case study: variable selection for aviation safety

In this section, we study a real problem coming from aviation safety. Airlines collect large amounts of information during flights using flight data recorders. For several years now, airlines are required to use these data for flight safety purposes. A large number of flight parameters (up to 1000) are recorded each second, Appendix A. Proof of Proposition 1

The vector X J and the vector X (J) are defined as in Section 2. We have

as X J and X J are independent and identically distributed.

Appendix B. Additional experiments on the grouped variable importance

In this section, we investigate the properties of the permutation importance measure of groups of variables with numerical experiments, in addition to the theoretical results given before. In particular, we compare this quantity with the sum of individual importances in various models. We also study how this quantity behaves in "sparse situations" where only a small number of variables in the group are relevant for predicting the outcome.

The general framework of the experiments is the following. For a fixed p ≥ 1, define X := (W , Z ) where W and Z are random vectors both of length p. Some of the components of W are correlated with Y whereas those in Z are all independent of Y . Let C w be the variance-covariance matrix of W. By incorporating the group Z in the model, we present a realistic framework where not all the X j have a link with Y . For each experiment, we simulate n = 1000 samples of Y and X and compute the importance I(W) of the group W, the rescaled grouped variable importance I nor (W) and the sum of the individual the signal X i . For any z ∈ Z, let ω jk (z) be the wavelet coefficient of level (j, k) of s(•, z), and ω jk := (ω jk (Z 1 ), . . . , ω jk (Z n )) . We define the common wavelet support of s by

If (j, k) ∈ L, then ω jk = (0, . . . , 0) almost surely and ξ jk 2 n has a centered chi-square distribution with n degrees of freedom. Otherwise, ω jk can be non null and in this case ξ jk 2 n has the distribution of a sum of n independent uncentered chi-square distributions. We can thus propose a thresholding rule for the statistic ξ jk n . For any j ∈ {0, . . . , J -1} and any k ∈ {0, . . . ,

where the threshold δ N,n depends on N, n and σ. Proving adaptive results in the spirit of [START_REF] Donoho | Wavelet shrinkage: asymptopia[END_REF] for this method is beyond the scope of the paper. However, an elementary consistent result can be proved. We would like ωjk to be a zero vector with high probability when (j, k) ∈ L. For some x ≥ 0, take

where we have used a deviation bound for central chi-square distributions from [START_REF] Laurent | Adaptive estimation of a quadratic functional by model selection[END_REF]Massart (2000, p. 1325). If the signal is exactly zero, it can be recovered with high probability by taking x log(N ). In particular, if we choose x = 2 log(N ), the threshold is δ 2 N,n = σ 2 (4 log(N ) + 2 2n log(N ) + n) and the convergence rate in (C.6) is of order O( 1 N ). In practice, σ can be estimated by a MAD estimator computed on the coefficients of the highest level of all n wavelet decompositions. Next, x and δ N,n can be chosen such that (C.6) is lower than a given probability q. Letting N e -x = q, we obtain the threshold

.

Assuming that ω j,k (Z) = 0 almost surely for some level (j, k) is a strong assumption that is difficult to meet in practice. It can be hoped that this method still works if the wavelet support of s(•, z) does not vary too much with z. In particular, it may be applied if there exists a common set S of indices (j, k) such that, for any z, the projection of s(•, z) on Vect(ψ jk | (j, k) ∈ S) is not too far from s(•, z) for the L 2 norm.


Introduction

Understanding and modeling cement based materials has always been a great challenge towards numerous issues, such as durability assessments or effective properties predictions. In this paper an original attempt to model cement paste microstructure is presented, along with an original numerical framework.

First, an original morphological tool is introduced, based upon level-set methods applied to Correlated Random Fields (hereafter CRFs) [START_REF] Roubin | [END_REF]). Unlike most of the existing models, such tools allow to yield "random" morphologies, with targeted volume / surface area fractions. Furthermore, this method allows the control of the entire morphology through a reduced set of parameters, and their evolution is driven by a unique one. In opposition to "particular" models (e.g. HYMOSTRUC, van Breugel (1991)) where a set of three parameters are needed to describe every inclusion (in 2D), or "voxelized" models (CEMHYD, [START_REF] Bentz | Digital-image-Based computer modeling of cement-based materials[END_REF]) where local rules governing the evolution through time ; here a unique threshold will control the entire morphology. This framework is applied at the cement paste scale, in 2D, in order to yield morphologies of hydrating cement paste.

The Lattice Boltzmann Method (LBM) provides a simple way to simulate transport phenomena in heterogeneous media such as cement paste (e.g the 3D modelling of [START_REF] Zhang | Lattice Boltzmann simulation of the ionic diffusivity of cement paste[END_REF]. A regular mesh defined once only allows for a simple implementation and its explicit formulation makes it attractive for the simulation of evolving geometries at different time steps, which is the final aim of this study. The parallelization capacity of the LBM also implies a competitive computational cost time compare to traditional numerical methods. This aspect is of particular interest for variability studies: given the statistical CRF-based reconstruction of cementitious material, a high number of simulations may be required for variability studies, such as presented in the final application of this work, which is the extension of the proof of concept provided in Walther et al. (2014).

Due to a higher percolation effects through pore connectivity in three dimensions, the two-dimensional calculation of effective diffusivity is indeed less relevant in the case open-cell porosity, however the originality of this work -compared to the extensive approach of [START_REF] Zhang | Lattice Boltzmann simulation of the ionic diffusivity of cement paste[END_REF] and [START_REF] Zhang | Modeling of ionic diffusivity in non-saturated cement-based materials using lattice Boltzmann method[END_REF] -relies on the stability analysis of the LBM as well as in the statistical determination of the equivalent diffusivity, based on numerous calculations of independently generated morphologies as in Roubin et al. (2014). Furthermore, this work is to be considered as a validation of the in-house developed LBM code and will be extended to three dimensions in the future.

Morphology creation and description

In this first section a morphological framework allowing to easily yield complex morphologies is presented, after [START_REF] Roubin | [END_REF] and Bogdan et al. (2014). A simplified hydration model is then introduced, and when combined with the morphological yielding process, it allows to generate hydrating cement paste morphologies.

Correlated Random Fields and Level-set methods

The presented morphological framework is based upon Gaussian CRFs thresholding. Those particular random fields are chosen because they can be described with very few parameters: a correlation length L c , a variance σ and an underlying correlation function. They represent respectively the average width and height of the peaks / valley on Figure 1(a). Lastly the mean value of the RF might be added, but in the following zero mean Gaussian CRFs are used.

Within this work, those CRFs are considered as random functions: γ(x, ω) = M ∈ R 2 → R, where M is a two dimensions Euclidean space, x is the space variable, and ω the random variable. Thus, the CRFs are defined over a square of size T , with values in R. An excursion A u of such a function is defined by "the set of points where the values of γ are above a given threshold u", as presented in Equation [1].

A u ≡ A u (γ, M ) = {x ∈ M |γ(x, ω) ≥ u} (1)
An example is shown on Figure 1. On the left side, a realization γ of a CRF over a square is presented, and on the right side the thresholding with the yielded excursion beneath. The key point in using this framework lies in the work of Adler (2008). It establishes a probabilistic link between the CRFs' quantities (L c and σ) and the threshold value (u) on one side, and the geometrical and topological quantities of the excursion sets on the other side (System [2], in the general 3D case, V being the volume fraction, S the surface area, and EC the Euler's characteristic that can be seen as a measure of the connectivity). It is thus possible, for a given threshold, to access statistical information about the surface area fraction, the average size of the inclusions, and to a measure of the connectivity for a realization of such morphologies. The predictions are never perfect since the System [2] is only made in terms of expectancies. This link makes this method very powerful, since it is possible to generate morphologies with a chosen surface area, or a given inclusion size, and then make them evolve through a unique parameter, the threshold. For more details about the probabilistic link and the explicit formulae (functions f , g, and h in System [2]), the reader is referred to Adler (2008) and [START_REF] Roubin | [END_REF].

E {V } = f (σ, L c ) (2a) E {S} = g(σ, L c , u) (2b) E {EC} = h(σ, L c , u) (2c)
In the following this process is used to yield cement paste morphologies, and thus numerous phases has to be represented. The same way a threshold is able to separate a CRF in two parts ; two thresholds yield three phase morphologies, with a new phase "concentric" to the initial inclusions (see Figure 2(c) and 2(d)). In order to describe cement paste in hydrated states as well as during hydration, a model is needed to describe their evolution. In a first approximation, the Powers model is implemented [START_REF] Powers | Studies of the physical properties of hardened Portland cement paste[END_REF] for its simplicity (few phases considered in the reactants as well as for the new hydrates, and a macroscopic description of the phases volume fraction evolutions).

Usually, this model considers five volume fractions: anhydrous cement, hydration products, gel water, free water and chemical shrinkage. It has been chosen to use only three for this first application, taking into account only anhydrous cement, hydration products (including gel water) and free water. A more detailed description of the cement paste is possible, but would require a more advanced hydration model (e.g. Jennings & Tennis (1994) with each hydrate treated separately). The chemical shrinkage is also neglected in this first implementation, and its volume fraction is allocated to porosity (i.e. to the free water phase).

Thus, the overall description is given by System [3]:

p = w/c w/c + ρ w /ρ c (3a) V anh =(1 -p)(1 -ξ) (3b) V h =2.12(1 -p)ξ (3c) V w =1 -V anh -V hyd ( 3d
)
where p is the initial porosity, ξ the hydration degree, and V anh , V h , V w respectively the volume fractions of anhydrous cement, hydration products and water. An example is given Figure 3.

Remark: The simplicity of the hydration description through a unique parameter (the hydration degree ξ) is reminded, as its morphological description through a set of two threshold governing the whole evolution of the microstructures. From a practical point of view, we will consider the following:

-The initial porosity p give us the initial surface area (volume fraction) of the cement grains in the fresh cement paste, -the evolution of the three phases, and thus two threshold, is governed by a unique parameter.

-A maximal hydration degree is considered -ξ max (not developed here) -to undermine the strong hydration model simplifications.

For more information about level-set methods applied to hydration modeling, and the extension to three dimensional simulations please refer to Bogdan et al. (2014). It is shown that the particle size distribution of the anhydrous cement is taken into account (using unions of excursion sets with various correlation lengths, somewhat visible on Figure 3, with various inclusion sizes), as well as more developed hydration

Calculation of the effective diffusivity of hydrated cement paste

In this section the interest of using the LBM is presented, in particular for effective diffusivity calculation in heterogeneous media. The aspects of the associated numerical stability are then reviewed. Eventually, the validation of the code versus an analytical solution and an example of result is presented.

Diffusive transport in porous media

Classically, the molecular diffusivity is given by Fick's law as per Equation [4], where D j is the diffusivity of species j and c its concentration.

∂c ∂t = -D j ∇c (4)
This formulation of diffusivity is valid locally on every phase j independently, which means it does not take into account the ionic interactions that may actually occur depending on the species diffusing in the considered materials (hence the effective diffusivity prediction, in contrast to the apparent one).

Regarding heterogeneous media, the homogenized effective diffusivity D eq is calculated by applying a constant macroscopic concentration gradient between two opposite sides (based on Hill's criterion [START_REF] Hill | [END_REF] -and the energy-dissipation equiv-alence between scales), as in [START_REF] Schulte | Bounding of effective thermal conductivities of multiscale materials by essential and natural boundary conditions[END_REF]. The macroscopic flux obtained by integration of the species flux on the boundaries [START_REF] Pouya | Définition de la perméabilité équivalente des massifs fracturés par des méthodes d'homogénéisation[END_REF]) is then related to the effective diffusion coefficient.

The Lattice Boltzmann Method

The LBM is a mesoscopic approach for partial derivative equations problems. The evolution of a fluid or a solid is deduced from the consideration of a population of particles interacting with well-defined rules of collision. Unlike the well-known Navier-Stokes approach that uses the equations describing the macroscopic behaviour of a fluid and then discretizes them to the microscopic level, the LBM is an "ascendant" method, where the particles' agitation defined by the Maxwell-Boltzmann distribution of statistical mechanics allows for a consistent extrapolation towards the macroscopic behaviour. The link between micro and macroscopic scales is reached thanks to a Chapman-Enskog expansion (Wolf-Gladrow, 2000). Practically, the LBM consists in considering a population of particles streaming along given directions in the space, systematically on the grid points. The particles' displacement is equal to unity in every direction, hence the time and space scales are strongly linked. At each time step, the particles' distributions stream orthogonally and diagonally to their neighbouring nodes and exchange their energy by collision. Figure 4 gives an example of distribution propagation for the D 2 Q 9 model, where subscript "2" represents the dimension of space, "9" being the number of distributions streaming on the grid (8 possible space directions, plus the probability of staying on the same point).

Boltzmann's equation ruling the particles' behaviour is given under following form, f being the density of probability of the particles, c the speed of propagation and ω the collision operator:

∂f ∂t + c∇f = Ω(f ) (5)
After discretizing the space into as previously defined, with a "lattice unit per time step in direction k" as the mesoscopic speed of the particle on the grid, one obtains:

∂f k ∂t + c k ∂f k ∂dx = Ω k (6)
In the Bhatnagar-Gross-Krook approach (Nourgaliev et al. 2003), the collision operator is approximated as a relaxation of the distribution function towards the equilibrium distribution:

Ω k = - 1 τ (f k -f eq k ) (7)
The equilibrium distribution function f eq k of Equation [7] is computed through the Taylor expansion of the Maxwell-Boltzmann distribution of a stationary perfect gas at the equilibrium, which simplifies for the diffusion process merely as Equation [7]:

f eq k = ω k ρ (8)
The ρ factor in Equation [8] is simply the sum of the f k distributions at the point considered and represents the computed macroscopic value, e.g. the concentration or temperature. The ω k factor is a constant weighting relative to the propagation direction, such as k ω k = 1. Long displacements are less likely than short ones, hence the decreasing weightings of zero displacement, orthogonal displacement and diagonal displacement.

Setting ω = ∆t/τ out and writing the Boltzmann equation in its discrete form, one obtains the explicit equation of the space-time evolution of the functions:

f k (x + ∆x, t + ∆t) = f k (x, t) -ω(f k (x, t) -f eq k (x, t)) (9)
For stability reasons of the LBM, the ω factor has to remain strictly between 0 and 2. In the case the case of two-dimensional diffusion in a single given medium, the relationship between ω and the "LBM diffusivity" d is as follows:

d = c 2 s ( 1 ω - 1 2 ) (10)
In Equation [10], c s = 1/ √ 3 is a parameter called the sound speed in the LBM spacetime, but is actually only a scaling factor in the case of diffusion (Nourgaliev et al. 2003[START_REF] Sai | Non-negativity and stability analyses of Lattice-Boltzmann Method for advection-diffusion equation[END_REF]. In order to make a LBM simulation with a given number of time steps and a mesh size with the evolution of a real domain , one has to make their Fourier numbers correspond:

F o = d × t LBM N 2 = D × t real L 2 (11)
Changing the "LBM diffusivity" or refining the mesh will hence affect the number of time steps to be computed. For more details about the LBM and its boundary conditions, please refer to Wolf-Gladrow (2000). Regarding the simulation of diffusion in two materials, the best practice is to define a relaxation factor for one of them and then determine the value of the other one depending on the ratio of real diffusivities, as developed hereinafter: The phase "a" is considered as the reference phase (D a , d a )

and is defined by its relaxation factor ω a . The other phase must satisfy the relative ration between LBM and real worlds such as:

d a d b = D a D b (12)
The other relaxation factor is hence given by following relationship:

ω b = ( d b c 2 s + 1 2 ) -1 (13)
The major interest of the LBM is its simple numerical implementation even for complex morphologies, the controlled change of the relaxation factor ω depending on the type of media considered on the regular grid allowing for a correct representation of the physics, without having to change the mesh at the interface. As we have seen, two parameters ω a and ω b have to be controlled in order ensure a consistent result, which will be detailed in following section.

Numerical stability

Another aspect of this study is the behaviour of the LBM-BGK algorithm depending on the variation of the relaxation factor ω, which affects the precision, computational cost and stability of the method. As mentioned in the previous section, the relaxation factor depends on the LBM-diffusivity for a chosen problem. It also affects the computational cost, i. e. the number of time steps required to reach a given real time and the precision of the method: taking the explicit Finite Differences Method (FDM) as a reference for comparison, one can see on Figure 5 the evolution of the error ratio between LBM and FDM, as well as the time steps ratio. The horizontal dotted line is the unity ratio for both methods' time steps and gives a measure of the time steps ratio required for the LBM compared to the explicit FDM. From Figure 5, one can conclude that the relaxation factor shall be small enough to be competitive versus classical methods, ideally below the unity ratio, yet big enough to allow for a sufficient precision of the method. The experience shows that for heterogeneous me-Figure 6. Oscillating behavior for several diffusivity ratios depending on the relaxation factor dia with a diffusivity ratio reaching several orders of magnitude, some perturbations might appear at the interfaces between the phases. A numerical campaign of experiment allowed for the determination of the loci where no oscillations appear, as per Figure 6.

Numerical results

Theoretical aspects of the method being set, the validation of the model is developed in next section, followed by an application to diffusion in numerous morphologies of hydrated cement paste.

Validation

The effective diffusivity of the considered materials is assessed with the by summing all local contributions to yield the macroscopic equivalent flux. [START_REF] Mori | [END_REF] analytic approach of homogenization allows to verify the aptitude of the code to predict the equivalent diffusivity of a matrix with inclusions. Cylindrical, regular inclusions were put in the center of a matrix and the effective diffusivity was computed for different volume fractions of inclusion f v and different diffusivity ratios d a /d b , resulting in a good accordance with Mori & Tanaka's reference scheme for this type

Cement paste application: Results

In order to get a statistically representative description of hydrating cement paste, the level-set morphology reconstruction technique is used over 300 realizations of CRFs, for each degree of hydration. The hydrating cement paste microstructures (up to 100×100 µm, with a resolution of 0.33 µm/pixel) are composed of water-saturated porosity (with a diffusive coefficient set to D w = 2.10 -9 m 2 /s), hydrates (D h = 2.10 -12 m 2 /s), and anhydrous cement. The latter being considered as a bulk matrix, thus non-diffusive. An example of the evolution of the morphology is shown Figure 3 where one can see two degrees of hydration. The effective diffusivity of cement paste with a w/c ratio of 0.4 is assessed during its hydration following the presented hydration framework.

RVE size for diffusion in cement pastes:

Several calculation were first done to assess a Representative Volume Element in regards to diffusion. In several numerical simulations [START_REF] Zhang | Lattice Boltzmann simulation of the ionic diffusivity of cement paste[END_REF], [START_REF] Zhang | Modeling of ionic diffusivity in non-saturated cement-based materials using lattice Boltzmann method[END_REF] among others), the RVE is considered to be of size T =100 µm. Results for smaller morphologies are presented now, with intermediate microstructures of side size16.6 µm, 33 µm, 66 µm and finally 100 µm. The discretization stays constant (0.33 µm / pixel), and all presented results are averaged over 200 realizations. Size dependencie results are drawn Figure 8.

As the microstructures size increases, the value of the homogenized normalized effective diffusivity tends to stabilize. For the morphologies at ξ = 0.1, the effective diffusivity increases slightly with the microstrucure size, and the asymptotic value of diffusion coefficient is the one found in free water ; which was expected since nor the inclusions or hydrates are percolated. It is only for sizes above 60 µm that the framework capture all the relevant morphological information (e.g. for smaller sizes, one could be simulating only a anhydrous residue, or a bulk hydrate with no porosity. For ξ = ξ max = 0.73, the results interpretations is less evident with first a decrease in terms of effective diffusivity, and then a slight increase.

In view of those results, and in the spirit of following the numerous previous attempts to simulates CP pastes morphologies, a side size of 100 µm is considered in the following.

Hydrating cement paste diffusivity:

The results in terms of averaged effective diffusivity are plotted on Figure 10 (and an example of a concentration distribution is shown Figure 9), with a 95% confidence interval for the 300 morphologies per hydration state. At the very beginning, the effective diffusivity is the one of water, which is in agreement with a well percolated pore-phase. Then, as the hydration evolves, the effective diffusivity decreases, with a significant change of slope around ξ = 0.45. This is due to the percolation of the hydrate's phase. Finally, the pore space is filled by the hydrates, and the effective diffusivity is closer to the hydrates one.

In addition, as the presented results are averaged over 300 realizations, the 95% confidence bounds are represented. They show an interesting feature as these bounds are slightly tighter for early hydration states as well as for final ones. At those states, one of the phases is "very" percolated (i.e. water is percolated or the hydrates are closing the diffusion paths -Figure 10 also illustrates this behavior), and thus the macroscopic response depends first on this phase properties. On the other hand, in the vicinity of the change of slope those bounds are looser, as there is no predominant percolated phase, and morphological variability is higher.

Conclusions

This study presents a method for creating hydrating cement paste morphologies thanks to the level-set method.Through its explicit formulation and regular meshing, the LBM offers a simple implementation which is of interest for the calculation of the effective diffusivity of heterogeneous media having a strong variability of diffusive properties over time.

Future developments of this work include the extension of the LBM model to 3D, allowing for a proper comparison with experimental results. The tortuosity and pore connectivity in three dimensions are indeed very different than the two-dimensional ones. As they play a major role in transport properties, we expect the effective diffusivity to change significantly. This will allow to compare the 2D and 3D modeling approaches, and possibly to find a link to reduce even further the computational costs of the effective properties predicting framework (Qiao, 2012). The calculation of instantaneous diffusion properties of cement paste during hydration (id est changing morphology), as well as including the effects of other chemical species in presence shall also be included in the model.

Figure 1 .

 1 Figure 1. CRF realization, thresholding, and excursion

  (a) One CRF realization, and the threshold u (b) Initial morphology: cement grains and water (c) Same CRF realization, with two thresholds: v ≥ u, and w ≤ u (d) Cement hydration: cement grains, water, and hydrates

Figure 2 .

 2 Figure 2. From the RF creation to cement paste morphology yielding -grey: unreacted cement; blue: water; orange: hydration products

Figure 3 .

 3 Figure 3. Evolution of a selected morphology during hydration -w/c = 0.4, size T =100 µm (0.5 µm / pix) -Yellow: anhydrous cement, Green: water, Grey: hydration products

Figure 4 .

 4 Figure 4. D2Q9 Scheme -Arrows representing the possible directions of propagation of the distributions

Figure 5 .

 5 Figure 5. Comparison of the LBM versus the finite-differences method for the diffusion problem

  e

Figure 7 .

 7 Figure 7. Inclusion with f v = 15% used for the benchmark (left) and results of the code validation versus Mori & Tanaka's analytical approach (right)

Figure 8 .

 8 Figure 8. Normalized diffusivity (with D pores =1 m 2 s -1 and D hydrates /D pores = 10 -3 ) for various sizes of microstructure

Figure 9 .

 9 Figure 9. Concentration distribution in steady state under a unitary concentration gradient -w/c = 0.4, size T =100 µm (0.5 µm / pix)

Figure 10 .

 10 Figure 10. Simulated average effective diffusivity of standard cement paste (w/c = 0.4 -ξ ∈ [0; ξ max ]) with 95% confidence intervals