Christophe Crambes
email: crambes.christophe@univ-montp2.fr

Yousri Henchiri
email: yousri.henchiri@univ-montp2.fr

Regression imputation in the functional linear model with missing values in the response

Keywords: Functional linear model, Missing data, Missing at random, Principal components regression, Mean square error prediction. 2010 MSC: 62G20, 62J05, 62F12, 62P12

We are interested in functional linear regression when some observations of the real response are missing, while the functional covariate is completely observed. A complete case regression imputation method of missing data is presented, using functional principal component regression to estimate the functional coefficient of the model. We study the asymptotic behaviour of the error when the missing data are replaced by the regression imputed value, in a 'missing at random' framework. The completed database can be used to estimate the functional coefficient of the model and to predict new values of the response.

The practical behaviour of the method is also studied on simulated data sets. A real dataset illustration is performed in the environmental context of air quality.

Introduction

Literature on functional data is really wide, as attested by the numerous books on this subject these last years. The estimation and forecasting theories of linear processes in function spaces are developed in [START_REF] Bosq | Linear Processes in Function Spaces: Theory and Applications[END_REF]. A comprehensive introduction to functional data analysis can be found in [START_REF] Ramsay | Functional Data Analysis[END_REF]. In the focus of [START_REF] Ferraty | Nonparametric functional data analysis: Theory and practice[END_REF] are nonparametric approaches. Computational issues are explained in [START_REF] Ramsay | Functional Data Analysis with R and MATLAB[END_REF].

Nonparametric statistical methods for functional regression analysis, specifically the methods based on a Gaussian process prior in a functional space are discussed in [START_REF] Shi | Gaussian Process Regression Analysis for Functional Data[END_REF]. In [START_REF] Horváth | Inference for Functional Data with Applications[END_REF] inferential procedures based on functional principal components are considered. [START_REF] Zhang | Analysis of Variance for Functional Data[END_REF] mainly focuses on hypothesis testing problems about functional data. Among this, the functional linear model has received a special attention (see [START_REF] Ramsay | Some tools for functional data analysis[END_REF][START_REF] Cardot | Functional linear model[END_REF][START_REF] Cardot | Spline estimators for the functional linear model[END_REF][START_REF] Cai | Prediction in functional linear regression[END_REF][START_REF] Hall | Methodology and Convergence Rates for Functional Linear Regression[END_REF][START_REF] Crambes | Smoothing splines estimators for functional linear regression[END_REF][START_REF] Cardot | Thresholding projection estimators in functional linear models[END_REF][START_REF] Yuan | A reproducing kernel Hilbert space approach to functional linear regression[END_REF] for main references).

In this paper, we are interested in the functional linear model

Y = θ, X + ε, ( 1
)
where θ is the unknown function of the model, Y is a real variable of interest, ε is a centered real random variable representing the error of the model, with

finite variance E(ε 2 ) = σ 2 ε , and X is a functional covariate belonging to some functional space H endowed with an inner product ., . and its associated norm . . Usually, H is the space L 2 ([a, b]) of square integrable functions defined on some real compact [a, b] and the corresponding inner product is defined by f, g = b a f (t)g(t) dt for functions f, g ∈ L 2 ([a, b]). Without loss of generality, we consider our work on [0, 1]. Moreover, we assume that X and ε are independent.

All the previously cited works are devoted to analyse complete data, however, this is not the case in many interesting applications including for example survival data analysis. For this reason, we focus in this work on the problem of missing data (see [START_REF] Little | Statistical analysis with missing data[END_REF][START_REF] Graham | Missing data analysis and design[END_REF] for a wide introduction in the multivariate framework). This subject has been widely studied, in particular the way to impute missing data and the accuracy of this imputation according to the types of missing data: Missing Completely At Random (MCAR), Missing At Random (MAR) and Missing Not At Random (MNAR). Even if this problematic has received a lot of attention in a multivariate framework, it is not the case for the functional data framework. Our objective is to study the problem of combining regression imputation, missing data mechanisms and functional data analysis.

As far as we know, few results are available for the moment. In MAR setting, [START_REF] He | A functional multiple imputation approach to incomplete longitudinal data[END_REF] have explored this area by developing a functional multiple imputation approach modeling missing longitudinal response under a functional mixed effects model. They developed a Gibbs sampling algorithm to draw model parameters and imputations for missing values. Besides, [START_REF] Ferraty | Mean estimation with data missing at random for functional covariables[END_REF] have considered two kinds of mean estimates of a scalar outcome, based on a sample in which an explanatory variable is observed for every subject while responses are missing (which is the closest to our context). A weak convergence result was proved. In MCAR setting, [START_REF] Preda | The NIPALS Algorithm for Functional Data[END_REF] have adapted a methodology based on the NIPALS (Nonlinear Iterative Partial Least Squares) algorithm, which provides an imputation method for missing data, which have affected the functional covariates. In MNAR setting, [START_REF] Bugni | Specification test for missing functional data[END_REF] adapts a specification test for functional data with the presence of missing observations. His method is able to extract the information available in the observed portion of the data while being agnostic about the nature of the missing observations. In MAR and MCAR setting, [START_REF] Chiou | A functional data approach to missing value imputation and outlier detection for traffic flow data[END_REF] have recently proposed a nonparametric approach to missing value imputation and outlier detection for functional data. To our knowledge, there is no existing theoretical result in the case of functional linear model under missing assumption operating on the response variable, this problem only being until now the subject of studies in the multivariate framework (see for instance [START_REF] Manski | Identification problems in the social sciences[END_REF], [START_REF] Manski | Partial identification of probability distributions[END_REF]).

We carefully distinguish the missing data problem from a simple prediction problem. Indeed, the missing data mechanism involves a random variable (which indicates whether the response is missing or not) which plays a central role when obtaining our asymptotic results. This random variable and the variable X are dependent in the MAR case. This is also highlighted in [START_REF] Ferraty | Mean estimation with data missing at random for functional covariables[END_REF]. In this paper, we first propose an imputation method, based on the completely observed cases, to replace missing values in the response of the functional linear model. We get mean square error rates for these imputed values. Secondly, once the database is completed, we are able to estimate the unknown function θ of the model with the whole sample. This estimator can then be used for predicting other values of the response on a test set.

Combining missing data and functional variables offers a very large field of applications. Among all possible applications, environment is a core issue interesting many people for the future of our planet, in particular in the study of pollution indexes. The dataset we study here deals with temperature curves in some French cities to predict a specific pollution atmospheric index. The atmospheric index is missing in some cities in the northwest of France, for which the corresponding temperature curves (the explanatory variable) are mild, and leads to consider MAR data. The main objective is to get a map of the atmospheric index on the whole French territory.

The rest of the paper is organized as follows. Section 2 introduces the problem of functional linear model under missing assumption operating on the response variable and formulates our main results of the imputation method and of the mean square error for prediction of a new observation using the complete dataset. A simulation study is performed in Section 3. An environmental data illustration is presented in Section 4. Some preliminary lemmas, which are used in the proofs of the main results, are collected in Section 5.

Imputation of a missing value of the response

Functional principal components regression

Let us consider a sample (X i , Y i ) i=1,...,n independent and identically distributed with the same distribution as (X, Y ). An estimation of θ based on principal components analysis of the curves X 1 , . . . , X n has been studied in many papers, see for instance [START_REF] Cardot | Functional linear model[END_REF]. We recall below the construction of this estimator. Considering the covariance operator of X defined under the condition E X 2 < +∞ (which is supposed to be satisfied in the following) by

Γu = E X, u X ,
for all u ∈ H and its empirical version

Γ n u = 1 n n i=1 X i , u X i ,
we call (λ j ) j≥1 resp. λ j j≥1 the sequence of eigenvalues of Γ resp. Γ n and (v j ) j≥1 resp. ( v j ) j≥1 the sequence of eigenfunctions of Γ resp. Γ n .

The identifiability of model ( 1) is ensured as long as we assume that λ 1 > λ 2 > . . . > 0 (see [START_REF] Cardot | Functional linear model[END_REF]). Moreover, assuming that λ 1 > . . . > λ kn > 0 for some integer k n depending on n, the estimator of θ is defined by

θ = 1 n n i=1 kn j=1 X i , v j Y i λ j v j . (2)
A consistency result of this estimator is given in [START_REF] Cardot | Functional linear model[END_REF], while more recent results can be found in [START_REF] Cai | Prediction in functional linear regression[END_REF][START_REF] Hall | Methodology and Convergence Rates for Functional Linear Regression[END_REF]. In particular, [START_REF] Cardot | Functional linear model[END_REF] give technical conditions on the decreasing rate to zero of the eigenvalues λ j 's in order to ensure the consistency of the estimator.

Operatorial point of view

We notice in this subsection that the model (1) can be seen from an operatorial point of view. Indeed, we can write the model

Y = ΘX + ε, (3)
where Θ : H -→ R is a linear continuous operator defined by Θu = θ, u for any function u ∈ H. Let us consider ∆ n the cross covariance operator defined by

∆ n u = 1 n n i=1 X i , u Y i , for all u ∈ H.
Then, it is easily seen that an estimator Θ of Θ, satisfying Θ = θ, . , is given by

Θ = θ, . = Π kn ∆ n Π kn Γ n -1 , (4)
where Π kn is the projection operator onto the subspace Span( v 1 , . . . , v kn ).

Imputation principle 110

Now, we present the context of missing data. There can be many reasons for which missing data can appear: breakdown in a measurement process, a person who is not willing to answer to some question of a questionnaire, . . . We consider that some of the observations Y 1 , . . . , Y n are not available. We define the real variable δ and we consider the sample (δ i ) i=1,...,n such that δ i = 1 if the value Y i is available and δ i = 0 if the value Y i is missing, for all i = 1, . . . , n.

The data we observe are

{(Y i , δ i , X i )} n i=1 .
We consider that the missing values are MAR. The MAR assumption implies that δ and Y are conditionally independent given X. That is,

P (δ = 1 | X, Y ) = P (δ = 1 | X) . (5)
Note that the MAR assumption is much weaker than MCAR (for which

P (δ = 1 | X, Y ) = P (δ = 1
)), as it allows the missing data to possibly depend on the observed data and may be reasonable for many practical problems. As a consequence of this MAR assumption, the variable δ (the fact that an observation is missing) is independent of the error of the model , conditionally on X.

In the following, the number of missing values in the sample is denoted

m n = n i=1 1 1 {δi=0} . (6)
Then, to impute a missing value, say Y (where is a given integer between 1 and n), a simple way is to consider complete case analysis (see for instance [START_REF] Little | Statistical analysis with missing data[END_REF][START_REF] Cheng | Nonparametric Estimation of Mean Functionals with Data Missing at Random[END_REF][START_REF] Wang | Semiparametric Regression Analysis with Missing Response at Random[END_REF][START_REF] Mojirsheibani | Nonparametric curve estimation with missing data: A general empirical process approach[END_REF][START_REF] Van Buuren | Flexible Imputation of Missing Data[END_REF]). This regression imputation method uses the pairs of observed data to define the estimator of the model coefficient. More precisely, we define

Y ,imp = 1 n -m n n i=1 i = kn j=1 X i , v j X , v j δ i Y i λ j . (7)
From the operatorial point of view, the imputation of the missing value Y comes back to

Y ,imp = Π kn,obs ∆ n,obs Π kn,obs Γ n,obs -1 X , (8)
where

Γ n,obs = 1 n-mn n i=1 X i , . δ i X i , ∆ n,obs = 1 n-mn n i=1 X i , . δ i Y i and Π kn,obs
is the projection operator onto the subspace span( v 1,obs , . . . , v kn,obs ) where v 1,obs , . . . , v kn,obs are the k n first eigenfunctions of the covariance operator Γ n,obs . Now we give our main results. We consider the following assumptions.

