Keywords: integer quadratic knapsack problem, separable objective function, direct expansion, binary expansion, piecewise interpolation

linear interpolation can not be applied to transform (QMKP) into an equivalent 0-1 piecewise linear problem.

Rewriting integer variables into zero-one variables: some guidelines for the integer quadratic multi-knapsack problem

Abstract: This paper is concerned with the integer quadratic multidimensional knapsack problem (QMKP) where the objective function is separable. Our objective is to determine which expansion technique of the integer variables is the most appropriate to solve (QMKP) to optimality using the upper bound method proposed by [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. To the best of our knowledge the upper bound method previously mentioned is the most effective method in the literature concerning (QMKP). This bound is computed by transforming the initial quadratic problem into a 0-1 equivalent piecewise linear formulation and then by establishing the surrogate problem associated. The linearization method consists in using a direct expansion initially suggested by [START_REF] Glover | Improved linear programming formulations of nonlinear integer problems[END_REF] of the integer variables and in applying a piecewise interpolation to the separable objective function. As the direct expansion results in an increase of the size of the problem, other expansions techniques may be utilized to reduce the number of 0-1 variables so as to make easier the solution to the linearized problem. We will compare theoretically the use in the upper bound process of the direct expansion (I) employed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF] with two other basic expansions, namely: (II) a direct expansion with additional constraints and (III) a binary expansion. We show that expansion (II) provides a bound which value is equal to the one computed by [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. Conversely, we provide the proof of the non applicability of expansion (III) in the upper bound method. More specifically, we will show that if (III) is used to rewrite the integer variables into 0-1 variables then a

Introduction

This paper deals with the integer quadratic multi-knapsack problem (QMKP) where the objective function is separable. Problems of this structure arise in numerous industrial and economic situations, for instance in production planning [START_REF] Stecke | Nonlinear integer production planning problems[END_REF], reliability allocation [START_REF] Nakagawa | Optimal reliability allocation by branch and bound techniques[END_REF] and finance [START_REF] Faaland | An integer programming algorithm for portfolio selection[END_REF]. These include the main application of (QMKP) which is in the portfolio management area where the investments are independent, see [START_REF] Djerdjour | A surrogate-based algorithm for the general quadratic multidimensional knapsack[END_REF] and [START_REF] Faaland | An integer programming algorithm for portfolio selection[END_REF]. Nevertheless, solving (QMKP) efficiently will constitute a starting point to solve the more general and realistic portfolio management problem where the investments are dependent, i.e. the objective function is non separable.

The integer quadratic multi-knapsack problem (QMKP) where the objective function is separable consists in maximizing a concave separable quadratic integer function subject to m linear capacity constraints. It may be stated mathematically as follows:

( ) ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ = = ≤ ≤ = ≤ = - = ∑ ∑ ∑ = = = n j x n j u x m i b x a x f x d x c x f QMKP j j j n j i j ij n j n j j j j j j j , , 1 integer, , , 1 , 0 , , 1 , t. s.

) ( max ) (

where the coefficients c j , d j , a ij , b i are nonnegative. The bounds u j of variables x j are pure integers, with ( )

j j j d c u 2 / ≤
. Indeed, the separable objective function is concave which implies that for all function f j , ( )

j j j d c x 2 / * ≤
, where x * j is the optimal solution of the program ( )

j j x x f j 0 max ≥ .
The problem (QMKP) which is a NP-hard problem [START_REF] Chaillou | Best network flow bounds for the quadratic knapsack problem[END_REF] is a generalization of both the integer quadratic knapsack problem [START_REF] Bretthauer | The nonlinear knapsack problem -algorithms and applications[END_REF] and the 0-1 quadratic knapsack problem where the objective function is subject to only one constraint [START_REF] Billionnet | A new upper bound for the 0-1 quadratic knapsack problem[END_REF].

Since, (QMKP) is NP-hard, one should not expect to find a polynomial time algorithm for solving it exactly. Hence, we are usually interested in developing branch-andbound algorithms. A key step in designing an effective exact solution method for such a maximization problem is to establish a tight upper bound on the optimal value. Basically, the available upper bound procedures for (QMKP) may be classified as attempting either to solve efficiently the LP-relaxation of (QMKP) (see [START_REF] Bretthauer | The nonlinear knapsack problem -algorithms and applications[END_REF] and [START_REF] Körner | Integer quadratic programming[END_REF]) or to find a good upper bound, of better quality than the LP-relaxation of (QMKP), transforming (QMKP) into a 0-1 linear problem easier to solve (see [START_REF] Djerdjour | A surrogate-based algorithm for the general quadratic multidimensional knapsack[END_REF] and [START_REF] Mathur | A branch and search algorithm for a class of nonlinear knapsack problems[END_REF]). To the best of our knowledge, the upper bound method we have proposed in a previous work [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF] is better than the existing methods (Djerdjour, Mathur and Salkin algorithm [START_REF] Djerdjour | A surrogate-based algorithm for the general quadratic multidimensional knapsack[END_REF], a 0-1 linearization method, a classical LP-relaxation of (QMKP)) from both a qualitative and a computational standpoint. We have first used a direct expansion of the integer variables, originally suggested by Glover [START_REF] Glover | Improved linear programming formulations of nonlinear integer problems[END_REF], and apply a piecewise interpolation to the objective function: an equivalent 0-1 linear problem is thus obtained. The second step of the algorithm consists in establishing and solving the surrogate relaxation problem associated to the equivalent linearized formulation.

Nevertheless the transformed linear formulation encounters numerous 0-1 variables because of the direct expansion used (denoted by expansion (I) in the following).

Consequently, other expansions techniques may be utilized to reduce the number of 0-1 variables so as to make easier the solution to the linearized problem. Let us consider the three basic expansions for rewriting the integer variables of (QMKP) into 0-1 variables:

• Expansion (I): direct expansion

∑ = = j u k j jk x y 1 , { } j j u x K 1 ∈ , { } 1 , 0 ∈ jk y n j K 1 = ∀ , j u k K 1 = • Expansion (II): direct expansion with additional constraint ∑ = = j u k j jk x y k 1 ' , ∑ = ≤ j u k jk y 1 ' 1 { } j j u x K 1 ∈ , { } 1 , 0 ' ∈ jk y n j K 1 = ∀ , j u k K 1 = • Expansion (III): binary expansion. ⎡ ⎤ ∑ = = ) log( 1 2 j u k j jk k x z , ⎡ ⎤ ∑ = ≤ ) log( 1 2 j u k j jk k u z { } j j u x K 1 ∈ , { } 1 , 0 ∈ jk z n j K 1 = ∀ , ⎡ ⎤ ) log( 1 j u k K =
The purpose of this note is to evaluate the impact of the use of the above expansion techniques, on the upper bound computation developed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. More specifically, we will determinate which expansion is the most appropriate to be used in the upper bound method for (QMKP) [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. We will compare theoretically the use of the direct expansion (I) with the direct expansion with additional constraints (II) and with the binary expansion (III). We will show that the use of (II) provides a bound which value is equal to the one computed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. Conversely, we provide the proof of the non applicability of both (III) and a linear interpolation to transform (QMKP) into an equivalent 0-1 piecewise linear problem.

The paper is organized as follows. The next section summarizes the upper bound method developed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF] detailing the direct expansion (I) of integer variables and the piecewise interpolation. In Section 3, the direct expansion with additional constraints (II) is applied to (QMKP) so as to compute the upper bound suggested by Quadri et al. [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. Section 4 is dedicated to the binary expansion (III). We finally conclude in Section 5.

In the remainder of this paper, we adopt the following notations: letting (P) be an integer or a 0-1 program, we will denote by ( P ) the continuous relaxation problem of (P). We let Z[P] be the optimal value of the problem (P) and Z[ P ] the optimal value of ( P ). Finally ⎡ ⎤

x (resp. ⎣ ⎦ x ) will denote the smallest (resp. highest) integer greater (resp. lower) than or equal to x.

Section 2. Direct expansion of the integer variables

In this section we summarize the upper bound method for (QMKP) proposed by Quadri et al. [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. First, an equivalent formulation is obtained by using a direct expansion (I) of the integer variables x j as originally proposed by Glover [START_REF] Glover | Improved linear programming formulations of nonlinear integer problems[END_REF] and by applying a piecewise interpolation to the initial objective function as discussed in [START_REF] Djerdjour | A surrogate-based algorithm for the general quadratic multidimensional knapsack[END_REF].

The direct expansion of the integer variables x j consists in replacing each variables x j by a sum of u j 0-1 variables y jk such that ∑ = =

( ) ( ) { } ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ = = ∈ = ≤ = ∑ ∑ ∑ ∑ = = = = 1 1 , 1 , 0 , , 1 , t. s. ) ( max ) ( 1 1 1 1 j jk n j i u k jk ij n j u k jk jk ,u , ,n, k , j y m i b y a y s y l MKP j j K K K where, ∑ = = j u k j jk x y 1 , 1 , - - = k j jk jk f f s and 2 k d k c f j j jk - = .
In the second step of the algorithm, a surrogate relaxation is applied to the LPrelaxation of (MKP). The resultant formulation (KP, w) is the surrogate relaxation problem of (MKP) and can be written as:

( ) [ ] { } ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ = = ∈ ≤ ∑ ∑ ∑ ∑ ∑ ∑ = = = = = = 1 1 , 1 , 0 t. s. max ) , ( 1 1 1 1 1 1 j jk n j m i i i u k jk m i ij i n j u k jk jk ,u , ,n, k , j y b w y a w y s w KP j j K K
As proved by Glover [START_REF] Glover | Improved linear programming formulations of nonlinear integer problems[END_REF], (KP, w) is a relaxation of (MKP). For any value of

0 ≥ w [ ] w KP Z ,
the optimal value of ( ) proposed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF].

Section 3. Direct expansion of integer variables with additional constraints

In this section we apply to the integer variables of (QMKP) the direct expansion with additional constraints (II). That is each variable x j is replaced by the following expression ∑ = ( ) Following the upper bound method developed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF] we then establish the surrogate problem (KP 2 ,w) associated to (MKP 2 ). The problem (KP 2 ,w) can be written as: Utilizing the optimal solution * ω of the dual of ) ( 2MKP as the surrogate multiplier, the optimal value of the LP-relaxation of (KP 2 , * ω ) provides an upper bound of (QMKP). The following proposition and its corollary show that this upper bound, computed in this section through the use of expansion (II), is equal to the one computed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. . We denote:

( ) { } ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ = = = ∈ ≤ = ≤ = ∑ ∑ ∑ ∑ ∑ = = = = = 1 1 ,..., 1 , 1 , 0 1 , , 1 , t. s.
( ) ( ) { } ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ = = = ∈ ≤ ≤ = ∑ ∑ ∑ ∑ ∑ ∑ ∑ = = = = = = = 1 1 ,..., 1 , 1 , 0 1 , t. s.
-

∑ = = j u k jk j y k x 1 ' . Obviously [ ]; , 0 j j u x ∈ - ⎣ ⎦ j j x p =
, where ⎣ ⎦

x stands for the highest integer lower than or equal to the real x. j p is the integer part of j x . Note that

j j u p < since j j u x ≤ ; - j j j p x - = ε , the fractional part of j x : j j j p x ε + = . Note that [ [. 1 , 0 ∈ j ε
We define the solution y as follows: . This implies that y is feasible for ( ) w KP, since we have:

{ } ⎪ ⎩ ⎪ ⎨ ⎧ ≤ ≤ + = = ≤ ≤ = ∈ ∀ + 2 for 0 1 for 1 , , , 1 1
∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ = = = = = = = = = ≤ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ m i i i m i u k jk n j ij i m i j n j ij i m i u k jk n j ij i ,w KP y b w y k a w x a w y a w j j 1 1 1 ' 1 1 1 1 1
{ } ∑ ∑ = = ≤ ∈ ∀ j j u k jk jk u k jk jk y s y f n j 1 1 ' , , , 1 K Let { } n j , , 1 K ∈
. We first note that:

( ) ( ) ( ) 1 , , , 1 , 1 1 , 1 , 1 1 1 + + = - + = = + - = - + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ - = + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ = ∑ ∑ ∑
K ε ( ) ( ) ⎪ ⎩ ⎪ ⎨ ⎧ ≥ = ≥ + + + 0 , , 1 t. s. minimize α β α β ε α j jk j j j u k f k p D K 1 0 1 t. s. ) ' ( maximize ' 1 ' 1 ' 1 '
The solution y ~ defined by (

{ }, , , 0 j u k K ∈ ∀ , jk f k ≥ + β α ~ Let us denote ( ) jk p j j p j j j f f p k f k p k j j - - + - + = Π +1 , , ) ( 1 ) (
. We thus have to prove that:

{ }, , , 0 j u k K ∈ ∀ , 0 ) ( ≥ Π k j Using the definition of jk f ( 2 k d k c f j j jk - =
), we find with some easy algebra:

) 1 ( ) 1 2 ( ) ( 2 + + + - = Π j j j j j j j p p d k d p k d k Thus, ) (k j Π
is a degree 2 polynomial in k. Its discriminant is equal to Let us now prove that ) ' ( ) ( y h y g = .

{ }, , , 0 j u k K ∈ ∀ 0 ) ( ≥ Π k j , which implies that ( )
-For 0 j j = , we have: We thus have proved in this Section that solving (QMKP) with a linearization technique and a surrogate relaxation may be equivalently done with expansion (I) or with expansion (II) of the integer variables. We proved that expansions (I) and (II) are equivalent in the sense that they provide the same upper bound. However, expansion (II) involves n added constraints without improving the quality of the upper bound in comparison with expansion (I).

( ) 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 , , ' , 1 , ' , , 1 ' , ,

Section 4. Binary expansion of integer variables

This section is dedicated to the use of a binary expansion (referred as expansion (III))

of this integer variables in the upper bound procedure developed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. Such expansion consists in rewriting each integer variable x j (j=1,…,n) as: A possible way to get a further decrease of the number of 0-1 variables would be the use of a classical linearization technique as suggested by [START_REF] Foret | Applications de l'algèbre de Boole en recherché opérationnelle[END_REF] [START_REF] Foret | Applications de l'algèbre de Boole en recherché opérationnelle[END_REF]. This linearization scheme would first consist in applying (III) to the integer variables and then using the basic linearization technique for 0-1 variables. Nevertheless, preliminary computational experiments have shown that the bound provided by this technique is of worst quality than the one computed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. Moreover, this process seems to be more time consuming than [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF].

(

  LAMSADE, Université Paris Dauphine, Place du Maréchal de Lattre de Tassigny, 75775 Paris Cedex 16, France) Eric Soutif (CEDRIC, Conservatoire National des Arts et Métiers, 292 rue Saint Martin, 75141 Paris Cedex 03, France)

1 .

 1 Since the objective function f is separable a linear interpolation can then be applied to each objective function term f j . Consequently, (QMKP) is equivalent to the 0-1 piecewise linear program (MKP):

  w KP, constitutes an upper bound for [ ] MKP Z . But the value of the bound [ ] w KP Z ,depends on the choice of the surrogate multiplier employed. It is proved in[START_REF] Nakagawa | Optimal reliability allocation by branch and bound techniques[END_REF] that if w* is chosen as the optimal solution of the dual of ( ) MKP then the optimal value of (KP,w*) is a tight upper bound for (QMKP) and it is obtained in a very competitive CPU time.Nevertheless, the direct expansion used in the above upper bound procedure results in an increase of the size of the linearized problem. Indeed, the number of 0-1 variables is equal to ∑ = well known that fewer variables are included in the program so less running time is consumed. The purpose of the next sections is to try to reduce the equivalent linearized problem size, using other expansion techniques and to evaluate the impact of such techniques on the computation of the upper bound

  ,…,u j and j=1,…,n. Since the integer variables are now replaced by 0-1 variables, we transform the resultant 0-1 quadratic program into a 0-1 linear problem as follows: we replace each objective function termf j (x j ) = c j k -d j k² .The problem (QMKP) is thus equivalent to the following problem (MKP 2 ):

=

  For a given index j, let us consider the following linear problem ( ) j P and its dual problem ( )

  continuous knapsack problem: for a given index j, all the variables jk y , coefficients of these variables in the objective function g(y) are in decreasing order (because of the concavity of function last equality holds for the same reason as in the previous paragraph (the coefficients of variables jk y in the unique knapsack constraint are identical and the coefficients of these variables in the objective function are in decreasing order). This implies that:

  the non applicability of both (III) and a linear interpolation to transform (QMKP) into an equivalent 0-1 piecewise linear problem through the use of only

As previously mentioned, the aim of this study is to compare other expansion techniques in the upper bound process developed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF]. Since, the variables are now binary the next step of the algorithm of Quadri et al. [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF] Proof 4.1 Assume that it is possible to replace each function f j (x j ) = c j x j -d j x j ², for all j from 1 to n, by a linear function g j (z j ) using exactly ⎡ ⎤ ) ( log 2 j u 0-1 variables. We denote by (H) this assumption.

If (H) is true then it should exist coefficients g jk such that:

The assumption (H) must be satisfied for all problem data. Let us set u j = 3, c j = 10 and d j = 4.

Consequently,

-if x j = 1 (cf. ( 2) and ( 3)) then g j1 = 6 (i), z j1 = 1 and z j2 = 0.

-if x j = 2 (cf. ( 2) and ( 3)) then g j2 = 4 (ii), z j1 = 0 and z j2 = 1.

-if x j = 3 then (cf. ( 3)) then z j1 = z j2 = 1 which implies (cf. ( 2)) that 3c j -9 d j = g j1 + g j2 = -6 (iii).

The equations system (I), (II) and (III) has clearly no solution. Consequently there is a contradiction with (H). Proposition 4.1 shows the non applicability of a binary expansion of the integer variables for (QMKP) so as to transform the initial problem into a 0-1 linear program.

Consequently, the upper bound method proposed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF] can not be applied together with expansion (III).

Section 5. Concluding remarks

In this paper we have theoretically compared the use of three techniques to rewrite integer variables into zero-one variables in an upper bound procedure for (QMKP) developed by Quadri et al. [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF], which provides, to the best of our knowledge a bound closer to the optimum than the existing methods. More specifically, we have compared a direct expansion of the integer variables, originally employed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF] with a direct expansion with additional constraints (II) and with a binary expansion (III).