(A.1) We assume that there exists a convex function λ such that λ(j) = λ j for all j ≥ 1 that continuously interpolates the λ j 's between j and j + 1.

(A.2) There exists a positive constant C such that

E X 4 ≤ C.
Our assumptions are quite classic in this context. Assumption (A.1) is similar to an assumption from [START_REF] Crambes | Asymptotics of prediction in functional linear regression with functional outputs[END_REF]. It is a mild condition that allows a large class of decreasing rate of eigenvalues for the covariance operator Γ, for example polynomial decay or exponential decay (see example 1 below, in page 7, for more details). Assumption (A.2) holds for many processes X (Gaussian processes, bounded processes) and can also be found for example in [START_REF] Cardot | Functional linear model[END_REF]. Then, we give our main results.

Remark 1. Notice that the assumptions (A.

E Y ,imp -θ, X 2 = +∞ j=kn+1 ΘΓ 1/2 v j 2 + σ 2 ε k n n -m n + o k n n -m n .
Moreover, for the aggregate mean square error of all the imputed values, we have

n =1 (1-δ )E Y ,imp -θ, X 2 = m n +∞ j=kn+1 ΘΓ 1/2 v j 2 + σ 2 ε k n m n n -m n + o k n m n n -m n .
In order to precise the convergence rate of the imputed value Y ,imp to the real one θ, X , we need an additional notation. For a function ϕ : R + -→ R + and a positive real number L, we define

C(ϕ, L) = T : H -→ R / ∀j ≥ 1, T v j ≤ L ϕ(j) .
Note that simple cases satisfy the fact that ΘΓ 1/2 belongs to C(ϕ, L). For example, consider the operator Θ expressed in the eigenfunctions basis (v j ) j≥1

such that Θu = +∞ j=1 θ j v j , u for any u ∈ H, with θ j going to zero as j goes to infinity. Hence there exists a bound L such that θ j ≤ L for any j ≥ 1 and

ΘΓ 1/2 v j = θ j λ j ≤ L λ j .
Remark 2. We introduce two notations to compare the magnitudes of two functions ũ(x) and ṽ(x) as the argument x tends to a limit ˜ (not necessarily finite).

The notation ũ(x)

∼ x→ ˜ ṽ(x), stands for lim x→ ˜ ũ(x) ṽ(x) = 1,
and the notation ũ(x)

x→ ˜ ṽ(x) denotes that |ũ(x)/ṽ(x)| remains bounded as

x → ˜ . Theorem 2.2. Let L = ΘΓ 1/2 ∞ and ϕ the function defined by ϕ(j) = (ΘΓ 1/2 vj ) 2 L 2
for all j ≥ 1 that continuously interpolates the ϕ(j)'s between j and j + 1. Under assumptions (A.1)-(A.2), the operator ΘΓ 1/2 belongs to C(ϕ, L) and

E (Y ,imp -θ, X ) 2 ∼ n→+∞ 2σ 2 ε k n n -m n ,
where k n is the solution of the equation in

x +∞ x ϕ(t) dt = σ 2 ε L 2 (n -m n ) x. (9)
Again, for the aggregate mean square error of all the imputed values, we have

n =1 (1 -δ )E (Y ,imp -θ, X ) 2 ∼ n→+∞ 2σ 2 ε k n m n n -m n .
Remark 3. Notice that the equation ( 9) has a unique solution (the left and right hand sides are descreasing and increasing in x, respectively). The practical resolution of equation ( 9) to get k n seems quite complicated due to the computation of L. In order to solve this problem, we will use other ways to select the optimal number of principal components (see Section 3 below).

The last result giving the convergence rate of the imputed value Y ,imp is similar to the convergence rate obtained in [START_REF] Crambes | Asymptotics of prediction in functional linear regression with functional outputs[END_REF] (who considered the case of a completely observed functional response). The rate is simply affected by the number m n of missing values. We precise the resulting rate of convergence in the following examples.

Example 1. We consider two different functions ϕ such that ϕ pol (j

) = C α j -(2+α)
and ϕ exp (j) = D α exp(-αj) where C α and D α are positive constants and α > 0.

Then the solution of equation ( 9) is

           k n,pol ∼ n→+∞ CαL 2 (1+α)σ 2 ε 1/(2+α) n 1/(2+α) , if ϕ = ϕ pol , k n,exp n→+∞ log n α , if ϕ = ϕ exp .
For ϕ = ϕ pol , the result of Theorem 2.2 becomes

E (Y ,imp -θ, X ) 2 ∼ n→+∞ 2 σ 2 ε (1+α)/(2+α) C α L 2 1 + α 1/(2+α) n 1/(2+α) n -m n ,
for a single imputation and

n =1 (1-δ )E (Y ,imp -θ, X ) 2 ∼ n→+∞ 2 σ 2 ε (1+α)/(2+α) C α L 2 1 + α 1/(2+α) n 1/(2+α) m n n -m n ,
for the aggregate error of all the imputed values.

For ϕ = ϕ exp , the result of Theorem 2.2 becomes

E (Y ,imp -θ, X ) 2 n→+∞ 2σ 2 ε log n α(n -m n ) ,
for a single imputation and

n =1 (1 -δ )E (Y ,imp -θ, X ) 2 n→+∞ 2σ 2 ε m n log n α(n -m n ) ,
for the aggregate error of all the imputed values.

Example 2. To precise in more specific cases our convergence rates, we consider three different levels of missing data: (i) when the number of missing data m n is negligeable compared to the sample size, that is m n = a n n with a n going to zero as n goes to infinity, (ii) when the number of missing values is proportional to the sample size, that is m n = ρn with 0 < ρ < 1, and (iii) when the number of observed values is negligeable compared to the sample size, that is

u n := n-m n = o(n).
We can sum up all the rates of convergence for the single imputation mean square error (Table 1) and for the aggregate mean square error (Table 2).

We can see that missing data do not affect the convergence rate for a single imputed value when there are not too many missing values (m n = o(n) or m n = ρn). The rate 1/n (1+α)/(2+α) matches the usual optimal rates in this context.

The rate log n/αn is not exact but obviously sharp since parametric up to a logarithm. It is no more the case when the number of missing values is high (m n ∼ n), the convergence rate is affected. For the aggregate error of several Single imputation mean square error convergence rates, where Kα := .

2 σ 2 ε (1+α)/(2+α) CαL 2 1+α 1/(2+α) . ϕ = ϕ pol ϕ = ϕ exp m n := a n n = o(n) ∼ n→+∞ K α n -(1+α)/(2+α) n→+∞ 2σ 2 ε log n αn m n = ρn ∼ n→+∞ K α (1 -ρ) 1/(2+α) n -(1+α)/(2+α) n→+∞ 2σ 2 ε log n α(1-ρ)n u n := n -m n = o(n) ∼ n→+∞ K α u -(1+α)/(2+α) n n→+∞ 2σ 2 ε log n αun
ϕ = ϕ pol ϕ = ϕ exp m n := a n n = o(n) ∼ n→+∞ K α a n n 1/(2+α) n→+∞ 2σ 2 ε an log n α m n = ρn ∼ n→+∞ K α ρ(1 -ρ) 1/(2+α) n 1/(2+α) n→+∞ 2σ 2 ε ρ log n α(1-ρ) u n := n -m n = o(n) ∼ n→+∞ K α nu -(1+α)/(2+α) n n→+∞ 2σ 2 ε n log n αun
imputed values, when there are not too many missing values (m n = o(n)), the number of missing values plays a crucial role, since the convergence depends on the fact that a n n 1/(2+α) or a n log n go to zero as n goes to infinity. In other cases (m n = ρn or m n ∼ n), missing data affect the convergence of the aggregate error term for several imputed values, since it cannot converge to zero.

Estimation of θ and prediction of future values

Once the database being reconstructed, we can use the full database to estimate the functional coefficient θ of the model (directly inspired from ( 2)) (see also [START_REF] Chu | Nonparametric regression estimation with missing data[END_REF]), namely

θ = 1 n n i=1 kn j=1 X i , v j Y i λ j v j , (10)
where

Y i = Y i δ i + Y i,imp (1 -δ i ) for all i = 1, . . . ,
Y new = X new , θ = 1 n n i=1 kn j=1 X i , v j X new , v j Y i λ j . ( 11
)
We give below a result allowing to control the mean square prediction error of Y new .

Theorem 2.3. Under the assumptions of Theorem 2.1, if we additionally assume that m n = o(n) (that is m n = a n n with a n going to zero as n goes to infinity) and

a 2 n n = o(1), then E Y new -θ, X new 2 = +∞ j=kn+1 ΘΓ 1/2 v j 2 + O k n n .
Remark 4. This result shows that, under the condition that there are not too many missing values, the convergence rate of the mean square error prediction of a new value of the covariate remains the same compared to the non missing values case.

Simulations

To observe the behavior of our estimator in practice, this section considers a simulation study.

Models

Two models are considered:

M odel 1 : Y = 1 0 sin(4πt)X t dt + , (12)
M odel 2 : Y = 1 0 log(15t 2 + 10) + cos(4πt) X t dt + , (13)
where the error is a Gaussian noise : ∼ N (0, σ ) and

• in equation ( 12), X := {X t } t∈[0,1] is the standard Brownian motion.

• In equation ( 13), X := {X t } t∈[0,1] is a Gaussian process where the covariance function is defined as

cov(X t , X t ) = exp(- |t -t | 2 0.2 ).
The simulation aims at considering processes X with different regularities (the standard Brownian motion being the case of the less smooth) in order to see if it has an impact on the results.

All the procedures described below were implemented by using the R software:

the trajectories of X i , 1 ≤ i ≤ n,

Criteria

The criteria we used are the following. Criteria 1 and 2 are related to the imputation step with the training samples, criteria 3 and 4 are related to the prediction step with the test samples, and criteria 5 is related to the estimation step with the reconstructed database.

• Criterion 1: the mean square errors (M SE) averaged over S samples

M SE = 1 S S j=1 M SE(j),
where

M SE(j) = 1 mn n =1 Y j ,imp -θ, X j 2 (1 -δ )
is the mean square error computed on the j th simulated sample, j ∈ {1, . . . , S}.

• Criterion 2: the ratio respect to truth between the mean square prediction error and the mean square prediction error when the true mean is known averaged over S samples

RT = 1 S S j=1 RT (j),
where

RT (j) = n =1 Y j ,imp -θ, X j 2 (1 -δ ) n =1 j 2 (1 -δ )
is the ratio between the mean square prediction error and the mean square prediction error when the true mean is known, computed on the j th simulated sample.

• Criterion 3: the mean square errors (M SE ) averaged over S samples

M SE = 1 S S j=1 M SE (j),
where M SE (j) = 1 n1 n+n1 =n+1 Y j -θ, X j 2 is the mean square error computed on the j th simulated sample, j ∈ {1, . . . , S}.

• Criterion 4: the ratio respect to truth between the mean square prediction error and the mean square prediction error when the true mean is known averaged over S samples

RT = 1 S S j=1 RT (j),
where

RT (j) = n+n1 =n+1 Y j -θ, X j 2 n+n1 =n+1 j 2
is the ratio between the mean square prediction error and the mean square prediction error when the true mean is known, computed on the j th simulated sample.

• Criterion 5: the mean square errors (M SE ) averaged over S samples

M SE = 1 S S j=1 M SE (j),
where

M SE (j) = θ j -θ 2
is the square error of estimation computed on the j th simulated sample. The M SE criterion is decomposed into variance and square bias in our results.

Notice that all the criteria tend to zero when the sample size tends to infinity.

RT and RT are rescaled versions of M SE and M SE if we substitute the denominator by its limit (specifically, M SE(j) = RT (j)σ 2 ).

Methodology

We use a smoothed version of the estimator (2) based on the Smooth Principal Components Regression (SPCR) [START_REF] Cardot | Spline estimators for the functional linear model[END_REF]. We use a regression spline basis with parameters: the number κ of knots of the spline functions, the degree q of spline functions and the number m of derivatives. Let us remark that, with appropriate conditions, all the theoretical results obtained in section 2 will also apply to the SPCR estimation. For example, we assume that the estimator θ has r derivatives for some integer r and θ (r ) satisfies, for some ν ∈]0, 1]

θ (r ) (t 1 ) -θ (r ) (t 2 ) ≤ C |t 1 -t 2 | ν , for all t 1 , t 2 ∈ [0, 1].
If we denote r = r + ν and if we assume that the degree q of the splines is such that q ≥ r, then sup

t∈[0,1] θ(t) -S κ,q ( θ)(t) = O κ -r ,
where S κ,q ( θ) is the spline approximation of θ (see [START_REF] De Boor | A practical guide to splines[END_REF]). In other words, any of the convergence results obtained in Section 2 can be transposed to the smoothed version of the estimators.

Here, we have fixed the number of knots to be 20, the degree has been chosen to be 3 and the number of derivatives was fixed to the moderate value of 2. The choice of these parameters is not the most important in our study, especially in comparison with the choice of the number of principal components.

In this subsection, we show firstly how to determine the number of missing data. Secondly, we present a procedure to choose the optimal tuning parameter (the best dimension k * n of the projection space for the SPCR).

Missing data simulation scenario

To determine the number of missing data in our simulations, we have adopted the following scenario. In the MAR case, we simulate δ according to the logistic functional regression. The variable δ follows the Bernoulli law with parameter p(X) such that log p(X)

1 -p(X) = α 0 , X + ct,
where α 0 (t) = sin(2πt) for all t ∈ [0, 1] and ct is a constant allowing to take different levels of missing data. We take ct = 2 for around 12.5% of missing data, ct = 1 for around 27.4% of missing data and ct = 0.2 for around 44.9% of missing data. Notice that, in the MCAR case, we simulate δ with the Bernoulli law with parameter p(X) := p = 0.9 (10% of missing data), p(X) := p = 0.75 (25% of missing data) or p(X) := p = 0.6 (40% of missing data).

Criteria for optimal parameter selection

We focus on the procedure allowing to select the optimal tuning parameter.

We consider a Generalized Cross Validation (GCV) criterion versus a Cross Validation (CV) criterion and K-fold Cross Validation (K-fold CV) criterion and we select the optimal tuning parameter k * n by minimizing these criteria.

The GCV procedure is known to be computationally fast. The CV, K-fold CV and GCV criteria are respectively given as follows for imputation

CV(k n ) = 1 n -m n n i=1 ( Ŷ [-i] i -θ, X i ) 2 δ i , K-fold CV(k n ) = 1 K K k=1 |B k | -1 i∈B k ( Ŷ [-B k ] i -θ, X i ) 2 δ i , GCV(k n ) = (n -m n ) n i=1 ( Ŷi -θ, X i ) 2 δ i ((n -m n ) -k n ) 2 .
The analogous criteria are given as follows for prediction

CV(k n ) = 1 n n i=1 ( Ŷ * [-i] i -θ, X i ) 2 , K-fold CV(k n ) = 1 K K k=1 |B k | -1 i∈B k ( Ŷ * [-B k ] i -θ, X i ) 2 , GCV(k n ) = n n i=1 ( Ŷ * i -θ, X i ) 2 (n -k n ) 2 , where Ŷ [-i] i and Ŷ [-B k ]
i respectively mean that the value of Y i is predicted using the whole sample except the i th observation or except the set of observations indexed in B k . In the same way Ŷ * [-i] i and Ŷ * [-B k ] i respectively mean that the value of Y i is predicted using the whole sample except the i th observation or except the set of observations indexed in B k . The data set is randomly partitioned into K equally sized (as equal as possible) subsets ∪ K k=1 B k such that B j ∩ B k = ∅ (j = k). In practice, often K = 5 or K = 10 are used. In our case, the K-fold CV splits are chosen in a special deterministic way. For imputation, we consider

K-fold CV(k n ) = 1 K K k=1 ((n -m n )/K) -1 nk/K i=(n(k-1))/K +1 ( Ŷ [-k] i -θ, X i ) 2 δ i .
The analogous criterion is given as follows for prediction

K-fold CV(k n ) = 1 K K k=1 (n/K) -1 nk/K i=(n(k-1))/K +1 ( Ŷ * [-k] i -θ, X i ) 2 .
In order to illustrate the advantage of the GCV criterion, we compared the computational times to obtain the tuning parameter with the three criteria on A.11 in the appendix. The GCV criterion shows a clear advantage with regard to computational time compared with the CV and K-fold criteria. In addition, we see that the three criteria behave in the same way and select the same optimal projection dimension (see Fig. 1 and2) for both models (under n = 1000 and p = 100). Notice that the GCV criterion (faster to compute) has been used in different simulations. We show on Fig. 3 dimension (k n = 4, 6, 8, 12, 16) and (k n = 2, 3, 5, 7, 8), respectively, by using the GCV criterion (used for its computational efficiency). We have chosen a percentage of missing values equal to 45.8518% for model 1 and equal to 46.8888% for model 2 (we obtain this rate with ct = 1 for both models).

Analysis of results

In this subsection, we analyse the results of the criteria presented in the previous subsection. Both MAR and MCAR context were considered. We only show the results for MAR and the results for MCAR are available on demand.

The different results given in Appendix A. Tables A.5, A.6 give the mean and 330 standard deviation errors for the imputed values on training samples for both models. Tables A.7, A.8 give the mean and standard deviation errors for the predicted values on test samples for both models. Tables A.9, A. mean and standard deviation errors for the estimation of θ using the fullfilled database with imputed values for both models. We can see that the errors increase when the rate of missing data increases. Similarly, the errors decrease as the size of the sample increases. When we compare the case of MAR and MCAR, we see that the error in case of MAR is slightly higher that in the MCAR case. Moreover, we can see that the regularity of the process X does not have a crucial impact on the results at least on these simulated examples.

All the results in these simulations are in accordance with what we can expect and confirm the theoretical results obtained in the previous section.

Illustration

In order to illustrate the contribution of our approach in functional prediction setting when the covariates are functions and some observations of the 20 real response are missing, we present in this section an environmental dataset application.

We start by describing the dataset. The functional covariate X is a daily temperature curve in some cities in France (from May 7, 2015 at 4 pm up to May 8, 2015 at 3 pm) obtained from www.meteociel.fr. This daily continu-350 ous curve is observed at some discretization points (here, at 24 discretization points, every hour). The graphical display of this daily temperature curves can be observed in Fig 5 . The response variable Y is an atmospheric index of air quality called ATMO (for a detailed description of this atmospheric index, see www.atmo-france.org). Its values range from 1 (very good quality of air) to 10 (very bad quality of air). Though these values are discrete, we will consider that Y is a continuous approximation. We obtained the values of the atmospheric index on May 8, 2015, for these same cities, from www2.prevair.org. Furthermore, we added some cities for which the temperature curve is available but the atmospheric index is missing. Notice that the response is missing for mild temperature curves cities: the fact that the value of the response variable Y is missing for these cities depends on the temperature curve X, and thus we consider the MAR case. We also refer the reader to the paper [START_REF] Junninen | Methods for imputation of missing values in air quality datasets[END_REF] for more discussions about missing data mechanism when dealing with air quality data.

In particular, this paper highlights the fact that air quality missing data can be considered as MAR. Fig 6 illustrates the selected cities in our study, the blue cities are given when the response variable Y is missing and the red cities are given when the response variable Y is observed. It is of primary importance to get a map of the atmospheric index on the whole French territory, and thus to impute missing data.

We have built a sample of 78 pairs {(Y i , X i )} 78 i=1 , where we have 8 missing values of the variable Y (the Y i s, i = 71, . . . , 78, are missing). Our goal is to impute these missing values {Y i } 78 i=71 .

We have fixed the number of knots to be 20, the degree of splines has been chosen equal to 3 and the number of derivatives was fixed to the moderate value of 2. Then, we use the GCV criterion to find the best parameter of projection dimension k n trying growing sequences: k n = 2, 3, . . . , 21, 22. In order to see the impact of missing data on this dataset, we have randomly drawn 700 tests samples in the initial sample and computed prediction errors on these tests samples, using the remaining of the sample as training sample. Results are given in Table 3. Here again, the more we have missing data in the training set, the more the prediction error on the test sample is. Now, we come back to the initial goal, imputing the missing data. The 4 gives the imputed values of the missing data. We see imputed values mainly around 4, which is a moderate value for the atmospheric index corresponding to a good quality of air. It is in accordance with the fact that these cities have moderate temperature curves. We can mention two particular 390 cases. The highest imputed value (4.161) corresponds to the city of Angers, and in parallel, we can see that the temperature curve of this city becomes high at the end of May 8. On the contrary, the lowest imputed value (3.491) corresponds to the city of Quimper, and the temperature curve of this city presents few variations along the 24 hours. 395

Proof of the results

Proof of Theorem 2.1

We begin with the following decomposition

∆ n,obs = 1 n -m n n i=1 X i , . δ i ΘX i + 1 n -m n n i=1 X i , . δ i ε i = Θ Γ n,obs + U n,obs , with U n,obs = 1 n-mn n i=1 X i , . δ i ε i .
Then, ε being independent from X and δ (MAR assumption), we deduce

E (Y ,imp -θ, X ) 2 = E Θ Π kn,obs X -ΘX 2 +E 1 n -m n n i=1 X i , Π kn,obs Γ n,obs -1 X δ i ε i 2 ≤ 2E Θ Π kn,obs X -ΘΠ kn,obs X 2 +2E (ΘΠ kn,obs X -ΘX ) 2 +E 1 n -m n n i=1 X i , Π kn,obs Γ n,obs -1 X δ i ε i 2 ,
where Π kn,obs is the projection onto the subspace span(v 1,obs , . . . , v kn,obs ) where v 1,obs , . . . , v kn,obs are the k n first eigenfunctions of the covariance operator Γ n,obs .

For a single imputation, the end of the proof of Theorem 2.1 is based on the following lemmas. For the aggregate error term of m n imputed values, it is just a sum of m n terms that behave like the term for single imputation.

Lemma 5.1. We have

E Θ Π kn,obs X -ΘΠ kn,obs X 2 = o λ kn k 2 n n -m n + k n n -m n .
Lemma 5.2. We have

E 1 n -m n n i=1 X i , Π kn,obs Γ n,obs -1 X δ i ε i 2 = σ 2 ε k n n -m n + o k n n -m n .
Lemma 5.3. We have

E (ΘΠ kn,obs X -ΘX ) 2 = +∞ j=kn+1 ΘΓ 1/2 v j 2 .

Proof of Lemma 5.1

Writing X in the basis (v j ) j≥1 , we obtain

E Θ Π kn,obs X -ΘΠ kn,obs X 2 = +∞ j=1 +∞ j =1
E X , v j X , v j Θ Π kn,obs -Π kn,obs v j Θ Π kn,obs -Π kn,obs v j .

Noticing that the variable X corresponds to the missing data Y hence independent of Π kn,obs , we get

E Θ Π kn,obs X -ΘΠ kn,obs X 2 = +∞ j=1 +∞ j =1
Γv j , v j E Θ Π kn,obs -Π kn,obs v j Θ Π kn,obs -Π kn,obs v j With these notations, denoting by ι the complex number such that ι 2 = -1, the difference between the projection operators Π kn,obs and Π kn,obs can be written Π kn,obs -Π kn,obs = 1 2πι