We have proved that (II) provides the same upper bound as the one computed in [START_REF] Quadri | Upper bounds for large scale integer quadratic multidimensional knapsack problems[END_REF] whereas it involves n additional constraints. We therefore do not expect an improvement of the upper bound computational time. Finally, we provide a proof of


Carole Sylvie Campagne
email: sylvie.campagne@irstea.fr

Philip Roche
email: philip.roche@irstea.fr

Frédéric Gosselin
email: frederic.gosselin@irstea.fr

Léïta Tschanz
email: l.tschanz@imbe.fr

Thierry Tatoni
email: t.tatoni@imbe.fr

Expert-based ecosystem services capacity matrices: Dealing with scoring variability

Keywords: Ecosystem services Capacity matrix Confidence index Uncertainty Expert panel size Participatory assessment ES, ecosystem service, ET, ecosystem types, SD, standard deviation, SE, standard errors, CP, cut point, RNP, Regional Natural Park, RNPBP, "Les Baronnies

Capacity matrices are widely used for assessment of ecosystems services, especially when based on participatory approaches. A capacity matrix is basically a look-up table that links land cover types to ecosystem services potentially provided. The method introduced by Burkhard et al. in 2009 has since been developed and applied in an array of case studies. Here we adress some of the criticisms on the use of capacity matrices such as expert panel size, expert confidence, and scoring variability.