kn m=1 Bm Λ(z) Γ -Γ n,obs Λ(z)dz, where Λ(z) = (zI -Γ) -1 . Noticing that Λ(z)v j = 1 z-λj v j , we deduce Θ Π kn,obs -Π kn,obs v j = 1 2πι kn m=1 Θ Bm Λ(z) Γ -Γ n,obs dz z -λ j = 1 2πι kn m=1 Θ Bm +∞ j =1 Γ -Γ n,obs v j , v j v j (z -λ j )(z -λ j ) dz.
Still using the results from [START_REF] Crambes | Asymptotics of prediction in functional linear regression with functional outputs[END_REF], we have

kn m=1 Bm dz (z -λ j )(z -λ j ) =                0, if j, j > k n , 0, if j, j ≤ k n , (λ j -λ j ) -1 , if j ≤ k n < j , (λ j -λ j ) -1 , if j ≤ k n < j.
hence we deduce 420

E Θ Π kn,obs X -ΘΠ kn,obs X 2 = E    1 4π 2 kn j=1 λ j   +∞ j =kn+1 Γ -Γ n,obs v j , v j λ j -λ j Θv j   2    +E    1 4π 2 +∞ j=kn+1 λ j   kn j =1 Γ -Γ n,obs v j , v j λ j -λ j Θv j   2    .
In the following, C corresponds to a generic constant. We denote E(A) and E(B) the above two terms. We start with the computation of E(A). Using the same technique as in [START_REF] Crambes | Asymptotics of prediction in functional linear regression with functional outputs[END_REF], we get the following bound

E Γ -Γ n,obs v j , v j Γ -Γ n,obs v j , v r ≤ C n -m n λ j λ j λ r ,
noticing that the n rate of convergence given in [START_REF] Crambes | Asymptotics of prediction in functional linear regression with functional outputs[END_REF] is here transformed into the n -m m rate because we use Γ n,obs with n -m m observed data. Hence we deduce

E   Γ -Γ n,obs v j , v j λ j -λ j Θv j   2 = +∞ j =kn+1 +∞ r=kn+1 E Γ -Γ n,obs v j , v j Γ -Γ n,obs v j , v r (λ j -λ j )(λ j -λ r ) Θv j Θv r ≤ Cλ j n -m n   +∞ j =kn+1 λ j λ j -λ j Θv j   2 .
Coming back to the computation of E(A), we can write (using Lemma 12 in [START_REF] Crambes | Asymptotics of prediction in functional linear regression with functional outputs[END_REF])

E(A) ≤ C n -m n kn j=1 λ 2 j λ kn+1 (λ j -λ kn+1 ) 2   +∞ j =kn+1 Θv j   2 ≤ Cλ kn+1 n -m n kn j=1 (k n + 1) 2 (k n + 1 -j) 2   +∞ j =kn+1 Θv j   2 ≤ Cλ kn+1 (k n + 1) 2 n -m n kn j=1 1 j 2   +∞ j =kn+1 Θv j   2 .
As θ ∈ L 2 ([0, 1]) (hence θ is integrable), we finally get

E(A) ≤ Cλ kn k 2 n n -m n a n ,
where (a n ) n≥1 is a sequence of real numbers going to zero as n goes to infinity.

430

We are now interested in the computation of E(B). Beginning in the same way as E(A) and still using Lemma 12 in [START_REF] Crambes | Asymptotics of prediction in functional linear regression with functional outputs[END_REF], we get

E(B) ≤ C n -m n +∞ j=kn+1 λ 2 j   kn j =1 λ j λ j -λ j Θv j   2 ≤ C n -m n +∞ j=kn+1 λ j   kn j =1 λ j λ j -λ j Θv j   2 ≤ C n -m n +∞ j=kn+1 λ j j j -k n 2   kn j =1 Θv j   2 .
Now, again with the integrability of θ and the fact that

+∞ j=kn+1 λ j j j -k n 2 ≤ Ck n b n ,
with (b n ) n≥1 going to zero as n goes to infinity (see [START_REF] Crambes | Asymptotics of prediction in functional linear regression with functional outputs[END_REF] Let us denote

T n = 1 n -m n n i=1 X i , Π kn,obs Γ n,obs -1 X δ i ε i .
We can write

T 2 n = 1 (n -m n ) 2 n i=1 X i , Π kn,obs Γ n,obs -1 X 2 δ 2 i ε 2 i + 1 (n -m n ) 2 n i=1 n i =1 i =i X i , Π kn,obs Γ n,obs -1 X X i , Π kn,obs Γ n,obs -1 X δ i δ i ε i ε i .
From the independence between ε and X and the MAR assumption, the expec-440 tation of the second term above is zero, hence

E T 2 n = 1 n -m n E X i , Π kn,obs Γ n,obs -1 X 2 δ 2 i ε 2 i = σ 2 ε n -m n E X i , Π kn,obs Γ n,obs -1 X 2 δ 2 i ,
the index i corresponding to an observed data in the sample (and consequently δ i = 1 for this observation). We finally get

E T 2 n = σ 2 ε n -m n E X i , Π kn,obs Γ n,obs -1 X 2 .
Following the same lines of the proof of Proposition 17 and Lemma 19 in [START_REF] Crambes | Asymptotics of prediction in functional linear regression with functional outputs[END_REF],

we obtain 445

E X i , Π kn,obs Γ n,obs -1 X 2 = k n + o (k n ) ,
which achieves the proof of the Lemma.

S 2 n = 1 n 2 n i=1 X i , Π kn Γ n -1 X new 2 (ε i ) 2 + 1 n 2 n i=1 n i =1 i =i X i , Π kn Γ n -1 X new X i , Π kn Γ n -1 X new ε i ε i .
We notice that, for i = i , we have

E (ε i ε i ) ≤ 4E(Y i,imp -Y i ) 2 ≤ 8 E(Y i,imp -θ, X i ) 2 + σ 2 ε .
This bound and the lines of the proof of Lemma 5.2 give

E 1 n n i=1 X i , Π kn Γ n -1 X new ε i 2 = O (n -m n )k n n 2 + m 2 n k 2 n n 2 . (16)
Now, combining relations ( 14), ( 15) and ( 16) and the fact that m n = o(n) and

m 2 n k n = O(n) (due to a 2 n n = o( 1
)), we get the desired result.

  in the two models are discretized in p = 100 equidistant points, values of Y are computed using integration by rectangular interpolation, the variability of noise is such that σ = τ * Var 1 0 θ(t)X(t)dt ≈ 0.2. Note that some Monte Carlo experiments are achieved to determine the values of τ : τ ≈ 21.726 for the model 1 (low level of noise) and τ ≈ 0.048 for the model 2 (high level of noise), the sample sizes are respectively n = 100, 300 and 1200 for the training sets (X 1 , Y 1 ), . . . , (X n , Y n ) and n 1 = 50, 150 and 600 for the test sets (X n+1 , Y n+1 ), . . . , (X n+n1 , Y n+n1 ).

a growing sequence of dimension k n = 2 ,

 2 . . . , 22. The characteristics of the computer used to perform these computations were McBook pro: Processor 2.66 GHz intel core 2 Duo, Memory 4 Gb 1067 MHz DDR3. The computational times are displayed in Table

Figure 1 :

 1 Figure 1: GCV, CV and K-fold criteria for different values of dimension kn in model 1 : best dimension k * n = 8 and M SE' (×10 4 ) = 1.6640 (in GCV criterion case), best dimension k * n = 6 and M SE (×10 4 ) = 2.3081 (in 5-fold CV criterion case), best dimension k * n = 8 and M SE (×10 4 ) = 1.9584 (in 10-fold CV criterion case), best dimension k * n = 8 and M SE (×10 4 ) = 1.6598 (in CV criterion case), for the model 1 .

  and Fig. 4 different estimates of the slope function of the M odel 1 and M odel 2 (under n = 1000 and p = 100) with different values of 0.000 0.005 0.010 0.015 0.020 0.025 0

Figure 2 :

 2 Figure 2: GCV, CV and K-fold criteria for different values of dimension kn in model 2 : best dimension k * n = 5 and M SE (×10 4 ) = 3.7589 (in GCV criterion case), best dimension k * n = 5 and M SE (×10 4 ) = 4.2132 (in 5-fold CV criterion case), best dimension k * n = 5 and M SE (×10 4 ) = 3.9758 (in 10-fold CV criterion case), best dimension k * n = 5 and M SE (×10 4 ) = 3.7270 (in CV criterion case), for the model 2 .

  325

Figure 3 :

 3 Figure 3: Plots of the true slope function (solid black) and estimates with different values of dimension kn in model 1 . The plots of estimates slope function with best dimension k * n = 8 (solid red), with dimension kn = 4 (dotted), with dimension kn = 6 (dashed), with dimension kn = 12 (dotdashed), with dimension kn = 16 (twodash).

Figure 4 :

 4 Figure 4: Plots of the true slope function (solid black) and estimates with different values of dimension kn in model 2 . The plots of estimates slope function with best dimension k * n = 5 (solid red), with dimension kn = 2 (dotted), with dimension kn = 3 (dashed), with dimension kn = 7 (dotdashed), with dimension kn = 8 (twodash).

Figure 5 :

 5 Figure 5: Plot of the 78 daily temperature curves (the blue curves are given when the response variable Y is missing).

Figure 6 :

 6 Figure 6: Map of France locating the selected cities of our study: the cities are red when the variable Y is observed and the cities are blue when the variable Y is missing.

λ 2 .

 2 j E Θ Π kn,obs -Π kn,obs v j Now, following the proof of Proposition 15 in [11], for any m ≥ 1 we denote B m the oriented circle of the complex plane with center λ m and radius ρ m /2 where ρ m = min (λ m -λ m+1 , λ m-1 -λ m ) for m ≥ 2 and ρ 1 = λ 2 -λ 1 . With the 415 convexity assumption (A.1), we actually have ρ m = λ m -λ m+1 for all m ≥ 1.

5 . 3 .

 53 p.19 in the proof of Proposition 15), we conclude 435 E(B) ≤ Ck n n -m n b n , and this achieves the proof of Lemma 5.1. Proof of Lemma 5.2

Table 1 :

 1

Table 2 :

 2 Aggregate imputation mean square error convergence rates, where Kα :=

	ε 2 σ 2	1+α (1+α)/(2+α) CαL 2	1/(2+α)

  [START_REF] Crambes | Smoothing splines estimators for functional linear regression[END_REF] give the

	4
		True slope function
	3	Estimed slope function with best dimension 8 Estimed slope function with dimension 4
		Estimed slope function with dimension 6
		Estimed slope function with dimension 12
		Estimed slope function with dimension 16
	2
	1
	0
	-1
	-2
	0	0.2	0.4	0.6	0.8	1

Table 3 :

 3 Real data set: prediction errors over 700 drawn samples.

		n = 78, 8 missing data, 70 observed data
	Test sets	n/4	n/3	n/2
	Rate of missing data (%)	13	15	20
	M SE × 10 2	24.5650	25.5172	29.7827
		(8.4750)	(8.1444)	(15.0889)
	minimum value of the GCV criterion is reached for k * n = 5 and M SE (×10 2 ) =
	20.791. Table

Table 4 :

 4 Imputed values of the missing response variable.

	Missing values of Y	Y 71	Y 72	Y 73	Y 74	Y 75	Y 76	Y 77	Y 78
	Imputed values	4.161 3.496 3.850 3.758 3.590 3.491 3.990 3.821

Table A .

 A 7: MAR (M odel 1 ): Predicted values mean errors and standard deviations for samples with different sizes discretized in p = 100 equidistant points based on 500 simulation replications.

			n + n 1 = 150
	Rate of missing data (%)
	Mean	12.520	27.420	44.882
	Median	13	27	45
	SD	3.307	4.515	5.038
	Criterion 3: [M SE × 10 3 ]	2.3383	2.7173	3.1939
		(1.4987)	(1.8390)	(2.0391)
	Criterion 4: [RT × 10 2 ]	5.9523	6.9769	8.2677
		(3.7338)	(4.9933)	(5.6516)
			n + n 1 = 450
	Rate of missing data (%)
	Mean	12.433	27.456	45.209
	Median	12.333	27.333	45.333
	SD	1.877	2.487	3.041
	Criterion 3: [M SE × 10 3 ]	0.8453	0.9984	1.3046
		(0.5530)	(0.6729)	(0.8897)
	Criterion 4: [RT × 10 2 ]	2.1534	2.5348	3.3255
		(1.3984)	(1.6629)	(2.2417)
			n + n 1 = 1800
	Rate of missing data (%)
	Mean	12.529	27.536	45.213
	Median	12.500	27.500	45.250
	SD	0.934	1.280	1.355
	Criterion 3: [M SE × 10 3 ]	0.2295	0.2746	0.3474
		(0.1282)	(0.1512)	(0.1982)
	Criterion 4: [RT × 10 2 ]	0.5756	0.6887	0.8699
		(0.3165)	(0.3753)	(0.4888)

  Table A.8: MAR (M odel 2 ): Predicted values mean errors and standard deviations for samples with different sizes discretized in p = 100 equidistant points based on 500 simulation replications.

			n + n 1 = 150
	Rate of missing data (%)
	Mean	12.912	28.026	45.472
	Median	13	28	45
	SD	3.524	4.493	5.118
	Criterion 3: [M SE × 10 3 ]	2.3556	2.9148	3.6204
		(1.6157)	(2.2111)	(2.7093)
	Criterion 4: [RT × 10 2 ]	6.0704	7.4692	9.2007
		(4.1999)	(5.6623)	(6.5708)
			n + n 1 = 450
	Rate of missing data (%)
	Mean	12.924	28.018	45.277
	Median	13	28	45.33
	SD	1.871	2.533	2.844
	Criterion 3: [M SE × 10 3 ]	0.8183	0.9882	1.2666
		(0.5391)	(0.6270)	(0.8146)
	Criterion 4: [RT × 10 2 ]	2.0977	2.5322	3.2364
		(1.3686)	(1.5836)	(2.0620)
			n + n 1 = 1800
	Rate of missing data (%)
	Mean	13.010	28.081	45.289
	Median	13	28.083	45.250
	SD	0.970	1.330	1.456
	Criterion 3: [M SE × 10 2 ]	0.1896	0.2360	0.2935
		(0.1216)	(0.1531)	(0.1812)
	Criterion 4: [RT × 10 2 ]	0.4856	0.6035	0.7492
		(0.3148)	(0.3959)	(0.4618)

  Table A.9: MAR (M odel 1 ): Estimation of θ mean square errors, variance and square bias for samples with different sizes discretized in p = 100 equidistant points based on 500 simulation replications. Table A.10: MAR (M odel 2 ): Estimation of θ mean square errors, variance and square bias for samples with different sizes discretized in p = 100 equidistant points based on 500 simulation replications.

			n + n 1 = 150 n + n 1 = 150
	Rate of missing data (%) Rate of missing data (%)
	Mean Mean	12.520 12.912	27.420 28.026	44.882 45.472
	Median Median	13 13	27 28	45 45
	SD SD	3.307 3.524	4.515 4.493	5.038 5.118
	Criterion 5: M SE × 10 2 Criterion 5: M SE × 10 2	20.33993 22.84329 25.59843 25.77594 30.94147 35.58789
	V ariance × 10 2 V ariance × 10 2	16.42143 17.02001 17.58919 17.87099 20.83862 21.5734
	Bias 2 × 10 2 Bias 2 × 10 2	3.918497 5.823277 8.009239 7.904949 10.10285 14.01449
			n + n 1 = 450 n + n 1 = 450
	Rate of missing data (%) Rate of missing data (%)
	Mean Mean	12.433 12.924	27.456 28.018	45.209 45.277
	Median Median	12.333 13	27.333 28	45.333 45.33
	SD SD	1.877 1.871	2.487 2.533	3.041 2.844
	Criterion 5: M SE × 10 2 Criterion 5: M SE × 10 2	8.923099 10.01299 12.37846 12.80462 14.15714 16.64587
	V ariance × 10 2 V ariance × 10 2	7.636041 8.680379 10.64885 6.696352 8.047992 10.44823
	Bias 2 × 10 2 Bias 2 × 10 2	1.287058 1.332613 1.729613 6.108267 6.109149 6.197638
			n + n 1 = 1800 n + n 1 = 1800
	Rate of missing data (%) Rate of missing data (%)
	Mean Mean	12.529 13.010	27.536 28.081	45.213 45.289
	Median Median	12.500 13	27.500 28.083	45.250 45.250
	SD SD	0.934 0.970	1.280 1.330	1.355 1.456
	Criterion 5: M SE × 10 2 Criterion 5: M SE × 10 2	3.268755 3.663376 4.294925 7.50709 8.091252 8.477034
	V ariance × 10 2 V ariance × 10 2	2.517848 2.870331 3.410527 1.746334 2.096911 2.495418
	Bias 2 × 10 2 Bias 2 × 10 2	0.7509066 0.793045 0.884398 5.760756 5.994341 5.981616

Proof of Lemma 5.3

The proof of this lemma is quite immediate, noticing that E (ΘΠ kn,obs X -ΘX ) 2 = E (Π kn,obs -I) X , θ 2 = (Π kn,obs -I) Γθ, θ

Proof of Theorem 2.2

From Theorem 2.1, the last term in the asymptotic development is negligible, 450 so we just have to achieve the usual trade-off between the square bias and the variance. Given that

we approximate this sum with the integral +∞ x L 2 ϕ(t) dt, which gives the desired result.

Proof of Theorem 2.3 455

First, if we follow the same lines of the proof of Lemmas 5.1 and 5.3 in Theorem 2.1 but with all the sample X 1 , . . . , X n , we get

and

Now, let us denote, for i = 1, . . . , n,

We immediately can write

Then, following the proof of Lemma 5.2 in Theorem 2.1, we denote


Mathieu Nespoulous
email: mathieu.nespoulous@univ-amu.fr

Renaud Denoyel

Mickaël Antoni

Microstructure Formation in Freezing Nanosuspension Droplets

ABSTRACT

The structural evolution of suspensions upon freezing is studied with optical microscopy in a suspended droplet configuration. Droplets have millimeter size and consist in an aqueous mixture of silica particles while the surroundings phase is hexane. Freeze-thaw cycles are applied to this system and a two-step freezing mechanism evidenced. A fast adiabatic growth of dendrites that invade the full droplets is first observed, and occurs within a few milliseconds. Then a slow process lasts for several seconds and corresponds to the release of solidification latent heat into the hexane phase. The striking feature of this work is to evidence that after the first freeze-thaw cycle flocculated microstructures are generated. When a second cycle is performed, microstructures further flocculate and generate, for dense silica suspensions, stable porous spheres of the size of the droplets. A phenomenological description based on repulsion or engulfment of particles by solidifying ice fronts is proposed.

TOC GRAPHICS

The solidification dynamics of undercooled liquids is of fundamental interest for many areas of science such as materials science [START_REF] Zhang | Aligned Two-and Three-dimensional Structures by Directional Freezing of Polymers and Nanoparticles[END_REF], metallurgy [START_REF] Li | Structural Evolution of Undercooled Ni-Cu Alloys[END_REF], cloud physics [START_REF] Bogdan | Formation of Mixed-phase Particles During the Freezing of Polar Stratospheric Ice Clouds[END_REF][START_REF] Gultepe | Fog Research: A Review of Past Achievements and Future Perspectives[END_REF] or food preservation [START_REF] Olivera | Effect of Freezing Rate in Textural and Rheological Characteristics of Frozen Cooked Organic Pasta[END_REF]. In the case of colloidal suspensions, it also has critical importance in engineering processes such as the so-called freeze-casting [START_REF] Deville | Metastable and Unstable Cellular Solidification of Colloidal Suspensions[END_REF][START_REF] Li | Freeze Casting of Porous Materials: Review of Critical Factors in Microstructures Evolution[END_REF] or freeze drying [START_REF] Han | Self-assembling Behavior of Cellulose Nanoparticles During Freeze-drying: Effect of Suspension Concentration, Particle Size, Crystal Structure, and Surface Charge[END_REF][START_REF] Vicent | Comparison of Freeze Drying and Spray Drying to Obtain Porous Nanostructured Granules from Nanosized Suspensions[END_REF]. In this process, the solid/liquid interface propagates through the colloidal suspension during solidification and its interaction with the particles can lead, after sublimation of the frozen phase, to the formation of porous materials [START_REF] Deville | Freezing as a Path to Build Complex Composites[END_REF][START_REF] Shanti | Particle Redistribution during Dendritic Solidification of Particle Suspensions[END_REF][START_REF] Thongprachan | Preparation of Macroporous Solid Foam from Multi-walled Carbon Nanotubes by Freeze-drying Technique[END_REF]. Considering the freezing of pure water, one particular aspect is the link between dendritic growth and latent heat release. It is known that dendrites can grow with velocities up to meters per seconds [START_REF] Hobbs | Ice Physics[END_REF]. This is far too fast to allow solidification latent heat to be evacuated. For undercooled droplets, two successive freezing regimes are therefore naturally expected and observed [START_REF] Hindmarsh | Experimental and Numerical Analysis of the Temperature Transition of a Suspended Freezing Water Droplet[END_REF][START_REF] Buttersack | Critical Radius of Supercooled Water Droplets: on the Transition toward Dendritic Freezing[END_REF]. Firstly, dendrites are formed. They can invade the whole medium but with a volume fraction limited by the undercooling temperature. The heat provided by the solidification is indeed reheating the droplets up to the melting point of water in a quasi-adiabatic process.

Secondly, driven by heat exchange with the cold surrounding medium, an ice front advances from the droplet interface to its center at a speed imposed by the heat conduction properties. As shown recently [START_REF] Buttersack | Critical Radius of Supercooled Water Droplets: on the Transition toward Dendritic Freezing[END_REF], the kinetics of this second regime depends upon the surface to volume ratio of the droplets and, for very small ones, the speed of the two regimes may become comparable.

In the present work, the objective is to investigate this problem when water is replaced by a colloidal suspension of silica particles. The following questions will be addressed: (i) are the aforementioned two regimes still observable, (ii) are the characteristic dynamic time scales modified and (iii) what is the incidence of freezing processes on the suspension stability and on its structural evolution?

To answer these fundamental questions, a system consisting in droplets pending at the end of capillary tip is used (see figure 1). Droplets are made of an aqueous silica nanoparticles suspension and immerged in hexane. When submitting this system to freeze-thaw cycles, an unexpected aggregation process where silica microstructures or even a full caking of the droplets can be evidenced. Scanning electron microscopy (SEM) experiments, performed on the freeze-dried nanofluid droplets after two cycles, reveal the formation of a porous spherical particle with size similar to that of the original droplet and exhibiting macropores. An explanation for this phenomenon based on the balance of the electrostatic interactions of the nanoparticles and ice surface charge is proposed. It is to be noted that a similar phenomenology has been evidenced when replacing hexane by air.

Hexane is mainly used for the optical purpose of refractive index matching with water.