Based on three case-study capacity matrices derived from expert participatory scoring, we used three different approaches to estimate the score means and standard errors: usual statistics, bootstrapping, and Bayesian models. Based on a resampling of the three capacity matrices, we show that central score stabilizes very quickly but that intersample variability shrinks after 10-15 experts while standard error of the scores continues to decrease as sample size increases. Compared to usual statistics, bootstrapping methods only reduce the estimated standard errors for small samples. The use of confidence scores expressed by experts and associated with their scores on ecosystem services does not change the mean scores but slightly increases the standard errors associated with the scores on ecosystem services. Here, computations considering the confidence scores marginally changed the final scores. Nevertheless, many participants felt it important to have a confidence score in the capacity matrix to let them express uncertainties on their own knowledge. This means that confidence scores could be considered as supplementary materials in a participatory approach but should not necessarily be used to compute final scores.

We compared usual statistics, bootstrapping and Bayesian models to estimate central scores and standard errors for a capacity matrix based on a panel of 30 experts, and found that the three methods give very similar results. This was interpreted as a consequence of having a panel size that counted twice the minimal number of experts needed. Bayesian models provided the lowest standard errors, whereas bootrapping with confidence scores provided the largest standard errors.

These conclusions prompt us to advocate when the panel size is small (less than 10 experts), to use bootstrapping to estimate final scores and their variability. If more than 15 experts are involved, the usual statistics are appropriate. Bayesian models are more complex to implement but can also provide more informative outputs to help analyze results.

Introduction

Ecosystem service (ES) is a popular and widely recognised concept [START_REF] Burkhard | Mapping ecosystem service supply, demand and budgets[END_REF][START_REF] Fisher | Defining and classifying ecosystem services for decision making[END_REF], and the term 'ecosystem services' has translated from scientific studies into the mainstream vocabulary of stakeholders and experts [START_REF] Jacobs | The Matrix Reloaded: a review of expert knowledge use for mapping ecosystem services[END_REF]. Increasing demand from policymakers like the European Commission has prompted the development of an array of ES http://dx.doi.org/10.1016/j.ecolind.2017.03.043 1470-160X/© 2017 Elsevier Ltd. All rights reserved.

assessment and mapping methods [START_REF] Martnez-Harms | Methods for mapping ecosystem service supply: a review[END_REF][START_REF] Willemen | Editorial: best practices for mapping ecosystem services[END_REF].

One such method that is gaining ground is the "capacity matrix approach", which was even touted as "the most popular ES assessment technique today" [START_REF] Jacobs | The Matrix Reloaded: a review of expert knowledge use for mapping ecosystem services[END_REF]. The capacity matrix is basically a look-up table that links land cover types to ecosystem services potentially provided [START_REF] Burkhard | Landscapes' capacities to provide ecosystem services-a concept for land-cover based assessments[END_REF].