When cooling pure water droplets below water melting temperature, the freezing occurs in two subsequent regimes at an undercooling temperature around -18°C in the present experiments. During the first regime, a dendritic growth phenomenon takes place as illustrated in figure 1a. It is fast since within a few milliseconds the droplet is completely filled with a three-dimensional network of entangled dendrites. As nucleation starts in most cases at the capillary tip, this network is growing top-bottom. This first regime will be denoted in the following as dendritic-freezing regime (DFR). The measured dendritic front velocity is about 0.15 m/s, in agreement with recent results [START_REF] Buttersack | Critical Radius of Supercooled Water Droplets: on the Transition toward Dendritic Freezing[END_REF]. The second regime, corresponds to the bulk-freezing regime (BFR) where a front of ice moving from the droplet interface to its center (see figure 1b) completely solidifies the droplet. It results from the quasi-equilibrium of isothermal bulk-water solidification and is mainly driven by the heat propagation in the entangled medium built up by the boundaries of the dendrites. In BFR, freezing front velocity is about 10 -4 m/s. This can be explained by the increase of the suspension viscosity when concentrated [START_REF] Azouni | Thermophysical Properties Effects on Segregation during Solidification[END_REF]. The data for v BFR are scattered but bounded within 5.0 10 -5 m/s and 1.5 10 -4 m/s. After the first freeze-thaw cycle, droplets recover their initial shape but surprisingly turn turbid. This is a typical signature of the presence of microstructures. When kept at room temperature for several hours after the first melting, a layer of sedimented solid material is observable at the bottom of the droplets. DLS analyses of the supernatant phase indicate the presence of particles with an average diameter of 120 nm and a dispersity parameter value of 0.25.

The formation of microstructures from the original nanoparticles suspension demonstrates that freezing triggers an irreversible aggregation process of the nanoparticles. When applying a second freeze-thaw cycle right after the first one, the same phenomenology shows up (first a DFR followed by a BFR). Unexpectedly, the final product, if silica concentration is 300 g/L, is a porous silica pearl with a diameter close to that of the original droplet, completely textured by water freezing. The mechanical strength of this pearl is actually low but high enough to resist the capillary pressure during drying at room temperature. Scanning electron microscopy analysis of the porous pearls of Figure 4 reveal a complex structure where nanoparticles define a three-dimensional macropore network. The morphology of the porous network seems to be different from what is obtained by freeze casting with a low unidirectional freezing rate [START_REF] Chino | Directionally Freeze-cast Titanium Foam with Aligned, Elongated Pores[END_REF][START_REF] Peppin | Experimental Verification of Morphological Instability in Freezing Aqueous Colloidal Suspensions[END_REF], which mainly leads to one dimensional structures. Here, the topology of the textures shows up as connected cells that resemble to the ones evidenced in cryo-SEM analysis of freeze-fractured nanofluid droplets [START_REF] Limage | Characterization of Solid-stabilized Water/Oil Emulsions by Scanning Electron Microscopy[END_REF]. The typical size of the visible pores in this micrograph is about 5 µm. They can be seen as signature of dendrite growth although their organization is not as regular as in directional freezing [START_REF] Deville | Freezing as a Path to Build Complex Composites[END_REF]. Beside this macroporosity, the walls are most probably mesoporous with a pore size of 8 nm given by the one of a compact arrangement of spherical nanoparticles. This was estimated from nitrogen adsorption at 77 K and application of the BJH method to a sample obtained by drying a large volume of the suspension. To understand the observed behavior, the interactions between the ice front and the nanoparticles have to be accounted for. Both van der Waals and electrostatic contributions must therefore be considered here [START_REF] Stahlberg | Electrostatic Interactions between a Charged Sphere and a Charged Planar Surface in an Electrolyte Solution[END_REF]. Four main interactions have to be accounted for in DFR and BFR: (i) The disjoining force between ice and nanoparticles which may be attractive or repulsive depending on the Hamaker constant. (ii) The electrostatic force between charged nanoparticles and charged ice front. This is a repulsive force in the present case since both silica and ice surface are negatively charged [START_REF] Kallay | Reversible Charging of the Ice-Water Interface. I. Measurement of the Surface Potential[END_REF]. (iii) The DLVO forces that operate between nanoparticles and create an osmotic pressure which acts against the increase of nanoparticle concentration when the ice invades the suspension.

(iv) Finally, the viscous forces that hinder the displacement of nanoparticles [START_REF] Azouni | Thermophysical Properties Effects on Segregation during Solidification[END_REF]. Buoyancy forces are negligible in this problem, moreover, freeze-thaw cycles are performed within short time scales (about 15 minutes), giving no time to sedimentation to set in. The surface potential of ice (point (ii)) can actually be split into two contributions. The first one is the reversible potential due to the ionization of OH groups at the surface of ice. It is considered to be around -25 mV in the experimental conditions of this work [START_REF] Kallay | Reversible Charging of the Ice-Water Interface. I. Measurement of the Surface Potential[END_REF][START_REF] Kallay | Reversible Charging of the Ice-Water Interface. II. Estimation of Equilibrium Parameters[END_REF]. The second is the freezing potential that develops at the liquid/solid interface for freezing water.

It appears when the propagation of liquid/solid fronts is fast and results from the fact that anions and cations are not engulfed at the same rate [START_REF] Workman | Electrical Phenomena Occurring During the Freezing of Dilute Aqueous Solutions and their Possible Relationship to Thunderstorm Electricity[END_REF]. This potential is known to decrease both with ionic strength and ice front speed [START_REF] Wilson | Workman-Reynolds Freezing Potential Measurements between Ice and Dilute Salt Solutions for Single Ice Crystal Faces[END_REF][START_REF] Wilson | Effect of Ice Growth Rate on the Measured Workman-Reynolds Freezing Potential between Ice and Dilute NaCl Solutions[END_REF].

Following the analysis of Azouni et al [START_REF] Azouni | Thermophysical Properties Effects on Segregation during Solidification[END_REF], a critical velocity, v C , can be introduced for the ice front: below this velocity, particles are repelled, above they are engulfed. v C is defined as the velocity at which the sum of the forces acting between particles and ice front is zero. Most of these forces are repulsive except the viscous force which is opposite to particle displacement. The equations of the model are detailed in the Supporting Information; and the values of the parameters corresponding to the silica suspension under focus here are shown in Table 1. The behavior of v C with respect to particle radius is displayed in Figure 5. One fundamental parameter in the model for the computation of v C is the distance h between the ice front and particles. Figure 5 displays the value of v C for h = 0.3 nm, 1.0 nm and 10 nm. Despite this broad range of h values, v C shows only a weak dependency with this parameter when r < 1 mm. This allows the definition of two distinct domains D 1 and D 2 (below and above the curve respectively). In D 1 particles are pushed along the propagating front while in D 2 they are engulfed [START_REF] Chino | Directionally Freeze-cast Titanium Foam with Aligned, Elongated Pores[END_REF].

Initially, the nanoparticles in the suspension have typical radius r = 15 nm. The critical velocity is then v C = 0.5 m/s and Figure 3 indicates that v C >v DFR >v BFR for the first freezethaw cycle. Both DFR and BFR are hence in D 1 and nanoparticles therefore repelled by the ice front. After DFR, the suspension concentration in the unfrozen part is increased by 20% [START_REF] Hindmarsh | Experimental and Numerical Analysis of the Temperature Transition of a Suspended Freezing Water Droplet[END_REF]. As a result, in the subsequent BFR, nanoparticles experience growing repulsive forces as they operate in more and more confined conditions. Dendritic network is indeed already fully developed and droplet volume most probably partitioned in ice free interstitial domains. As BFR goes on, these latter are reached by the freezing front and the suspension within them concentrates more and more until capillary pressure becomes comparable to osmotic pressure [START_REF] Olivera | Effect of Freezing Rate in Textural and Rheological Characteristics of Frozen Cooked Organic Pasta[END_REF]. Conditions are then reached to overcome the repulsive electrostatic barrier and flocculation to be triggered. The result here is the production of aggregates in the micrometric range which leads to a suspension with completely different granulometry.

After this first freeze-thaw cycle, the typical size of the particles, that are now microparticles, has been increased by two orders of magnitude. According to Figure 5, this means that v C drops down to about 10 -2 m/s. Figure 3 demonstrates that v DFR and v BFR are almost unchanged when performing a second cycle. This means that the mechanisms involved in particle flocculation can be studied along a horizontal line (resp. horizontal strip) located at 0.15 m/s for the new DFR (resp. (1.0 ± 0.5) 10 -4 for the new BFR) in figure 5. DFR situation is now drastically changed. As the suspension is micrometer sized, DFR occurs in D 2 and particles are, depending on their size, engulfed or partially engulfed by dendrites but not pushed as in the first cycle. In BFR however, v BFR < v C showing that engulfment is still not favored. Hence, similarly the first freeze-thaw cycle, microparticles are compressed against each other in a pre-existing complex frozen dendrite network.

However, conversely to what happened in this first cycle, the compression between microstructures now occurs in a 3-dimensional dendrite network that itself contains a relevant fraction of engulfed aggregates. These latter, as BFR takes place, are brought into contact with the microstructures already trapped in the dendrite free interstitial domains and create a new generation of flocculated objects. For large SiO 2 concentrations, this scenario yields a massing of the suspension that finally turns into a single porous pearl (see Figure 4).

In the same configuration of single, isolated droplet, methods are now being developed to produce large amounts of such pearls with the goal to carry out experimental campaigns of mercury porosimetry and gas adsorption. This well controlled protocol could also be seen here as new and original route for the synthesis of porous materials. This article focuses on the freezing dynamics of silica nanoparticle suspensions in a suspended droplet configuration. A two-step scenario is evidenced where the solidification first exhibits a fast dendritic growth followed by bulk freezing. Solidification velocity is measured and shown to depend upon the nanoparticle concentration. When performing freeze-thaw cycles the structure of the dispersion is strongly modified since from nanodispersed it turns microdispersed. To explain this scale change, the interplay between surface charge of solidification fronts in the liquid/solid coexisting system and interparticle repulsive forces is addressed. Flocculation mechanisms generated by particle confinement in a dendrite network are shown to occur. After two freeze-thaw cycles of the droplets, initially dense suspensions are massing and form a single flocculated pearl. This work demonstrates that granulometry of suspensions strongly evolve under freezing and raises the fundamental question of the consequence of this phenomenon in the microphysics of atmospheric clouds when air contains fine particles. Density of water @ 10°C r L (kg.m -3 ) 1000

Density of water @ 0°C r S (kg.m -3 ) 917

Density of np r P (kg.m -3 ) 2648

Viscosity of water @ 0°C µ (Pa.s) 1.79x10 -3 [START_REF] Lide | Handbook of Chemistry and Physics[END_REF] Mass concentration of np c (kg.m -3 ) 0.3

Debye length l D (m) 9.1x10 -9 Inverse Debye length k = 1/l D (m -1 ) 1.1x10 8 Surface potential Y 0 (V) 0.03 np-water-np Hamaker const. A p-p (J) 7.46x10 -21 [START_REF] Wilen | Dispersion-Force Effects in Interfacial Premelting of Ice[END_REF] np-water(ice) Hamaker const. A p-i (J) -8.2x10 -22 [START_REF] Azouni | Thermophysical Properties Effects on Segregation during Solidification[END_REF] *np stands for silica nanoparticle.

Experimental Methods

Sample preparation -Aqueous solutions of amorphous silicon oxide (SiO 2 ) nanoparticles (Levasil 200, AkzoNobel) were used. The SiO 2 nanoparticles have a diameter of 30 nm.

They are negatively charged and show a good stability when dispersed into water (with a zeta potential of about -40 mV). Stock solutions of 30% mass fraction in colloids are used.

A microliter droplet of the colloidal dispersion is formed within a vial previously filled with hexane. The alkane acts as a continuous surrounding phase for the pending droplet.

Freezing experiment -The pending droplet solidification dynamics is analyzed with optical microscopy. A Peltier cooling system allows to cycle temperature at a given cooling rates (usually 1°C/mn). The microscope is equipped with a CMOS camera (Mikrotron MC1310). The typical freeze-thaw cycle applied on the droplet is: from room temperature to -20 °C to room temperature, at a constant 1 °C/mn rate.

SEM experiment -A previously described procedure was followed to analyze the pearls Electric constant of vacuum e 0 (C 2 .J -1 ) 8.85x10 -12 [START_REF] Lide | Handbook of Chemistry and Physics[END_REF] Dielectric constant of water @ 0°C e W (C 2 .J -1 ) 87.9 [START_REF] Lide | Handbook of Chemistry and Physics[END_REF] Surface charge of np-water s S (C.m -2 ) -0.02 [START_REF] Barisik | Size Dependent Surface Charge Properties of Silica Nanoparticles[END_REF] Surface charge of ice-liq. interface s P (C.m -2 ) -0.01 [START_REF] Kallay | Reversible Charging of the Ice-Water Interface. I. Measurement of the Surface Potential[END_REF]

Figure 1 .

 1 Figure 1. a) Optical microscopy image of DFR of a pure water droplet (cooling rate

Figure 2 .

 2 Figure 2. Optical micrograph of suspended nanofluid droplets in DFR (cooling rate

Figure 3 .

 3 Figure 3. (a) (resp. (b)) v DFR (resp. v BFR ) as a function of [SiO2]. Diamonds (resp.

Figure 4 .

 4 Figure 4. a) Scanning electron microscopy pictures of a freeze-dried droplet ([SiO2] = 300

Figure 5 .

 5 Figure 5. v C as a function of particle radius. In domain D 1 , v<v C , and particles are pushed

Table 1 :

 1 Parameters and constants used in the model.

	Description	Parameter (unit)	Value [ref.]
	Radius of np*	R (m)	1.5x10 -8
	np -ice front distance	h 0 (m)	1x10 -9
	Boltzman constant	k B (J.K -1 )	1.38x10 -23
	Gravity constant	g (m.s -2 )	9.81
	Working temperature	T (K)	273.15
	Curvature parameter	a	0 [16]

 [START_REF] Limage | Characterization of Solid-stabilized Water/Oil Emulsions by Scanning Electron Microscopy[END_REF]. The scanning electron microscope used is a Philips XL30 SFEG STEM associated to an Oxford EDS analysis system. Supporting Information. Detailed expression of the forces used in the model to compute the critical velocity.


Elise Buisson
email: elise.buisson@univ.u-3mrs.fr

Thierry Dutoit

Franck Torre

Christine Römermann

Peter Poschlod

Calamintha * Conyza

Filago Pyramidata Memi

Rumex * Cynoglossum

Thymus

Trigonella

Logfia * Linaria

The implications of seed rain and seed bank patterns for plant succession at the edges of abandoned fields in Mediterranean landscapes

Keywords: ants, dry grassland, mosaic of herbaceous communities, seed bank, seed rain, steppe

Some environmental variables and above-ground vegetation, seed rain, ant-borne seeds and seed banks were studied on three abandoned fields, at the margins between these fields and a remnant patch of a native steppe vegetation of a high value for nature and species conservation. While the fields were the same size, were adjacent to the same patch of remnant steppe and were cultivated with the same crop, site-specific environmental variables contributed to 23% of the vegetation patterns; each site was characterised by its unique historical trajectory and thus, by a particular set of species. Distance from boundaries contributed to 10% of the vegetation patterns. Species characterising the surrounding steppe were found close to boundaries; species characterising abandoned fields were found further away. Winter seed banks and summer deep seed bank did not contribute much to either effect and were characterised by species dating back from past cultivation. Conversely, summer surface seed bank greatly contributed to (83%) the differences in species composition between the three fields. Seed rain contributed to differences in species composition between fields (91%) and distance (76%). Ant-borne seeds largely contributed to the differences between fields (87%). The colonisation of steppe species on field margins occurs mainly through seed rain and is very slow and incomplete. In a semi-arid fragmented openlandscape, patches of native vegetation do not play a great role in colonisation processes, and itinerant sheep grazing is insufficient to initiate recovery.

Introduction

While many studies on grazing, mowing and burning regimes have been carried out to plan conservation management for individual grasslands [START_REF] Hillier | Calcareous grassland: ecology and management[END_REF][START_REF] Poschlod | The historical and socio-economic perspective of calcareous grasslands-lessons from the distant and recent past[END_REF], only few studies plan for long-term management within fragmented open-landscapes and studies on secondary succession and edges have yet to be relevant to long-term management. Studies focusing on edges assessed either the influence of agriculture-related factors on field-margin vegetation [START_REF] Kleijn | Patterns in species composition of arable field boundary vegetation[END_REF] or the influence of field-margin vegetation on crops [START_REF] Dutoit | Field boundary effects on soil seed banks and weed vegetation distribution in an arable field without weed control (Vaucluse, France)[END_REF][START_REF] Von Arx | Land-use intensity and border structures as determinants of vegetation diversity in an agricultural area[END_REF].

Investigating colonisation processes at edges of abandoned fields is important [START_REF] Wilson | The distribution of dicotyledonous arable weeds in relation to distance from the field edge[END_REF] because the field potential to recover via long-term seed bank or longdistance seed rain has often been found to be relatively low in Northern Europe [START_REF] Graham | A field investigation of germination from the seed bank of a chalk grassland ley established on former arable land[END_REF]Hutchings and Booth, 1996a;[START_REF] Kalamees | The role of the seed bank in gap regeneration in a calcareous grassland community[END_REF]. Studying edges is even more important in dry climates, where drought induces conditions that dramatically slow down successional processes [START_REF] Blondel | Biology and Wildlife of the Mediterranean Region[END_REF]. Results may have important conservation implications, particularly in the Mediterranean basin where grasslands harbour a high proportion of endangered taxa [START_REF] Médail | Biodiversity hotspots in the Mediterranean Basin: setting global conservation priorities[END_REF].

The hypothesis of this research was that remnant patches of vegetation are a source of seeds for the colonisation of abandoned fields either by seed rain or ant dispersal. The aims of the present study were (1) to describe the colonisation processes on abandoned field margins over a complete growing season in order to determine how steppe vegetation reestablishes at the margins between abandoned fields and a remnant patch of steppe in order (2) to provide information for the long-term management of abandoned fields.

Methods

Study area

The study was carried out in the plain of La Crau, located ~50 km north-west of Marseille in south-eastern France. The steppe vegetation of the plain evolved with: 1) a dry Mediterranean climate; 2) a shallow soil (<400 mm) and an impermeable bedrock; 3) itinerant sheep grazing since 3000 BP [START_REF] Buisson | Colonisation by native species of abandoned farmland adjacent to a remnant patch of Mediterranean steppe[END_REF]. The steppe vegetation is dominated by the stress-tolerant species Brachypodium retusum and Thymus vulgaris, which represent ~50% of the biomass, and is composed of a great diversity of annual species (Buisson and Dutoit, in press). The remnant central patch of steppe was fragmented by row crop and melon cultivations between 1965 and 1985, after which all patches of steppe and fields were grazed by itinerant sheep from February to June.

Three abandoned fields (A, B, C) were selected on the sheepfold of Peau de Meau (43°33'E 4°50'N, elevation 10m), in the centre of the plain. All fields were adjacent to one same remnant patch of undisturbed steppe (D) in order to avoid confounding changes in the floristic composition on the margins of abandoned fields with differences due to species composition of steppe patches [START_REF] Römermann | Influence of former cultivation on the unique Mediterranean steppe of France and consequences for conservation management[END_REF]. However, fields had varied cultivation periods, dates of abandonment and locations given in Table 1.

Sampling

On each studied field (A, B, C), three transects (1, 2, 3) were set out perpendicular to the field boundary. Each transect started at the boundary and stretched towards the centre of abandoned fields and measured 10 m because the visible gradient of the typical steppe perennial species T. vulgaris % cover never exceeded 10 m. The three transects per field were laid 10 m apart because most herbaceous species do not disperse by wind over more than a few meters [START_REF] Verkaar | On the ecology of short lived forbs in chalk grasslands: dispersal of seeds[END_REF], especially in the plain were vegetation grow very low. Collection points were set 1 m apart along the transects starting at the field edge (point 0) between steppe and fields (0 to 10 m; 11 samples). All data were collected in 2001, a year with an average climate.

Above-ground vegetation was recorded in 10 sub-quadrats of 40  40 cm at each sampling point totalling 990 sub-quadrats, in May 2001. T. vulgaris and B. retusum % cover, as well as stones, vegetation and bare ground % cover were visually estimated by one observer in each quadrat. All other species were sampled in each quadrat using the presence/absence method. To avoid disturbing the vegetation, other sampling was taken on the other side of the transects.

200g of soil was taken at sampling points 0, 5 and 10 m along each transect in February 2001 (27 samples total). Samples were dried and sieved through a 200 µm mesh sieve. Total nitrogen was measured using the Kjeldahl method, phosphorus using the Olsen method and total carbon using the Anne method [START_REF] Baize | Guide des analyses en pédologie[END_REF]. The concentrations of calcium, potassium and magnesium available to plants were measured using the Shollenberger and Dreibelbis method [START_REF] Baize | Guide des analyses en pédologie[END_REF]. The soil pH was also measured in H2O [START_REF] Baize | Guide des analyses en pédologie[END_REF].

The persistent seed bank, seeds that persist in the seed bank >1year, [START_REF] Thompson | The soil seed bank of North West Europe: methodology, density and longevity[END_REF] was sampled in February 2001 before the input of fresh seeds and after the germination of most seeds in the fall and winter. To deduce the transient seed bank, the summer seed bank, seeds that persist in the seed bank <1year, [START_REF] Thompson | The soil seed bank of North West Europe: methodology, density and longevity[END_REF], was sampled at the beginning of July 2001 after seeds were shed. Cores could not be taken because the soil was extremely rocky. For each sampling point, two 500 cm³ containers were filled, one for the 0-10 cm depth and one for the >10 cm depth (198 samples for summer and 198 for winter seed bank). The volume of soil sampled was 1500 cm³/soil layer/field/distance from the boundary; the volume recommended by [START_REF] Bakker | Seed banks and seed dispersal: important topics in restoration ecology[END_REF] was 1200 cm³. The 396 samples of soil seed bank were set to germinate using the standard procedure recommended by Ter [START_REF] Ter Heerdt | An improved method for seed-bank analysis: seedling emergence after removing the soil by sieving[END_REF]. Soil samples were washed with water on sieves of 4 mm and 200 µm mesh sizes to reduce bulk and clay. The concentrated soil was spread in a thin layer in trays filled with 20mm of vermiculite and topped with medical compresses of 100 µm mesh size. Trays were watered every day from below. Emerged seedlings were identified, counted and removed weekly.

Seed rain was sampled from April 15 th to October 15 th 2001. Sticky traps, made of paper filters coated with grease and placed in 140 mm diameter Petri dishes (Hutchings and Booth, 1996a), were fixed with pegs at all sampling points along the transects. Sticky papers were changed every two weeks (1287 traps total). Seeds were identified to the species or genus level under a dissecting microscope using a reference collection.

To assess ant seed dispersal from the steppe to the abandoned fields, ant pitfall traps were set out at all sampling points along the transects. These traps, consisting of 50 mm diameter  100 mm high containers, were filled with a non-attractive liquid of 50% glycol to conserve seeds and reduce evaporation, 49% water, 1% acetic acid to act as an anti-fungic and a few drops of tensioactive detergent to make the ants sink. Traps were changed every two weeks from April 15 th to October 15 th 2001 (1287 traps total). Seeds were identified to the species or genus level and harvesting ants (Messor barbarus and Messor sanctus) were identified and counted.

Data analyses

Multivariate analyses were used to study simultaneously the complex relationships among species, and between the species and their environment. A partial canonical correspondence analysis pCCA [START_REF] Ter Braak | The analysis of vegetation-environment relationships by canonical correspondence analysis[END_REF] was used to determine 1) the relationship between the above-ground vegetation and a set of environmental variables once the impact of distance variables was statistically removed; and 2) the relationship above-ground vegetation and the distance from boundaries once the impact of the environmental variables was statistically removed [START_REF] Mcintyre | How environmental and disturbance factors influence species composition in temperate Australian grasslands[END_REF]. The above-ground vegetation data and some environmental variables (pH, P2O5, C:N, K2O and % cover of vegetation) were used at points 0, 5 and 10 m from boundary.