Since the "matrix" first introduced by [START_REF] Burkhard | Landscapes' capacities to provide ecosystem services-a concept for land-cover based assessments[END_REF], the method has been developed and applied in an array of case studies (e.g. [START_REF] Hermann | Assessment framework for landscape services in European cultural landscapes: an Austrian Hungarian case study[END_REF][START_REF] Kroll | Rural-urban gradient analysis of ecosystem services supply and demand dynamics[END_REF][START_REF] Stoll | Assessment of ecosystem integrity and service gradients across Europe using the LTER Europe network[END_REF][START_REF] Vihervaara | Ecosystem services-A tool for sustainable management of human-environment systems. Case study Finnish Forest Lapland[END_REF]etc.). Based on experts' knowledge, it gives a quick assessment of ES potentially provided in an area [START_REF] Vihervaara | The use of detailed biotope data for linking biodiversity with ecosystem services in Finland[END_REF][START_REF] Stoll | Assessment of ecosystem integrity and service gradients across Europe using the LTER Europe network[END_REF]. The ES concept makes the matrix method mobilize the ES concept in way that makes it easy for stakeholders to understand and appropriate. It is a pedagogical tool that has proven its utility by targeting priorities and highlighting management hotspots. The approach can be applied at different scales (e.g. [START_REF] Stoll | Assessment of ecosystem integrity and service gradients across Europe using the LTER Europe network[END_REF][START_REF] Hermann | Assessment framework for landscape services in European cultural landscapes: an Austrian Hungarian case study[END_REF]. A lack of quantitative data and their spatial heterogeneity raises issues that can be bypassed by asking experts to estimate scores. Depending on the concertation process applied, data can be developed by consensus among the different experts of a territory. The method is also flexible enough to integrate all kinds of data-from models or measurements alike [START_REF] Burkhard | Ecosystem service potentials, flows and demands-concepts for spatial localisation, indication and quantification[END_REF].

Some researchers have started to study the limits of this method, like [START_REF] Jacobs | The Matrix Reloaded: a review of expert knowledge use for mapping ecosystem services[END_REF] who point out its poor methodological transparency, lack of reproducibility and lack of appropriate factoring on uncertainty. [START_REF] Hou | Uncertainties in landscape analysis and ecosystem service assessment[END_REF] also discuss the uncertainties related to the matrix method. The uncertainty of experts judgements is often cited as a limit, but few have analyzed it or integrated it in ES studies [START_REF] Seppelt | A quantitative review of ecosystem service studies: approaches, shortcomings and the road ahead[END_REF][START_REF] Hou | Uncertainties in landscape analysis and ecosystem service assessment[END_REF][START_REF] Vihervaara | The use of detailed biotope data for linking biodiversity with ecosystem services in Finland[END_REF]. Based on the combination of experts' judgments, the scoring in the capacity matrix may carry two sources of uncertainties:

-Variability among experts: variability of the expertise within the chosen experts and of more general knowledge (professional or personal knowledge depending on their experiences) [START_REF] Hou | Uncertainties in landscape analysis and ecosystem service assessment[END_REF]. -Variability of each expert: the confidence the expert has in his/her own scores [START_REF] Jacobs | The Matrix Reloaded: a review of expert knowledge use for mapping ecosystem services[END_REF].

The objective of our study is twofold. First, we aim to integrate the different variabilities in the final score of the capacity matrix and present and compare different approaches to computing them. Second, we aim to identify a minimal size of the expert pool needed to obtain a reliable estimate of the mean of the scores and a small SE of this mean.

In order to meet these objectives, we begin by defining our capacity matrices, the experts, and how the scoring was done. The confidence score we added on the capacity matrix is also detailled. In the second section, we present three sets of approaches to the final scores on a matrix: raw parametric approaches (mean and weighted mean), bootstrap models, and Bayesian models. For each approach, we have two calculations: means of scores that experts expressed, and means integrating a metric of the expert's confidence on his/her own scores. The final scores in the final matrix that incorporates all experts scores is thus estimated with 6 calculations. Most existing capacity matrices (e.g. [START_REF] Stoll | Assessment of ecosystem integrity and service gradients across Europe using the LTER Europe network[END_REF] using expert knowledge are expressed as mean of scores of the expert panel, so we start by presenting the raw parametric approach. The bootstrap model enables to estimate different statistics by assuming an independent sampling from an unknown distribution and to integrate uncertainties. The Bayesian methods that we used are elaborate parametric statistical models that enable to integrate and estimate the different kinds of uncertainties in the statistical anal-ysis. We restrict ourselves to these three sets of statistical methods. This work is the first comparison of three calculations applied on capacity matrix scores. We are not setting out to identify the best calculation of final score but to highlight the various possibilities and their related advantages and disadvantages. In the third section, we present the results of the three calcultations on one matrix and look at the final scores and their variabilities on three capacity matrices with a growing number of experts. Finaly, we conclude with recommandations on using the capacity matrix approach.

The data used in this paper came from three ES assessments: ES provided by land-cover types in the 'Baronnies Provenc ¸ ales' Regional Natural Park (RNP) (associated scores noted RNP-BP) and two ES assessments in the 'Scarpe-Escaut' RNP in northern France-one on ES provided by wetlands (associated scores noted RNP-SEwet) and another on ES provided by all land-cover types (associated scores noted RNP-SEall).

Data

Study sites

The Baronnies Provenc ¸ ales RNP (http://www.baronniesprovencales.fr/) is a sub-mountainous rural area in Southern France located at the crossroad between the Alps and Provence influences. Created in 2015, it is the latest RNP in France, taking the total to 51. The capacity matrix made in 2014 was based on the Park project of 2350 km 2 and 130 municipalities. This nature-preserve territory is recognised nationally for its unique landscapes, rich "terroir", built heritage (terraces in dry stone, hilltop villages) and agriculture (orchards, olive groves, linden, lavender, thyme, rosemary, and more), as well as its remarkable geology and biodiversity.

The Scarpe-Escaut RNP (http://www.pnr-scarpe-escaut.fr/en) in northern France, near the Belgian border, extends over 430 km 2 crossed by the Scarpe and the Escaut rivers with 55 municipalities. It is the oldest of the 51 French RNP. It is also the largest European park, as together with its Belgian neighbour, the Plaines de l'Escaut Natural Park, they form the Hainaut cross-border Nature Park. The Scarpe-Escaut RNP is especially marked by the wet lowland plain around the Scarpe and the Escaut rivers. As a peri-urban area, urban pressure is high (use of space) in a landscape formed by a mosaic of agricultural and natural environments (crops, grasslands, woodlands, marshes, ponds. . .) and urbanized areas. Water is everywhere, and man has been managing it for centuries to develop key activities (drainage, land use, channeling of rivers. . .). For decades, wetlands have been considered a less attractive landscape, and wet meadows have been declining under urban pressure or exploited as profitable sites for agricultural production such as livestock and as landscaped recreational ponds. Perceptions of wetlands today are either bad for certain local stakeholders or nonexistent for the wider community, despite their importance as ES provided to the territory. This perception deficit prompted a study of the ES provided by the wetland types in 2015. After the positive local feeback on the initiative, the method was applied in 2016 to all land cover types.

Data

The capacity matrices

We define ES as goods or services provided by ecosystems that directly or indirectly benefit humans (Millennium Ecosystem Assessment, 2005). For provisioning services and regulating services, the ES list has been based on the European CICES classification [START_REF] Haines-Young | CICES V4.3 -Report Prepared Following Consultation on CICES Version 4[END_REF]. We considered provisioning services, regulating services and cultural services.

Besides the score value for each ES/ecosystem type combination, the experts were asked to provide an index of confidence in their score for each ES and each ecosystem type for the last two datasets. This confidence index is used to estimate the impact of expert confidence on the capacity score and their variations.

As this paper focuses not on capacity scores to estimate ES capacity per se but on the methodology used to obtain expert-based capacity matrices and evaluate score confidence levels, detailed ES and ecosystem types for the datasets are given as supplementary materials. Our dataset is composed of three sets of expert scoring exercices used to produce three different capacity matrices.

The first set was completed in 2014 and consists of a panel of 23 experts in the RNP-BP. This capacity matrix counts 33 ecosystem types (ET) (15 agriculture habitats, 12 forest habitats and 6 aquatic and urban habitats) related to land cover types used in the regional land cover map derived from the National Forest Inventory (IFNv2), the agricultural typology of the Ministry of Agriculture, Food and Forestry, and the Agency for Services and Payment (ASP; 2010), called the Graphic Parcel Register (2011), as well as the CRIGE-PACA regional land cover classification called Ocsol (2006) (a regional implementation of Corine Land Cover with a higher spatial resolution of 1/50,000) by 22 ES (8 provisioning services, 9 regulating services and 5 cultural services), giving 726 scores [START_REF] Tschanz | Services écosystémiques. Pour une gestion durable des territoires, pensons en bouquet[END_REF]. This total dataset is thus 726 ES/ET expert scores × 23 experts.

The second set was completed in 2015 and consisst of a panel of 30 experts in the Scarpe-Escut RNP dealing exclusively with wetland (RNP-SEwet). In this typology, ET is associated with wetland environments based on a local wetlands map. We identified nine ET (2 forest types, 3 agricultural types, 4 semi-natural types) and 17 ES (5 provisioning services, 7 regulating services and 5 cultural services), giving 153 scores. This total dataset is 153 ES/ET expert scores × (9 + 17) confidence scores × 30 experts.

The third set was completed in 2016 and consists of a panel of 17 experts in the Scarpe-Escaut RNP as a whole (noted RNP-SEall) intended to be more exhaustive as it integrates all land cover types. This capacity matrix counts 33 ET (6 aquatic habitats, 13 agriculture habitats, 7 forest habitats and 7 anthropized habitats, corresponding to the land cover types in the ARCH map, www.archnature.fr) by 25 ES (9 provisioning services, 11 regulating services and 5 cultural services), giving 726 final scores. The total dataset is 726 ES/ET expert scores × (33 + 25) confidence scores × 17 experts.

Theses three matrices were produced in different contexts and with differents objectives. While the RNP-SEwet matrix was completed in a context of experts' concertation and outreach education to local citizens, the other two were completed for a more exhaustive spatial assessment, so we did not use identical ET or ES typologies but instead ran local contextual adaptations based on the land cover maps available on each application site.

Expert-based participatory method

The capacity matrix aims to produce estimates of the biophysical capacity of the ecosystem to provide ES. We restricted the panel of elligible experts to those having both local and global ecological knowledge in order to take into account all major ecosystem types and all major activities applied on them. We generally tried to follow the recommendations of [START_REF] Jacobs | The Matrix Reloaded: a review of expert knowledge use for mapping ecosystem services[END_REF] on forming a relevant sample of experts with specific affinity to their territory. Our definition of 'expert' was a person with extensive knowledge or skills based on research, experience, or occupation in a particular field.

The experts considered included researchers with expertise in ecology and/or ES, project or site managers, technicians working on environmental or ecological fields, and heads of territorial organizations such as the water agency, regional chamber of agriculture, regional professional centre for forest owners, regional environ-mental science council, regional conservatories of natural areas, national botanical conservatory, local or regional departments of environmental affairs. The departmental federation of hunters and fishers was also represented, along with associations for environmental protection or naturalists.

For each of the ES assessments, the participants were invited to a workshop dedicated to filling out capacity matrix scores. Any participant who was unable to attend the workshop was interviewed individually. During the workshop and the personal interviews, participants were informed on the state-of-the-art in ES, the study, the methods, and the list of ES and ET. We took time to discuss and clarify all the definitions involved in the matrix and the scoring. We proposed to complete the matrix by columns and to give a score by comparing the different ES capacities of each ET. Our experience has proven that this was an effective way to fill in the matrix. The workshops lasted one day, with the morning session dedicated to presenting all definitions and the afternoon left for participants to fill in the matrices and discuss their understanding of the scorings. In both situations (workshop or personal interviews), we let people give their own score independently. The difference between the two approaches lies in the dialogue initiated in the workshop after the scoring, where everyone was given time to voice their chosen score and open a dialogue on any divergences. At the beginning of the workshop, we defined and specified the rules of dialogue: freedom of speech and to hold divergent opinions, and the possibility of constructive criticism.

Scoring

The capacities of the different ET to provide an individual ES were quantified following the scale given in [START_REF] Burkhard | Landscapes' capacities to provide ecosystem services-a concept for land-cover based assessments[END_REF], i.e. "1 = low relevant capacity, 2 = relevant capacity, 3 = moderate relevant capacity, 4 = high relevant capacity, 5 = very high relevant capacity" and "0 = no relevant capacity". For the RNP-SEwet matrix, "0 = no relevant capacity" was not considered since all the selected ES were provided to some level by all the wetland habitats discussed.

In the rest of the paper, "expert score" refers to the potential ES capacity value that individual experts provided and is noted Y k,i,o , where i is the ES, k is the ET and o is the experts (the total of all experts involved is noted n), and "final score" refers to the average potential capacity value computed from a subset or all expert scores and is noted k,i,o .

The confidence score

In any approach based on experts' knowledge, one of the difficulties is to assess the degree of confidence of each individual expert in their own score. As [START_REF] Jacobs | The Matrix Reloaded: a review of expert knowledge use for mapping ecosystem services[END_REF] stressed "confidence reporting is paramount for communication of results and quality comparison". We propose to quantify a score of the confidence the expert has on their scores, where "confidence refers to the degree of confidence in being correct" [START_REF] Jacobs | The Matrix Reloaded: a review of expert knowledge use for mapping ecosystem services[END_REF]. Each expert was asked to state their confidence in their knowledge on ET and ES using a confidence score ranging from: "1 = I don't feel comfortable on my score, 2 = I feel fairly comfortable on my score and 3 = I feel comfortable on my score".

Thus, for each expert (noted o) of the RNP-SEwet matrix, we had 17 confidence scores for each ES noted V ES i,o and 9 confidence scores for each ET noted V ET k,o . We considered these ET and ES confidence scores as margins of a table, and computed a confidence score for each ET and ES by multiplying them:

V k,i,o = V k,o × V i,o .
The confidence score obtained could take values of 1,2, 3, 4, 6 or 9 that we recoded on a scale of 1-6 (by changing the original "6" to "5" and the original "9" to "6") to faciliate the analysis.

For the bootstrapping and the Bayesian approach incorporating the confidence score, we first estimated the variation of ES × ET score associated with each confidence score value based on an additional direct survey with 10 experts of the RNP-SEwet matrix (allocated within the different ecological specificities/expertises and knowledge) in order to express the uncertainty underlying their confidence scores (cited as levels of uncertainty and noted U k,i,o ). Through concrete examples (combinations of ES and ET with different confidence values), we asked them to express a range of ES × ET scores (0; 0.33; 0.5; 1 or 2 up and down their own score) that they could have considered on a set of specific ET × ES combinations with different associated confidence scores. We tried to cover as many of the experts' questions of possible combinations of confidence scores V k,i,o . From this dataset we computed the SE of ET × ES scores associated with each confidence score.

Statistical analysis

Usual statistics

Most capacity matrices using expert knowledge are expressed as mean of the score values of the expert panel.

First, a final capacity matrix can be obtained with usual direct statistics through the mean, also cited as average score (noted 1 k,i )

and its SD (noted 1 k,i ). From the SD, the SE of the mean score is estimated by dividing by √ n calculated as ŝe 1 k,i and as follows for the SEs of the mean:

ŝe 1 k,i = ˆ 1 k,i √ n
We also used the confidence index score (1-6) as a weighting coefficient for the usual statistics. The weighted mean (noted

1 k,i,o
) can be calculated with different values depending on study, context, etc. Here we used confidence score as weighting factor for the expert scores, so the weighted score is the expert's score multiplied by confidence score, as follows:

Y k,i,o = Y k,i,o × V k,i,o
Accordingly, the weighted mean of the scores is calculated as follows:

1 k,i = i,o Y k,i,o i,o V k,i,o
Then the SD of the weighted means (noted ˆ 3 k,i )is calculated as follows:

ˆ 1 k,i = i,o Y k,i,o -3 i 2 * V k,i,o i,o V k,i,o
From the SDs, the SEs of the mean of scores is estimated by dividing by √ n and thus calculated as ŝe

1 k,i = ˆ 1 k,i
√ n and as follows for the SEs of the weighted mean:

ŝe 1 k,i = ˆ 1 k,i √ n

Bootstrapping statistics

Assuming an independent sampling from an unknown distribution, the bootstrap method enables to estimate different statistics in a better way. The Monte Carlo simulation approach provides approximations of different statistics by independently repeating N times a resampling with replacement from a given set of n observations. Each resampling provides an independent realization from the initial dataset from which N realizations of the statistics are produced, providing a distribution of the statistics considered. The boostrapped Monte Carlo simulations were done using R3.3.1, which features the "boot" package but we preferred, here, to write our own script for greater flexibility.

In our case, n, the number of observations, is equal to the number of experts. We estimated two statistics: 2 k,i the mean of scores and ŝe

2 k,i = ˆ 2 k,i
√ n the SE of the mean of scores, with N the number of repetitions being equal to 500. Number of repetitions is not usually crucial once it passes 100 [START_REF] Efron | Nonparametric estimates of standard error: the jacknife, the bootstrap and other methods[END_REF].

With the bootstrapping statistics, the variation in ES × ET score associated with each confidence value was used in a bootstrap Monte Carlo approach to estimate 2 k,i and its SE (noted ŝe 2 k,i ) taking into account the confidence score associated to Y k,i . We added in the bootstrap script a random error ε k,i,o drawn from a normal distribution of errors associated with each level of the confidence scores, which gives ε k,i,o ∈ N 0, ˆ k,i,o,c with c being the value of the confidence score (0,. . .,6).

Bayesian approach

For the Bayesian analyses, we used two versions of a heteroscedastic ordered probit model [START_REF] Litchfield | An analysis of life satisfaction in Albania: an heteroscedastic ordered probit model approach[END_REF]; see Supplementary Materials) for Y k,i,o . Indeed, the probabilities associated with each observation Y k,i,o were calculated as if the observations stemmed from a latent random variable ω k,i,o following a normal distribution according to the following formulas:

P(Y k,i,o = 1) = P(ω k,i,o < CP 1,i ) P Y k,i,o = 2 = P CP 1,i < ω k,i,o < CP 2,i P Y k,i,o = 3 = P CP 2,i < ω k,i,o < CP 3,i P Y k,i,o = 4 = P CP 3,i < ω k,i,o < CP 4,i P Y k,i,o = 5 = P ω k,i,o > CP 4,i
where the cut-points CP j,i were (i) fixed to 0 and 1 for CP 1,i and CP 4,i ; and (ii) estimated between 0 and 1 within the model for CP 2,i and CP 3,i , via the following formulas:

CP 2,i = 1/ 1 + exp -CP0 2,i and CP 3,i = 1/ 1 + exp -CP0 2,i -exp -CP0 3,i with both
CP0 2,i and CP0 3,i having a Gaussian prior distribution of ℵ (0, 5).

CP j,i cut-points were considered to vary from one ES to another because we conjectured that the span between score values would vary from one ES to another. We checked the variations of the estimators of CP j,i with i and found that, for some ES at least, the cut-point values were very significantly different from other values.

The two parameters describing the distribution of ω k,i,o were the location or mean parameter ˆ 3 k,i,o and the dispersion parameter ŝe where: (i) k,i is a random effect used to parametrize SD level; it followed a Gaussian prior distribution of ℵ ( , ) where had a Gaussian prior centered on 0 with 5 as SD and had a relatively non-informative prior uniform distribution between 0 and 100;

(ii) ␣ k,i is the mean random effect estimated for the latent variable for each ES i and ET k, with a non-informative Gaussian prior distribution of ℵ ( ˛, ˛) where ␣ and ␣ had the same priors as and ; (iii) o is a random observer o effect on the mean, chosen here to be constant across ETs and ESs, and drawn from a common Gaussian distribution centered on 0 and with SD obs :

( o) o ∼ℵ (0, obs )
where obs has the same prior as the k,i s.

(iv) confidence score

• V k,i,o is the product of the confidence scores of each ET and ES (V ET k,o V ES i,o
), restricted to its order (from 1 to 6);

(v) the confidence multipliers ( i ) 1≤i≤6 in model B2 were constructed with no constraint on their order, so that the vector summed to one (to prevent convergence issues with the k,i s). This was done by first setting ' 1 = 1/S, i = i /S for 2 ≤ i ≤ 6 with S = 1 + 6 i=2 i and the i chosen from a uniform prior distribution between 0 and 100. We first tried a version that constrained the order of these parameters to decrease with i, but we met strong convergence issues that forced us to adopt this unconstrained structure.

Furthermore, we also integrated in model B2 the data related to the calibration of confidence scores U k,i,o for specific scores Y k,i,o corresponding to cases (k,i) with controlled levels of confidence V ET k,o and V ES i,o . U k,i,o could take predefined values (0, 0.33, 0.5, 1, and 2) and was assumed to represent the standard variation of Y k,i,o . We introduced the information provided by U k,i,o through a discrete distribution on the 5 possible values proportional to Gaussian distribution ℵ k,i,o + , U , where: (vi) k,i,o is the estimated SD of Y k,i,o as calculated from the above equations for P Y k,i,o = y 1≤y≤5 ;

(vii) is a positive mean SD multiplier chosen to multiply the underlying SD from a relatively non-informative prior uniform distribution between 0.1 and 5;

(viii) is a negative number with a relatively informative uniform prior between -2.0 and 0.0 used to scale k,i,o to make it fit better to the scale of U k,i,o ;

(ix) U is a mean SD with an informative prior distribution between 0 and 1.

We mostly used informative priors here because lessinformative priors gave much noisier estimators, which was contrary to expectations (as we did not expect the mean of U to be above 10.0 given the variation of the Us, whereas we observed such values with non-informative priors).

Model B2 therefore incorporated variations of precision and mean between ES and ET, as well as potential systematic biases between observers. This model also incorporated additional information, i.e. the confidence score provided by experts together with their raw score, as well as a calibration of this confidence in an independent survey. In the rest of the paper, for Model B2, ˆ 3 k,i,o is mean parameter and ŝe 3 k,i,o is the SE. These two models plus three additional models (see Supplementary Materials) were fitted with Stan from R 3.3.0 through the rstan library (2.9.0-3) (Stan Development Team, 2015). We used 3 Markov Chains with a total number of iterations of 100,000 each, a burn-in length of 20,000 and a thinning parameter of 20. Mixing was correct and had 12,000 output values for each parameter, with a number of effective values (CODA R library) above 11,000 for Model 1 and above 3500 for Model 2 (for the worst parameters).

Based on the posterior samples of parameters, we calculated the posterior mean and SEs of the mean across observers of several ES (on the original scale of Y) provided by the different land-cover types studied.

Resampling

In order to estimate the impact of expert panel size on the central and dispersion parameters, we used a bootstrap subsampling method. We simulated different expert panel sizes by subsampling each dataset. We fixed the size of the panel "n" between 2 and the actual number of experts (30 for the RNP-SEwet matrix, 23 for the RNP-BP matrix and 17 for the RNP-SEall matrix), and for each size we generated samples with replacement from the actual data. At each iteration, several statistical parameters were computed. This approach enables us to evaluate the variation of the statistics given smaller-than-actual samples. We used it to try to identify a minimal size of the expert pool needed to obtain a reliable estimate of the final scores and a small SE of this mean.

Results

Central values, capacity matrix scores

Based on the RNP-SEwet dataset that we analyzed in more detail than the other datasets, we observed that the different statistical methods used to compute central parameters (mean) provided very similar results. Table 1, which presents mean and SD of all expert scores (153 scores × 17 expert scores; first column) and of the mean scores of our approaches (means of the 153 final scores with the six different calculations), shows the proximity of their mean results clearly visible at around 3.3 with smiliar values to the third decimal place for all experts scores, Mean of averages scores, Mean of bootstrap scores and Mean Bayesian scores. The mean of all expert scores (calculated with normal theory statistics) is 3.363, and SD at 1.368 (Table 1) with median at 4. For the mean of average scores (mean of the 30 expert scores for each of the 153 ET × ES), the SD is 0.967.

The SD on all means' scores has the same proximity of results as the means except for all experts scores which show higher variation.

Specific ET/ES scores (see Table S.2 in Supplementary Materials) show for example the mean 1 i,o which ranks from 1.40 for ES10 "Contribution to freshwater" in ET8 "Crops" to 4.93 for ES17 "Knowledge, scientific and educational" in ET6 "Marsh, bogs and reed beds". Examination of means between approaches and specific ES and ET finds no distinction between means of the 6 calculations.

Scores variability

We examined means of SEs of the scores (based on the 153 ES × ET combinations) estimated from each method as presented in Fig. 1 and found that variability in mean scores was lowest for weighted scores (0.168) and highest for the Bootstrap scores with confidence score (0.195). The three approaches without confidence score shared similar variability in mean scores, at around 0.173. The bootstrap SE with integration of expert confidence index provided the highest values (14% higher than without confidence index), as expected by including another source of variability (i.e. the individual expert range of values based on confidence on the raw score). This does not mean that this bootstrap is less precise but that, by construct, it integrates more sources of variability than the other method. The Bayesian model with confidence scores also incorporates intra-expert variability and provided higher variability in mean scores than without confidence score (0.1794 vs 0.1746).

As examining the whole capacity matrix in detail is such a cumbersoem task, we propose to focus on the SE of ES9 "Domestic food and their outputs" for each ET. The SE of Bootstrap scores with confidence scores has the highest values of 8-9 of the SE presented in Fig. 1. However, there are no clearly visible patterns to emerge between the other five approaches. The results presented in Fig. 2 for ET8 ("Crops") are more dispersed than for other ET. For this specific ET, the SE of Bayesian models are lower and the SE of weighted scores are higher than with other approaches. When we look closer at the ET8 × ES9 scores (available in Supplementary Materials Table S.2), the mean of average score is 4.73 and the mean of confidence index is 4.37.

Effect of expert confidence

During the different ES assessments, several experts accepted to fill the matrix only when offered a confidence score giving the opportunity to express theirs self-confidence on their scores. In the RNP-SEwet matrix, mean confidence score per expert shows a trend towards high confidence scores for 10 experts out of 30 (median at 3) while two experts tended to have low confidence score with a median at 1. The table S.3 (Supplementary Materials) presents V k,i,o in the table and the means of each V ES i,o and V ET k,o in margin. The mean confidence score of all experts for all ET and ES (noted

V ET k,o andV ES i,o
, respectively, and rating between 1 and 3) is 2.23 with each confidence score (1, 2 and 3) represented at 20%, 37% and 43%, respectively. Mean score of all ES was 2.24 and mean score of all ET was 2.21. A more detailed analysis shows there is no uniform profile for ET or ES except for cultural services where each confidence score (1, 2 and 3) represented at 12%, 28% and 60%, respectively.

The results from additional direct survey presented in Fig. 3 show the range and means of variation in expert scores (the levels of uncertainty U k,i,o ) related to each confidence score (with V k,i,o = V k,o × V i,o recoded on the 1-to-6 scale used in the different models).

Note the gradual decrease in score variation as confidence score increases. The variability of the ES score errors associated with confidence scores was higher for lower confidence values (1-2-3) than for higher confidence values (4-5-6). The high values of V k,i,o are related to high confidence of the experts on theirs scores and thus to lower variation in their final scores. There was no significant difference between scores errors associated with confidence score 5 and 6.

Effect of expert panel size

Expert panel size had an effect on two aspects of variability. First, it had an effect on the stability of the final score obtained from different sets of resampled same-size tables: SD of resampled final scores drops sharply from 2 to 10 (Fig. 4) but then remains relatively stable before flattening after a size of 15 for the three matrices. This means that above an expert size panel of 15 experts, the variation of final score obtained from different panel compositions is stable. The total reduction of SD is approximately 0.2 (Fig. 4), i.e. about 16% lower when we have 15 experts involved compared to 2.

Second, we tested the effect of expert panel size on the SEs of final mean score (the precision of estimates). As expected, the SEs decrease with the sample size. The bootstrap estimates show lower values than usual statistics when the sample size is below 10-15 experts (Fig. 5). All three datasets show similar trends. However, even with an expert panel size at 30, we still do not reach stabilization.

Discussion

Final scores and variability

The mean is the most widely used statistic for computing a final score in the capacity matrix approach, but some authors advocate using median scores instead [START_REF] Kopperoinen | Using expert knowledge in combining green infrastructure and ecosystem services in land use planning: an insight into a new place-based methodology[END_REF]. Different kind of weights can also be applied to obtain the final score. Depending on the approach, weighting can be applied during the extrapolation at larger geographic scale, as in [START_REF] Hermann | Assessment framework for landscape services in European cultural landscapes: an Austrian Hungarian case study[END_REF] who chose to weight scores by the area of base unit. [START_REF] Koschke | A multi-criteria approach for an integrated land-cover-based assessment of ecosystem services provision to support landscape planning[END_REF] elected to apply "explicit weights as the importance of the various ecosystem services might differ with respect to the context, the included stakeholders, and the investigated region". These weights are linked to specific methodological choices. Here, we did not apply ES or ET dependent weightings as we wanted to equally consider every ES and ET. We used a weight based on the confidence scores, but the upshot is that a participant with more confidence than another gives a better or a more realistic estimate of the habitat's capacity to provide the ES.

With 30 experts involved in the RNP-SEwet dataset, the different statistical methods used to compute parameters provided very similar results. The results show only marginal differences in final scores when we integrate confidence score or the uncertainty underlying confidence scores in general observation of the data (mean of means) or at an ES level. Indeed, with 30 experts involved, the final score did not change much after taking the confidence score into account. The variabilities of final scores are also quite similar, except for bootstrap scores with confidence scores which had higher variability. The increase of variability is coherent with the addition of inner confidence or doubt of each expert in terms of variability of scores.

Regarding the Bayesian approach, the variability is similar for the two models considered, which may be due to the fact that both integrate an observer effect. The observer effect may integrate some systematic bias of each expert when providing their scores. Nevertheless, the final results of the two Bayesian models are close to those of the other methods. However, the Bayesian framework has the potential to deliver substantially more information than just central values and errors-it can clearly integrate all the information and provide estimates on these other pieces of information. To illustrate, the Bayesian method can reach an estimate of the variation of precision of the latent process as a function of confidence score. If we had observed strong systematic biases between observers-which was not the case here, as the mean estimate of obs was as low as 0.10-Bayesian methods would have made it possible to provide estimates of SEs that would have removed these biases, which other methods are unable to do.

Given the amount of information contained in our capacity matrices, it was not possible to separately detail all the ES and ET combinations, so we have detailed one example combination here, and provided all the tables as Supplementary Materials for interested readers. For the specific ES × ET combinations looked at here, the results of ET8 × ES1 scores can be explained by the nature of this score. "Crops" have a far higher potential capacity to provide "Domestic food and their outputs" than other ES types. The variability of its score is low due to the logical general agreement on this potential capacity. The Bayesian models reflect this general agreement "better" as the lowest SE. Only the weighted score has a surprising result compared to other methods, whichmay be explained by the fact that the confidence score is less a reflection of the variability of each participant than a reflection of the uncer-tainty underlying confidence scores expressed through the results of the additional survey. et al. (2014) recommended following the IPCC (2014) or Millennium Ecosystem Assessment (2005) methods with participant confidence reported with five levels of confidence or assessing the reliability of outcomes in a similar standardized manner. In our case, asking for such specific confidence scores was too much of a burden given the time and thought already invested to fill out the capacity matrix. Additional information like a confidence score has to be suitable, understandable, and on a different scale to the participants' score (5 levels).

Importance of experts' confidence in their scores

Jacobs

Many participants expressed the importance of having a confidence score in the capacity matrix to let them moderate their self-reliability on their own knowledge. Analysis of the confidence score showed that 80% of the participants felt "comfortable or moderately comfortable on their scores". It is important to add a confidence score in a participatory process, not just as supplementary materials but also to be able to recalibrate it for integration into the analysis. We consequently chose to use three levels of confidence for ET and three levels of confidence for ES, which when combined provide a 6-level confidence score for each ET and ES. This approach proved to be flexible, workable, and easy for all participants to understand.

The additional survey supports the hypothesis of a correlation between confidence score and SD of the final score, and shows the negative relationship between the confidence score and the mean of score errors. Our confidence index can this be used as a proxy of the errors in each expert's estimate scores. A limitation here is that we have not explored all the possibilities of the Bayesian models. We were able to analyze the noise around expert scores with other approaches but we chose not to report the results in this paper.

How large should an expert panel be?

Without distinction between the three matrices, the variability of the final scores is stable after 15 experts involved and the SE in the final scores decreased with increasing expert panel size. The bootstrap approach also increases the precision by reducing the variation when panel sizes are small.

These results serve to explain our low difference between means of the three approaches and with or without taking into account the confidence score. With 30 experts involved, we reach a stable mean that may thus iron out disparities.

This result is interesting in practice as it is the first time that a recommendation regarding panel size is provided for participatory ES assessment. Having a large panel is good, but smaller panels can still provide reliable estimates. It should be stressed that the case configuration here is panel composition based on experts and scientists with expertise in ecosystems and/or ES in the studied region. Other studies are needed to address the question of panels or surveys based on non-specialists.

The stability of the scores expresses a consensus of the panel on the potential capacity of a given ET to provide a given set of ES. There is still a need to further explore the link between this consensus score and actual potential capacity.

Conclusion

The two parts of this paper bring complementary analysis to help identify what methods to use or how many experts are needed. These questions also depend on the level of error accepted and of course the number of experts that can viably be involved in the participatory approach.

For a restricted number of participants, it is better to use bootstrap estimate parameters in order to increase the stability of the final score and its variability. If more than 15 experts are involved, the usual statistics is sufficient. Indeed, with 30 experts involved in our matrix, final scores showed little distinction between the three approaches, as we had reach a stable mean.

Bayesian methods are more demanding but they also give more complete results. Bayesian methods effectively require a degree of expertise in how to write, fit and interpret Bayesian models. However, they do enable to integrate different kinds of information (here, the expert scores, confidence scores and theirs uncertainties were highlighted) in an integrated statistical modelling framework [START_REF] Gimenez | Statistical ecology comes of age[END_REF] as well as to estimate all the parameters associated with these variables. In particular, Bayesian models hold great potential if the aim is to gain a more inferential point of view. For the sake of consistence and concision, this paper did not explore all the possibilities of the Bayesian models.

For participatory processes, we recommend adding a confidence score, on the caveat that incorporation of a confidence score or the uncertainty it carries to the final score in the capacity matrix depends on the SE accepted and the method used. We found here that considering confidence score in final score induced only marginal differences.

  3k,i,o corresponding to the SE of the latent normal distribution. After comparing four models (see Supplementary Materials), we retained two different models that contained different parametrizations of the heteroscedasticity (i.e. in parameter ŝe 3 k,i,o ). In the first model (model B1), no other information than basic score Y k,i,o was provided to the model, whereas in the second model (model B2), both the confidence scores V k,i,o and the levels of uncertainty U k,i,o were taken into account. The formula for the mean parameter k,i,o of the latent variable ω k,i,o was the same for all models

Fig. 1 .

 1 Fig. 1. Means of standard errors (SE) of the different methods.

Fig. 2 .

 2 Fig. 2. Standard errors (SE) of ES9 "Domestic food and their outputs" for different land-cover types ET1: Wet forest, ET2: Broadleaf monospecific afforestation in wetlands, ET3: Wet meadow, ET4: Ponds and pools, ET5: Mining subsidence ponds and quarry waters, ET6: Marsh, bogs and reedbed, ET7: Wetland fallows, ET8: Crops, ET9: Green space.

Fig. 3 .

 3 Fig. 3. Mean and standard errors of the variation in expert scores (U k,i,o ) related to each confidence score (V k,i,o ).

Fig. 4 .

 4 Fig. 4. Effect of expert panel size (number of participants) on the variation (SD) of means of average scores.

Fig. 5 .

 5 Fig. 5. Effect of expert panel size (number of participants) on the mean of SEs of ES × ET scores (average of the whole capacity matrix) using usual statistics (SE) and bootstrapped statistics (SE Boot).

Table 1

 1 Statistics (mean and Standard Deviation SD) using 3 different methods: Usual, Bootstrap and Bayesian. Usual refers to the use of standard parametric statistics. Bootstrap is based on resampling and Monte-Carlo methods. Bayesian is based on heteroscedastic ordered probit models.

			Usual		Bootstrap		Bayesian
		On all expert	Mean of average	Mean of weighted	Mean of	Mean of Bootstrap	Mean Bayesian	Mean Bayesian scores
		scores Y k,i,o	scores 1 k,i	scores 1 k,i	Bootstrap scores 2 k,i	scores with confidence scores 2 k,i	scores ˆ 3 k,i	with confidence scores ˆ 3 k,i
	Mean	3.363	3.363	3.371	3.363	3.298	3.363	3.365
	SD	1.368	0.967	1.006	0.967	0.903	0.932	0.935

Acknowledgements

We thank Jean-Michel Salles for his input and comments on the paper. This study was supported by the Natural Regional Park of Scarpe-Escaut. (PNRSE), and we thank Gerald Duhayon and Tangui Lefort from the PNRSE for their support. The RNP-BP matrix is part of a 2017 PhD project by L. Tschanz titled "Towards a socio-ecological observation for a global approach to biodiversity challenges of in regional planning" at the IMBE, Aix-Marseille University. This work was funded by Irstea, the PNRSE, the PNRBP, and the Institut INSPIRE.

Appendix A. Supplementary data

Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.ecolind.2017. 03.043.