To identify which of the six data-sets best contributed to the composition of the above-ground vegetation, a K-sets principal component analysis was used (K-sets PCA, Pagès and Escofier,1994). Classical multivariate approaches consider alternatively: i) as many separate PCAs as there are data-sets, but the interpretation is doubtful as it is based on empirical comparisons of separate analyses; ii) a unique overall PCA of all six data-sets, although some may weigh more than others in the analysis. A K-sets PCA allowed us to analyse all six data-sets by weighting differently each variable so that they played a balanced role in the analysis. Separate PCAs were first conducted for each data-set and then the inverse of the first eigenvalue was used as a weighting coefficient for the global K-PCA analysis [START_REF] Pagès | Multiple factor analysis (AFMULT package)[END_REF]. Computations for the pCCA and the K-sets PCA were performed using ADE4 software [START_REF] Thioulouse | ADE-4: a multivariate analysis and graphical display software[END_REF] To explain plant distribution at field margins, the means per distance (n=9 transects) were calculated for (1) soil variables (2) total and steppe-only species richness and density of above-ground vegetation, seed rain, seeds in ant-traps, summer and winter surface and deep seed banks. Data were directly compared or log10 transformed [START_REF] Sokal | Biometry: the principles and practice of statistics in biological research[END_REF] before being compared with one-way ANOVAs (Statistica software, version 6.0), followed by post-hoc tests (Least Significant Difference). Non-parametric multiple comparisons by STP tests (Simultaneous Test Procedure, based on Mann-Whitney U test, [START_REF] Sokal | Biometry: the principles and practice of statistics in biological research[END_REF] were performed on non-normal data.

Results

The total number of species found in the vegetation, seed bank, seed rain and ant-traps are given for each field in Table 2. All the species found during this study are listed in the Appendix and their frequency of occurrence in each of the various types of data are recorded.

The partial CCA showed that the environmental variables alone explained 23% of the variance of the above-ground vegetation matrix (no figure). The partial CCA also showed that distances from sampling points to boundary contributed to 10% of the above-ground vegetation matrix (Fig. 1). This distance effect showed that points close to the boundary were associated with species such as Asphodelus ayardii, B. retusum, Filago pyramidata, Linaria arvensis and T. vulgaris, whereas points far from the boundary were associated with species such as Calamintha nepeta, Lepidium graminifolium, Rumex pulcher and Cynoglossum officinale (Fig. 1). Non-parametric multiple comparisons by STP tests confirmed that T. vulgaris % cover was significantly higher at 0, 1 and 2 m from boundaries than at 10 m (U=34.7 p<0.001); and phosphorus contents significantly increased with the distance from boundaries (U=11.7 p<0.001).

Axis 1 of the K-sets PCA was correlated with the first principal components of all seven matrices and particularly with above-ground vegetation (89%), seed rain (91%), antborne seed (87%) and summer surface seed bank (83%) matrices (Fig. 2), showing correlation between these variables. This axis (6%) separated abandoned field A, associated with species, such as Bromus sp., Diplotaxis tenuifolia, Polycarpon tetraphyllum and Vulpia sp., from abandoned field C points associated with species, such as Aegilops sp., Bellis sylvestris, Carthamus lanatus and Senecio vulgaris. Abandoned field B points were grouped between A and C. A site effect was thus observable not only on above-ground vegetation but also on most of the types of data studied except on winter seed bank. Axis 2 of the K-sets PCA was particularly correlated with the second components of above-ground vegetation (80%) and seed rain (76%) matrices (Fig. 2). This axis (4%) separated the points close to the boundary from the points far from the boundary. Points close to the boundary were associated with species, such as Erodium cicutarium, Euphorbia exigua, Galium sp. and P. tetraphyllum for both matrices as well as Linum trigynum, Sideritis romana and T vulgaris in the vegetation and Brachyposium distachyon and Plantago bellardii in the seed rain (Fig. 3). Points far from the boundary were associated with species, such as Conyza sp. in the vegetation and Lobularia maritima in the seed rain (Fig. 3).

ANOVA showed that steppe species richness of above-ground vegetation and steppe species richness of seed rain were significantly higher at 0, 1, 2, and 3 m as well as 4 m for vegetation than at 7 to 10 m from boundaries (respectively F=9.5 and F=6.1, P0.001 and LSD test P<0.05).

Discussion

Independently from the site effect and from the differences in species number found on each field (Table 2) which are well-known phenomena [START_REF] Buisson | Colonisation by native species of abandoned farmland adjacent to a remnant patch of Mediterranean steppe[END_REF]), original patterns were observed which were comparable on the three field margins. The composition of above-ground vegetation changed with distance from boundary: there were more steppe species close to boundaries and more arable weed species and mesophilous species further away. Soil pH and soil concentrations in phosphorus increased with distance from boundary because more fertiliser was applied on the fields than on the field margins that were less deeply ploughed. [START_REF] Grime | Floristic diversity in a model system using experimental microcosms[END_REF] and [START_REF] Gough | A comparison of soil fertility between semi-natural and agricultural plant communities: implications for the creation of species-rich grasslands on abandoned agricultural land[END_REF] have shown that ruderal species out-compete dry grassland species when phosphorus is in excess. Although a phosphorus gradient existed on the field margins in this study, [START_REF] Buisson | Colonisation by native species of abandoned farmland adjacent to a remnant patch of Mediterranean steppe[END_REF] have shown that concentrations observed in La Crau were not high enough on field margins (mean phosphorus = 0.019g/kg) to explain the degree of vegetation change [START_REF] Janssens | Relationship between soil chemical factors and grassland diversity[END_REF].

On field margins, the deep seed bank mainly consisted of arable weed species (Kickxia elatine and Portulaca oleracea), ruderal and nitrophilous species (Calamintha nepeta, L. maritima, Polygonum aviculare, S. vulgaris, Solanum nigrum) dating back to melon cultivation. The surface seed bank consisted not only of recently established steppe species (seed bank <1 year), but also of the species cited above, some of which (e.g. P. oleracea and S. nigrum) have not been found in above ground vegetation since the abandonment of melon cultivation because favourable conditions no longer exist (ploughing, fertiliser, irrigation etc.). Re-establishment of steppe species from seed bank appeared to be minimal, as most steppe species were shown to have a transient seed bank (see also [START_REF] Römermann | Influence of former cultivation on the unique Mediterranean steppe of France and consequences for conservation management[END_REF], although margins offered better environmental conditions for stress-tolerant steppe species than the middle of the fields.

The stress-tolerant steppe species now growing on field margins, but absent in the seed bank, came from seed rain or ant-borne seeds. Seed dispersal by ants on the margins, while worth noting (about 10000 seeds and 10000 ants trapped throughout the study), was probably limited to small distances, as observed by [START_REF] Harrington | The effect of fire and ants on the seed-bank of a shrub in a semi-arid grassland[END_REF]. Seed rain thus contributed to most steppe species colonisation. Even sampled over a complete fructification and dispersal season, seed rain and seed dispersal by ants did not suffice for the re-establishment of steppe species on the fields since only a few steppe species were found on margins while typical steppe species such as Stipa capillata, A. ayardii, Plantago holosteum, Hyssopus canescens, and Fumana procumbens were not. These results contrast strongly with Etienne et al.'s (1998) assessment of potential steppe vegetation reestablishment in La Crau and with results found in north-western Europe [START_REF] Gibson | Grazing and vegetation change: deflected or modified succession[END_REF].

The results of this study under Mediterranean climate show that the role of the seed bank in the re-establishment of former steppe on abandoned fields is minimal and thus corroborate those found elsewhere [START_REF] Graham | A field investigation of germination from the seed bank of a chalk grassland ley established on former arable land[END_REF]Hutchings and Booth, 1996a). While in many landscapes, remnant patches of native vegetation have a positive role on the secondary succession of abandoned fields, in a Mediterranean fragmented openlandscape, remnant patches of steppe adjacent to abandoned fields do not seem to have much of a role in the colonisation processes of steppe species as seed dispersal by wind or ants is limited.

Plant and insect dispersal in fragmented landscape has been studied in open landscapes and it has been concluded that livestock [START_REF] Poschlod | Plant species richness in calcareous grasslands as affected by dispersability in space and time[END_REF] are an excellent vector to disperse seeds between plots by epizoochory and endozoochory [START_REF] Fischer | Experimental studies on the dispersal of plants and animals on sheep in calcareous grasslands[END_REF] in Northern Europe. The steppe of La Crau and abandoned fields are grazed by itinerant sheep and their daily ranging patterns could easily be planned so that they graze on the steppe first, then on fields, in order to disperse steppe-species seeds and insects onto fields. However, grasshoppers thrive and most plants set seeds in summer [START_REF] Bourrelly | Dynamique annuelle et production primaire nette de l'écosystème steppique de Crau[END_REF]) when all flocks are in transhumance in the Alps. Also, [START_REF] Fischer | Experimental studies on the dispersal of plants and animals on sheep in calcareous grasslands[END_REF] have shown that most of the seeds dispersed came from plants bigger than 80 cm high and [START_REF] Dutoit | Rôles des ovins dans le transport de graines d'espèces messicoles: le cas d'une exploitation agricole du Parc Naturel Régional du Luberon[END_REF] have shown that, in a Mediterranean area where sheep graze year around, sheep were not a good vector of rare arable weeds. In La Crau, most dicotyledons grow rather small and have small and smooth seeds resembling those of rare arable weeds. Therefore, although sheep grazing is required to maintain steppe vegetation [START_REF] Fabre | Intégration des mesures agri-environnementales aux systèmes d'élevage ovin de la Crau[END_REF], it is not sufficient for the long-term management of abandoned fields.

The long-term conservation of the steppe has to include a plan to re-connect the larger patches of steppe. Future experiments should focus on sowing seeds or transplanting nurse species. This management tool has been used to enhance species richness of grasslands in northern Europe. Its success depends on the seeds sown (Hutchings and Booth, 1996b), on facilitation processes between transplanted nurse species and others native species [START_REF] Pywell | Restoration of species-rich grassland on arable land: assessing the limiting processes using a multi-site experiment[END_REF] or on concurrent vegetation management [START_REF] Warren | Effects of sowing and management on vegetation succession during grassland habitat restoration[END_REF].

Acknowledgements

This study was financed with a grant from "Ministère du Développement Durable et de l'Ecologie", program 'Espaces Protégés'. We acknowledge the CEEP Ecomusée de Crau, Mr Bernhardt, Fadda, Gerbaud, and Mrs Bottone, Rolando, and Trivelly for technical support and field assistance, as well as Mr Hayes and Mrs Sweetko for proof reading the manuscript. We also thank several anonymous reviewers for helpful comments on the manuscript.

Appendix

Total frequency of occurrence (%) of species in all three fields in the vegetation (990 quadrats of 40  40 cm), in the 1287 seed rain traps, in the 1287 ant traps and in the various seed banks Nomenclature according to [START_REF] Kerguélen | Index synonymique de la flore de France[END_REF].

Species

Typical steppe vegetation [START_REF] Römermann | Influence of former cultivation on the unique Mediterranean steppe of France and consequences for conservation management[END_REF] Above-ground vegetation

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. Partial CCA performed on distance and vegetation variables once the impact of environmental variables was statistically removed. All sampling points are shown. Point abbreviations: first number = distance to boundary (0 to 10 m); letter = field A, B or C; second number = transect number (1, 2 or 3). Only the species with high contributions are shown. The nine sampling points far from the boundaries are circled.

Table 1 .

 1 Cultivation history of the three abandoned fields A, B, C and the steppe on the sheepfold of Peau de Meau. After abandonment, grazing occurred on all fields, as on the steppe, by itinerant sheep from February to June.

	Seed rain	Ant-borne-seeds	Summer surface seed bank	Summer deep seed bank	Winter surface seed bank	Winter deep seed bank