Gérard Biau
email: gerard.biau@upmc.fr

Aurélie Fischer
email: aurelie.fischer@upmc.fr

Parameter Selection for Principal Curves

Keywords: Principal curves, parameter selection, model selection, oracle inequality, penalty calibration, slope heuristics. 2010 Mathematics Subject Classification: 62G08, 62G05

Principal curves are nonlinear generalizations of the notion of first principal component. Roughly, a principal curve is a parameterized curve in R d which passes through the "middle" of a data cloud drawn from some unknown probability distribution. Depending on the definition, a principal curve relies on some unknown parameters (number of segments, length, turn. . . ) which have to be properly chosen to recover the shape of the data without interpolating. In the present paper, we consider the principal curve problem from an empirical risk minimization perspective and address the parameter selection issue using the point of view of model selection via penalization. We offer oracle inequalities and implement the proposed approaches to recover the hidden structures in both simulated and real-life data.

Introduction

Principal curves

Statisticians use various methods in order to sum up information and represent the data by simpler quantities. Among these methods, Principal Component Analysis (PCA) aims at determining the maximal variance axes of a data cloud, as a means to represent the observations in a compact manner revealing as well as possible their variability (see, e.g., Mardia, Kent and Bibby [START_REF] Mardia | Multivariate Analysis[END_REF]). This technique, initiated at the beginning of the last century by Pearson [START_REF] Pearson | On lines and planes of closest fit to systems of points in space[END_REF] and Spearman [START_REF] Spearman | General intelligence, objectively determined and measured[END_REF], and further developed by Hotelling [START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF], is certainly one of the most famous and most widely used procedure of multivariate analysis. Whether in the context of dimension reduction or feature extraction, PCA often provides a first important insight in the data structure.

However, in a number of situations, it may be of interest to summarize information in a nonlinear manner instead of representing the data by straight lines. This approach leads to the notion of principal curve, which can be thought of as a nonlinear generalization of the first principal component. Roughly, the purpose is to search for a curve passing through the middle of the observations, as illustrated in Figure 1. Principal curves have a broad range of applications in many different areas, such as physics (Hastie and Stuetzle [START_REF] Hastie | Principal curves[END_REF], Friedsam and Oren [START_REF] Friedsam | The application of the principal curve analysis technique to smooth beamlines[END_REF]), character and speech recognition (Kégl and Krzyżak [START_REF] Kégl | Piecewise linear skeletonization using principal curves[END_REF], Reinhard and Niranjan [START_REF] Reinhard | Parametric subspace modeling of speech transitions[END_REF]), mapping and geology (Brunsdon [START_REF] Brunsdon | Path estimation from GPS tracks[END_REF], Stanford and Raftery [START_REF] Stanford | Finding curvilinear features in spatial point patterns: principal curve clustering with noise[END_REF], Banfield and Raftery [START_REF] Banfield | Ice floe identification in satellite images using mathematical morphology and clustering about principal curves[END_REF], Einbeck, Tutz and Evers [START_REF] Einbeck | Exploring multivariate data structures with local principal curves[END_REF][START_REF] Einbeck | Local principal curves[END_REF]), natural sciences (De'ath [START_REF] De'ath | Principal curves: a new technique for indirect and direct gradient analysis[END_REF], Corkeron, Anthony and Martin [START_REF] Corkeron | Ranging and diving behaviour of two 'offshore' bottlenose dolphins, Tursiops sp., off eastern Australia[END_REF], Einbeck, Tutz and Evers [START_REF] Einbeck | Exploring multivariate data structures with local principal curves[END_REF]) and medicine (Wong and Chung [START_REF] Wong | Principal curves to extract vessels in 3D angiograms[END_REF], Caffo, Crainiceanu, Deng and Hendrix [START_REF] Caffo | A case study in pharmacologic colon imaging using principal curves in single photon emission computed tomography[END_REF]).

The definition of a principal curve typically depends of the principal component property one wants to generalize. Most of the time, this definition is first stated for an R d -valued random variable X = (X 1 , . . . , X d ) with known distribution, and then adapted to the practical situation where one observes independent draws X 1 , . . . , X n distributed as X.

x 1 f (t 1 ) x 2 f (t 2 ) x 3 f (t 3 ) x 4 f (t 4 ) x 5 f (t 5 ) x 6 f (t 6 ) x 7 f (t 7 )
Figure 2: The projection index t f . For all i, t i stands for t f (x i ).

The original definition of a principal curve goes back to Hastie and Stuetzle [START_REF] Hastie | Principal curves[END_REF] and relies on the self-consistency property of principal components. In words, a smooth (infinitely differentiable) parameterized curve f (t) = (f 1 (t), . . . , f d (t)) is a principal curve for X if f does not intersect itself, if it has finite length inside any bounded subset of R d , and if it is self-consistent. This last requirement means that

f (t) = E[X|t f (X) = t], (1)
where the so-called projection index t f (x) is the largest real number t minimizing the squared Euclidean distance between x and f (t), as depicted in Figure 2. More formally, t f (x) = sup t : xf (t) = inf t ′

xf (t ′ ) .

The self-consistency property may be interpreted by saying that each point of the curve f is the mean of the observations projecting on f around this point. Hastie and Stuetzle discuss in [START_REF] Hastie | Principal curves[END_REF] an iterative algorithm, alternating between a projection and a conditional expectation step, which yields an approximate principal curve. As this approach exhibits different types of bias, Banfield and Raftery [START_REF] Banfield | Ice floe identification in satellite images using mathematical morphology and clustering about principal curves[END_REF] and Chang and Ghosh [START_REF] Chang | Principal curves for nonlinear feature extraction and classification[END_REF] propose a modification of the algorithm, whereas Tibshirani, tackling the model bias problem, adopts in [START_REF] Tibshirani | Principal curves revisited[END_REF] a semiparametric strategy and defines principal curves in terms of a mixture model. For more references on principal curves and related points of view, we refer the reader to Verbeek, Vlassis and Kröse [START_REF] Verbeek | A soft k-segments algorithm for principal curves[END_REF] (k-segments algorithm), Delicado [START_REF] Delicado | Another look at principal curves and surfaces[END_REF] (principal curves of oriented points), Einbeck, Tutz and Evers [START_REF] Einbeck | Local principal curves[END_REF] (local principal curves) and Genovese, Perone-Pacifico, Verdinelli and Wasserman [START_REF] Genovese | The geometry of nonparametric filament estimation[END_REF], who recently discussed a closely related approach, called nonparametric filament estimation.

In the present paper, we will adopt the principal curve definition of Kégl, Krzyżak, Linder and Zeger [START_REF] Kégl | Learning and design of principal curves[END_REF], which is slightly different from the original one. The main advantage of this definition, which is recalled in the next paragraph, is that it avoids the implicit conditional expectation requirement (1) and, consequently, turns out to be more easily amenable to mathematical analysis.

Constrained principal curves

In the definition of Kégl, Krzyżak, Linder and Zeger [START_REF] Kégl | Learning and design of principal curves[END_REF] ( ÐÞ hereafter), a principal curve of length (at most) L for X is a parameterized curve minimizing the least-square criterion

∆(f ) = E inf t X -f (t) 2
over a collection F L of curves of length not larger than some prespecified positive L. We note that, in this context, a principal curve always exists provided E X 2 < ∞, but that it may not necessarily be unique. In practice, as the distribution of X is unknown, ∆(f ) is replaced by its empirical counterpart

∆ n (f ) = 1 n n i=1 inf t X i -f (t) 2
based on a sample X 1 , . . . , X n of independent random variables distributed as X. Considering the minimum fk,n of ∆ n (f ) over the subclass F k,L ⊂ F L of all polygonal lines f k,n with k segments and length not larger than L, Kégl, Krzyżak, Linder and Zeger [START_REF] Kégl | Learning and design of principal curves[END_REF] prove that, whenever X is almost surely bounded, and for the choice k ∝ n 1/3 , ∆( fk,n )min

f ∈F L ∆(f ) = O(n -1/3 ).
As the task of finding a polygonal line with k segments and length at most

L minimizing ∆ n (f ) is computationally difficult,
ÐÞ propose an approximate iterative algorithm that they call the Polygonal Line Algorithm. This algorithm is initialized using the smallest segment included in the first principal component containing all projected data points. Then, at each step, a vertex-and thus, a segment-is added to the current polygonal line, and the vertices are updated in a cyclic manner during an inner loop alternating between a projection and an optimization step. Performing the projection step is similar to constructing a Voronoi partition, with respect to both the vertices and segments. To optimize a vertex, a local version of ∆ n (f ) is used, involving only the data projecting to this vertex and to the adjacent segments. The criterion is penalized to avoid sharp angles, which in turn amounts to penalizing the length of the curve.

Working out the angle penalty in the Polygonal Line Algorithm, Sandilya and Kulkarni (× hereafter) propose in [START_REF] Sandilya | Principal curves with bounded turn[END_REF] a closely related definition, by imposing a constraint on the turn (Alexandrov and Reshetnyak [START_REF] Alexandrov | General Theory of Irregular Curves[END_REF]) of the curve f . This approach consists in replacing the class F L by F K , where K stands for the maximal turn. Thus, denoting by F k,K ⊂ F K the subclass of all polygonal lines f k,n with k segments and turn not larger than K, × prove that, whenever X is almost surely bounded, and for the choice

k ∝ n 1/3 , ∆( fk,n ) -min f ∈F K ∆(f ) = O(n -1/3 ).
Whether in the ÐÞ definition or in the × one, selecting the various smoothness parameters (the number k of segments, the curve length ℓ, the turn κ) is an essential issue, as illustrated in Figure 3. A good choice of these parameters is critical, since a principal curve obtained with a poor class will be too rough, whereas a class containing too many curves may lead to severe interpolation problems.

In practice, the Polygonal Line Algorithm stops when k is larger than a certain threshold, chosen heuristically and tuned after carrying out several experiments. The stopping condition involves the number n of observations and the actual value of the criterion ∆ n . However, to our knowledge, this empirical procedure is not supported by any theoretical argument and leads to variable results, depending on the data set. Besides, note that assessing the complexity of Hastie and Stuetzle [START_REF] Hastie | Principal curves[END_REF] principal curve estimates by cross-validation has often been observed to fail. As put forward by Duchamp and Stuetzle [START_REF] Duchamp | Extremal properties of principal curves in the plane[END_REF], these principal curves are saddle points of the distance between a random vector and a curve, and therefore cross-validation is not a well-suited technique in the principal curve framework. As far as we know, the issue of an automatic (i.e., data-dependent) choice of the parameters k, ℓ and κ has not been addressed in the literature. Thus, to fill the gap, we propose in the present contribution to focus on this question both from a theoretical and practical point of view. Our approach will strongly rely on the model selection theory by penalization introduced by Birgé and Massart [START_REF] Birgé | From model selection to adaptive estimation[END_REF] and Barron, Birgé and Massart [START_REF] Barron | Risk bounds for model selection via penalization[END_REF], as well as on a recent penalty calibration approach proposed by Birgé and Massart [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF] and Arlot and Massart [START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF].

[A] [B] [C]
The paper is organized as follows. First, we consider in Section 2 principal curves with bounded length and show that the polygonal line obtained by minimizing some appropriate penalized criterion satisfies an oracle-type inequality. Section 3 provides a similar result in the context of principal curves with bounded turn. Our theoretical findings are illustrated on both simulated and real data sets in Section 4. For the sake of clarity, proofs are collected in Section 5.

Principal curves with bounded length

Let • be the standard Euclidean norm over R d . A parameterized curve in R d is a continuous function

f : I → R d t → (f 1 , . . . , f d ),
where I = [a, b] is a closed interval of the real line. The length of f is defined by

L (f ) = sup m j=1 f (t j ) -f (t j-1 ) ,
where the supremum is taken over all subdivisions a = t 0 < t 1 < • • • < t m = b, m ≥ 1 (see, e.g., Kolmogorov and Fomin [START_REF] Kolmogorov | Introductory Real Analysis[END_REF]). Throughout the document, it is assumed that E X 2 < ∞ and that

P {X ∈ C} = 1, (2)
where C is a convex compact subset of R d , with diameter δ. By Lemma 1 in Kégl [START_REF] Kégl | Principal Curves: Learning, Design, and Applications[END_REF], the requirement (2) implies that, for any given positive length L, there exists a principal curve for X with length at most L in C, that is a (non necessarily unique) parameterized curve f ⋆ with length not larger than L and support in C achieving the minimum of E [inf t∈I Xf (t) 2 ]. Consequently, in the sequel, we will restrict ourselves to curves whose support is included in C and denote by F the set of all parameterized curves f = (f 1 , . . . , f d ) belonging to C.

Let X 1 , . . . , X n be a sample of independent random variables distributed as X, and consider the contrast

∆(f , x) = inf t∈I x -f (t) 2 , f ∈ F , x ∈ R d .
The associated empirical risk based on the sample X 1 , . . . X n is defined as

∆ n (f ) = 1 n n i=1 ∆(f , X i ) = 1 n n i=1 inf t∈I X i -f (t) 2 .
For some prespecified length L > 0, we set

f ⋆ ∈ arg min f ∈F ,L (f )≤L E[∆(f , X)].
Next, let L be a countable subset of ]0, L] and Q a grid over C, that is

Q = C ∩ Γ,
where Γ is a lattice of R d . For every k ≥ 1 and ℓ ∈ L, the model F k,ℓ is defined as the collection of all polygonal lines with k segments, with length at most ℓ, and with vertices belonging to Q. We note that each model F k,ℓ as well as the family of models {F k,ℓ } k≥1,ℓ∈L are countable. For k ≥ 1 and ℓ ∈ L, let fk,ℓ ∈ arg min

f ∈F k,ℓ ∆ n (f )
be a curve achieving the minimum of the empirical criterion ∆ n (f ) over the polygonal line class F k,ℓ .

At this stage of the procedure, we have at hand a family of estimates { fk,ℓ } k≥1,ℓ∈L and our goal is to select the best principal curve f among this collection. To this aim, we make use of the model selection approach of Barron, Birgé and Massart [START_REF] Barron | Risk bounds for model selection via penalization[END_REF], which allows to assess the adjustment quality by controlling the loss

D(f ⋆ , f) = E[∆( f , X) -∆(f ⋆ , X)]
between the target f ⋆ and the selected curve f . (For a comprehensive introduction to the area of model selection, the reader is referred to the monograph of Massart [START_REF] Massart | Concentration Inequalities and Model Selection[END_REF].) More formally, let pen : N * × L → R + be some penalty function and denote by ( k, l) a pair of minimizers of the criterion crit(k, ℓ) = ∆ n ( fk,ℓ ) + pen(k, ℓ).

In order to obtain the desired principal curve f = fk , l, we have to design an adequate penalty pen(k, ℓ). This is done in the following theorem, which is an adaptation of a general model selection result of Massart [START_REF] Massart | Concentration Inequalities and Model Selection[END_REF]Theorem 8.1]. However, for the sake of completeness, it is proved in its full length in Section 5.

Theorem 2.1. Consider a family of nonnegative weights {x

k,ℓ } k≥1,ℓ∈L such that k≥1,ℓ∈L e -x k,ℓ = Σ < ∞,
and a penalty function pen :

N * × L → R + . Let f = fk , l. If for all (k, ℓ) ∈ N * × L, pen(k, ℓ) ≥ E sup f ∈F k,ℓ E[∆(f , X)] -∆ n (f ) + δ 2 x k,ℓ 2n , then E[D(f ⋆ , f)] ≤ inf k≥1,ℓ∈L D(f ⋆ , F k,ℓ ) + pen(k, ℓ) + δ 2 Σ 2 3/2 π n ,
where

D(f ⋆ , F k,ℓ ) = inf f ∈F k,ℓ D(f ⋆ , f ).
Theorem 2.1 offers a nonasymptotic bound, expressing the fact that the expected loss of the final estimate f is close to the minimal loss over all k ≥ 1 and ℓ ∈ L, up to a term tending to 0. Thus, in order to apply this theorem to the principal curve problem, we now have to find an upper bound on the quantity

E sup f ∈F k,ℓ E[∆(f , X)] -∆ n (f ) . ( 3
)
This is achieved by Proposition 2.1 below, which is proved by showing that the expected maximal deviation (3) may be bounded by a Rademacher average (see Bartlett, Boucheron and Lugosi [START_REF] Bartlett | Model selection and error estimation[END_REF] and Koltchinskii [START_REF] Koltchinskii | Rademacher penalties and structural risk minimization[END_REF]) and by resorting to a Dudley integral (Dudley [START_REF] Dudley | The sizes of compact subsets of Hilbert space and continuity of Gaussian processes[END_REF]).

Proposition 2.1. Let F k,ℓ be the set of all polygonal lines with k segments, length at most ℓ, and vertices in a grid Q ⊂ C. Then there exist nonnegative constants a 0 , . . . , a 2 , depending on the maximal length L, the dimension d and the diameter δ of the convex set C, such that

E sup f ∈F k,ℓ E[∆(f , X)] -∆ n (f ) ≤ 1 √ n a 1 √ k + a 2 ℓ + a 0 .
Finally, by combining Theorem 2.1 and Proposition 2.1, we are in a position to state the main result of this section.

Theorem 2.2. Consider a family of nonnegative weights {x

k,ℓ } k≥1,ℓ∈L such that k≥1,ℓ∈L e -x k,ℓ = Σ < ∞,
and a penalty function pen :

N * × L → R + . Let f = fk , l.
There exist nonnegative constants c 0 , . . . , c 2 , depending on the maximal length L, the dimension d and the diameter δ of the convex set C, such that, if for all (k, ℓ)

∈ N * × L, pen(k, ℓ) = 1 √ n c 1 √ k + c 2 ℓ + c 0 + δ 2 x k,ℓ 2n , then E[D(f ⋆ , f)] ≤ inf k≥1,ℓ∈L D(f ⋆ , F k,ℓ ) + pen(k, ℓ) + δ 2 Σ 2 3/2 π n ,
where

D(f ⋆ , F k,ℓ ) = inf f ∈F k,ℓ D(f ⋆ , f ).
Some comments are in order.

Firstly, we see that the penalty shape involves a term proportional to k/n and a term proportional to ℓ/ √ n. This penalty form, which vanishes at the rate 1/ √ n, seems relevant insofar as the number k of segments and the length ℓ of the curves measure the complexity of the models.

Observe next that the proof of Proposition 2.1 provides possible values for the constants c 0 , . . . , c 2 . However, these values are not very helpful since they are upper bounds which are probably far from being tight. Nevertheless, the proof also reveals that

c 1 = c ′ 1 δ 2 , c 2 = c ′ 2 δ and c 0 = c ′ 0 δ 2
, where c ′ 0 , c ′ 1 and c ′ 2 are constants without dimension, so that the penalty is in fact homogeneous to a squared length, just like the criterion ∆ n (f ) is.

Finally, an important practical issue is how to choose the weights {x k,ℓ } k≥1,ℓ∈L . These weights should be large enough to ensure the finiteness of Σ, but not too large at the risk of overpenalizing. If the cardinality of the collection of models is not larger than n 2 (this will be the case in all our practical examples), we may set x k,ℓ = 2 ln n for every (k, ℓ). This choice does not affect the penalty shape, though modifying the rate, and leads to Σ = 1 in the risk bound.

Remark 2.1. Clearly, when the length ℓ of polygonal lines is fixed, and the aim is to select the number k of segments, the dominant term reflecting the complexity of the models in the penalty is k/n. In this particular context, the weights may be taken equal to ln n, or, by analogy with the Gaussian linear model selection framework, proportional to k. Indeed, in this framework, each model S m , m ∈ M, has dimension D m and an interesting choice for x m is then x m = x(D m ), where x(D) = cD + ln |{m ∈ M : D m = D}| and c > 0. When there is no redundancy in the models dimension, this strategy amounts to choosing x m proportional to D m . In our problem, this means setting x k = ck for every k, where the constant c > 0 ensures that k≥1 e -ck = Σ < ∞. Thus, in this somewhat restrictive situation, the penalty is of the order k/n.

Principal curves with bounded turn

As it was already mentioned in the Introduction, Sandilya and Kulkarni [START_REF] Sandilya | Principal curves with bounded turn[END_REF] (× ) suggest an alternative approach for principal curves, based on the control of the turn. Recall that the turn K (f ) of a curve f :

I → R d , I = [a, b], is given by K (f ) = sup m-1 j=1 f (t j ),
where f (t j ) denotes the angle between the vectors --------→ f (t j-1 )f (t j ) and --------→ f (t j )f (t j+1 ), and the supremum is taken over all subdivisions a = t 0 < t 1 < • • • < t m = b, m ≥ 1 (Alexandrov and Reshetnyak [START_REF] Alexandrov | General Theory of Irregular Curves[END_REF]). Thus, the turn of a polygonal line with vertices v 1 , . . . , v k+1 is just the sum of the angles at v 2 , . . . , v k (see Figure 4 for an illustration). As a logical continuation to Section 2, we propose in the present section to analyse the × definition from a model selection point of view. To this aim, we use the fact that a curve with bounded turn also has bounded length, as shown in Lemma 3.1 below.

φ 2 φ 3 v 1 -→ s 1 v 2 -→ s 2 v 3 -→ s 3 v 4
We still assume that P {X ∈ C} = 1, where C is a convex compact subset of R d with diameter δ. By Proposition 1 in × , this requirement ensures the existence of a curve f ⋆ with bounded turn minimizing the criterion E[∆(f , X)]. More formally, for some prespecified turn K ≥ 0, we set

f ⋆ ∈ arg min f ∈F ,K (f )≤K E[∆(f , X)],
where K (f ) denotes the turn of f . Proceeding as in Section 2, we let K be a countable subset of [0, K] and define a countable collection of models {F k,κ } k≥1,κ∈K as follows. Each F k,κ consists of polygonal lines with k segments, with turn at most κ, and with vertices belonging to some grid Q over C. For k ≥ 1 and κ ∈ K, define fk,κ ∈ arg min

f ∈F k,κ ∆ n (f )
to be a polygonal line minimizing the empirical criterion ∆ n (f ) over F k,κ . We wish to design an appropriate penalty function pen :

N * × K → R + and minimize the criterion crit(k, κ) = ∆ n ( fk,κ ) + pen(k, κ)
in order to obtain a suitable principal curve. As before, we let f = fk ,κ , where ( k, κ) is a minimizer of the penalized criterion crit(k, κ), and intend to control the loss

D(f ⋆ , f ) = E[∆( f , X) -∆(f ⋆ , X)].
To get a result of the form of Theorem 2.2, we already know that it suffices to find an upper bound on the quantity

E sup f ∈F k,κ E[∆(f , X)] -∆ n (f ) .
As a first step towards this direction, we will need the following lemma, which establishes an interesting link between the length of a curve and its turn.

ζ(x) =                1 cos(x/2) if 0 ≤ x ≤ π 2 2 sin(x/2) if π 2 ≤ x ≤ 2π 3 x 2 - π 3 + √ 3 if x ≥ 2π 3 .
The graph of the function ζ is shown in Figure 5.

E sup f ∈F k,κ E[∆(f , X)] -∆ n (f ) ≤ δ 2 a 1 √ k + a 2 ζ(κ) + a 3 ζ(κ) √ k 1 { ζ(κ) 3k <1} + a 4 k ln ζ(κ) k 1 { ζ(κ) 3k ≥1} + a 0 .
Putting finally Theorem 2.1 and Proposition 3.1 together, we obtain:

(k, κ) ∈ N * × K, pen(k, κ) ≥ δ 2 √ n c 1 √ k + c 2 max ζ(κ) √ k , k ln ζ(κ) k + c 0 + x k,κ 2 , then E[D(f ⋆ , f)] ≤ inf k≥1,κ∈K D(f ⋆ , F k,κ ) + pen(k, κ) + δ 2 Σ 2 3/2 π n ,
where

D(f ⋆ , F k,κ ) = inf f ∈F k,κ D(f ⋆ , f ).
The expression of the penalty shape involves a term of the order k/n-just like in the case of curves with bounded length-, whereas the length ℓ is replaced by ζ(κ), which is an increasing function of the turn κ. This is relevant, since the number of segments k and the turn κ characterize the complexity of the models. Moreover, the additive term max ζ(κ)/ √ kn, k ln ζ(κ)/kn shows that k and κ should be cleverly chosen relatively to each other in order to get a nice principal curve. Roughly, a greater curvature implies more segments.

Experimental results

This section presents some simulations and real data experiments, carried out with the software MATLAB, to illustrate the model selection procedure suggested by Theorem 2.2. The penalty shapes in the theorem involve constants which have to be practically determined. To this end, a possible route is to use the so-called slope heuristics, introduced by Birgé and Massart [START_REF] Birgé | Minimal penalties for Gaussian model selection[END_REF] and further developed by Arlot and Massart [START_REF] Arlot | Data-driven calibration of penalties for least-squares regression[END_REF] (see also Lerasle [34], Saumard [START_REF] Saumard | The slope heuristics in heteroscedastic regression[END_REF], and the overview by Baudry, Maugis and Michel [START_REF] Baudry | Slope heuristics: overview and implementation[END_REF]). In short, this calibration method allows to tune a penalty known up to some multiplicative constant. The slope heuristics assumes that the empirical contrast decreases when the complexity of the models increases, which is clearly the case in our principal curve context. The procedure is based on the fact that the graph of the empirical contrast as a function of the penalty shape decreases strongly at the beginning and more slowly later, with a linear trend. At the end, the heuristics specifies that the desired constant is equal to -2s, where s is the slope of this line. Our approach consists in adapting this method to the bivariate case.

Hence, in the sequel, the number k of segments and the length ℓ of the principal curve are chosen according to the following strategy, denoted hereafter by MS:

Algorithm MS 1.
For each number k of segments, k = 1, . . . , 80, and for a range of values of the length ℓ, compute fk,ℓ by minimizing the empirical criterion ∆ n (f ) and record

∆ n ( fk,ℓ ) = 1 n n i=1 ∆( fk,ℓ , X i ).
2. Set x k,ℓ = 2 ln n and consider a penalty of the form

pen(k, ℓ) = c 1 √ k + c 2 ℓ.
3. Select the constants ĉ1 and ĉ2 using a bivariate version of the slope heuristics.

4. Retain the curve fk , l obtained by minimizing the penalized criterion

crit(k, ℓ) = ∆ n ( fk,ℓ ) -2(ĉ 1 √ k + ĉ2 ℓ).
Throughout this experimental section, the maximal values of the parameters have been chosen to be reasonably large without increasing the computation time uselessly. The maximal length and the step defining the range of values of ℓ depend on the scale of the considered data set.

The minimization of the criterion ∆ n (f ) (step 1 of the algorithm) is achieved through a MATLAB optimization routine.

The weights x k,ℓ were all set to 2 ln n. We realize that this choice is somewhat arbitrary. However, as mentioned in the discussion after Theorem 2.2, as soon as the number of models is not larger than n 2 -which is clearly the case in our examples-, this is a convenient choice, which, moreover, does not modify the penalty shape. Besides, the calibration of c 0 is a challenging question which has been given little consideration in the literature so far, even in the standard slope heuristics context. Note that a possible route to take the constant term into account was proposed by Lebarbier [START_REF] Lebarbier | Detecting multiple change-points in the mean of gaussian process by model selection[END_REF]. Here, in our bivariate framework, we deal, to a first approximation, with a penalty of the form

c 1 √ k + c 2 ℓ.
Assessing the values of the constants via the slope heuristics rests upon the assumption that, for large values of k and ℓ, ∆ n ( fk,ℓ ) behaves like c 1 √ k + c 2 ℓ. The constants ĉ1 and ĉ2 are then chosen via an ordinary least square regression and we compute the corresponding R 2 coefficient to measure the quality of the regression. We also tried a robust regression, whose results were observed to be very similar, and thus, are not reported here.

Finally, the results of the algorithm MS were systematically compared to the outputs of the Polygonal Line Algorithm of Kégl, Krzyżak, Linder and Zeger [START_REF] Kégl | Learning and design of principal curves[END_REF]. In short, this procedure optimizes the vertices of the curve one after the other, using a local version of the criterion ∆ n (f ), which relies on a local angle penalty. To our knowledge, this heuristic technique is not supported by any theoretical result. However, it is known to perform well and should in our context be understood as a benchmark.

Simulated data

In this first series of experiments, we considered two-dimensional data distributed with some noise around a reference curve. More formally, observations were generated from the model

X = Y + ε,
where Y is uniformly distributed over some planar curve f and ε is a bivariate Gaussian noise, independent of Y. Even if the generative curve f is not a principal curve stricto sensu-because of the model bias-, this Gaussian model is considered as a benchmark for simulations in the literature on principal curves.

The union of the generative curve and the estimated curve can be seen as a selfintersecting polygon, the area of which may be used to quantitatively asses how far the estimated curve is from the true one. In the sequel, we compute for each simulated example an error criterion corresponding to the average area over 20 trials, normalized with respect to the scale of the data.

In a first example, we let f be a half-circle with radius 1. The noise variance is set to 0.004 and the number n of observations to 100 (see Figure 6). Recall that the algorithm MS computes the criterion ∆ n ( fk,ℓ ) for a table of values of √ k and ℓ and selects the best constants according to a bivariate slope heuristics. Figure 7 shows the contour plot of ∆ n ( fk,ℓ ) as a function of √ k and ℓ, which supports the idea that this function is linear in √ k and ℓ when k and ℓ become large. The irregularities reflect the fact that the criterion ∆ n ( fk,ℓ ) is not decreasing continuously when increasing the parameters, though decreasing on the whole. This phenomenon, which also appears in the Polygonal Line Algorithm (PL hereafter), is due to a convergence problem related to the optimization function.

[A]

[B] It can be noted that the outputs of both algorithms have approximately the same quality, despite a few irregularities on the MS principal curve, not visible on the PL result. The methods MS and PL were also tested on a larger sample n = 250. We observed that both principal curves obtained with this sample size are very accurate.

[A] [B] [C] [D]
In a second set of numerical examples, we took handwritten-type digits as generative curves, with noise variance 0.04. As depicted in Figure 10, 150 observations were sampled around the digit 2 and the digit 3 and 250 observations around the digit 5.

[A]

[B]

[C] With respect to the digit 2 data, the MS principal curve follows the observations more closely than what is expected. On the other hand, the PL output looks smoother, but a comparison with the generative curve shows that the loop at the top and the angle at the bottom of the digit 2 are not recovered precisely. For the digit 3, we note again that the algorithm MS slightly overfits the data, whereas the smoother curve PL misses the angle. The same comment holds for the digit 5, but to a lesser degree. On this last example, both algorithms performed quite similarly and the resulting principal curves are visually satisfactory. For these simulated digits, the average R 2 coefficients equal respectively 0.87, 0.91 and 0.93 and the error areas 0.032, 0.026 and 0.021.

This small simulation study reveals, as expected, that a good automatic choice of the parameters k and l is crucial to obtain a suitable principal curve. On the whole, the visual quality of MS is fully acceptable, even if the principal curves fitted by this algorithm often follow the data quite closely, in particular when the sample size is not very large. In return, the global shape of the digit is better recovered than using PL.

Real data sets 4.2.1 NIST database digits

The first real-life data set used in this second series of experiments originated from NIST Special Database 19 (http://www.nist.gov/srd/nistsd19.cfm), containing handwritten characters from 3600 writers. The data consists in binary images scanned at 11.8 dots per millimeter (300 dpi), which uniformly fill the area corresponding to the thickness of the pen stroke. Skeletonization, which consists in reducing foreground regions in such an image without affecting the general shape of the handwritten character, often constitutes a preliminary step to perform character recognition (see, e.g., Deutsch [START_REF] Deutsch | Preprocessing for character recognition[END_REF] and Alcorn and Hoggar [START_REF] Alcorn | Preprocessing of data for character recognition[END_REF]).

Algorithms MS and PL were applied to the three NIST database digits visible in Figure 11.

[A]

[B]

[C] We observe that both results for the digit 2 are similar and completely satisfactory.

Regarding the digit 5, MS seems to better recover the overall shape. Finally, the principal curve fitted by MS for the digit 7 is suitable, whereas the output of PL looks definitely not satisfactory. In the MS regression step, we obtained R 2 coefficients equal to 0.99, 0.99 and 0.95 respectively.

As a general conclusion on these NIST digit data sets, we found that MS performs well. Here, the algorithm does not seem to overfit, probably because the sample size is large enough.

Seismic data

Together with satellite images, the localization of earthquakes is an essential source of information in geology for the study of seismic faults, whether in accretion or subduction regions. As an illustration, Figure 12 depicts seismic impacts in the world-the map is drawn using Miller's projection-, as well as a world map from the USGS (United States Geological Survey) showing the various lithospheric plates. The data set, which can be downloaded on the USGS website (http://earthquake.usgs.gov/research/data/centennial.php), is part of the "Centennial Catalog", listing the major earthquakes registered since 1900 (Engdahl et Villaseñor [START_REF] Engdahl | Global seismicity[END_REF]). In this subsection, we employ algorithm MS as a means to recover the borders of lithospheric plates using the earthquake localization data of Figure 12. Again, the PL output is given as a benchmark.

We decided to focus on two particularly representative seismic active zones. The first one (Z1 hereafter) is located in the Atlantic Ocean, to the west of the African continent (about 60 • S 50 • W to 40 • N 0 • ), and the second one (Z2 hereafter) extends from the south of Africa to the south of Australia (about 65

• S 0 • to 25 • S 160 • E).
The localization of these two regions on the world map is visible in Figure 13. The results for Z1 are shown in Figure 14 and for Z2 in Figure 15. In Figure 14, we see, for the seismic zone Z1, that the method MS again yields a principal curve following the data points quite closely. On the contrary, the PL algorithm provides a smoother curve, which at first sight seems a better result. However, the border of the lithospheric plate is probably more likely to look like the more irregular MS principal curve, as suggested by Figure 12 [B]. The same observation holds for Z2 (Figure 15). Moreover, in this case, the PL output does not recover the shape of the plate border, which certainly passes through the most northern points and not several degrees south. Apparently, the local penalty on the angles leads here to overpenalization. Thus, on this seismic data set, MS results seem to be more relevant.

It is noteworthy that using this type of earthquake data to draw faults could be especially useful to locate some faults which cannot be easily spotted and necessitate monitoring for seismic risk prevention. With this respect, Harding and Berghoff [START_REF] Harding | Fault scarp detection beneath dense vegetation cover: airborne lidar mapping of the Seattle fault zone, Bainbridge Island, Washington State[END_REF], employing a method based on airborne laser mapping, study for instance seismic hazards in a zone densely covered by vegetation, located in the Puget Lowland of Washington State, USA. Using a principal curve approach to solve this kind of problems is undoubtedly an interesting project for future research.

Proofs

Proof of Theorem 2.1

Theorem 2.1 is an adaptation of Theorem 8.1 in Massart [START_REF] Massart | Concentration Inequalities and Model Selection[END_REF]. We first recall the following lemma, which is a consequence of McDiarmid's inequality [START_REF] Mcdiarmid | On the method of bounded differences[END_REF] (see Massart [START_REF] Massart | Concentration Inequalities and Model Selection[END_REF]Theorem 5.3]).

Lemma 5.1. If X 1 , . . . , X n are independent random variables and G is a finite or countable class of real-valued functions such that a ≤ g ≤ b for every function g ∈ G, then, setting Z = sup g∈G n i=1 (g(X i ) -E[g(X i )]), we have, for every ε ≥ 0,

P {Z -E[Z] ≥ ε} ≤ exp - 2ε 2 n(b -a) 2 .
Proof of the theorem. Let ∆n (f ) = ∆ n (f ) -E[∆(f , X)] denote the centered empirical process. For all k ≥ 1 and ℓ ∈ L, for any f k,ℓ ∈ F k,ℓ , we have, by definition of f ,

∆ n ( f) + pen( k, l) ≤ ∆ n (f k,ℓ ) + pen(k, ℓ).
Equivalently,

∆ n ( f ) -∆ n (f k,ℓ ) ≤ pen(k, ℓ) -pen( k, l). Since ∆ n ( f ) = E[∆( f , X)] + ∆n ( f ) and ∆ n (f k,ℓ ) = E[∆(f k,ℓ , X)] + ∆n (f k,ℓ ), this inequality becomes E[∆( f , X)] -E[∆(f k,ℓ , X)] ≤ ∆n (f k,ℓ ) -∆n ( f ) + pen(k, ℓ) -pen( k, l). (4)
Moreover, for every f ∈ F ,

D(f ⋆ , f ) = E[∆(f , X) -∆(f ⋆ , X)],
so that

E[∆( f , X)] -E[∆(f k,ℓ , X)] = D(f ⋆ , f ) -D(f ⋆ , f k,ℓ ). (5)
Therefore, combining (4) and ( 5),

D(f ⋆ , f ) ≤ D(f ⋆ , f k,ℓ ) + ∆n (f k,ℓ ) -∆n ( f) + pen(k, ℓ) -pen( k, l). (6)
Consider now a family of nonnegative weights {x k,ℓ } k≥1,ℓ∈L such that k≥1,ℓ∈L e -x k,ℓ = Σ < ∞, and let z > 0. Applying Lemma 5.1, we get, for all k ′ ≥ 1, ℓ ′ ∈ L and ε ≥ 0,

P    sup f ∈F k ′ ,ℓ ′ (-∆n (f )) ≥ E sup f ∈F k ′ ,ℓ ′ (-∆n (f )) + ε    ≤ exp - 2nε 2 δ 4 .
This may be rewritten, for

ε = δ 2 x k ′ ,ℓ ′ +z 2n , P    sup f ∈F k ′ ,ℓ ′ (-∆n (f )) ≥ E sup f ∈F k ′ ,ℓ ′ (-∆n (f )) + δ 2 x k ′ ,ℓ ′ + z 2n    ≤ e -x k ′ ,ℓ ′ -z . Setting E k ′ ,ℓ ′ = E sup f ∈F k ′ ,ℓ ′ (-∆n (f ))
, we thus have, for all k ′ ≥ 1 and ℓ ′ ∈ L,

sup f ∈F k ′ ,ℓ ′ (-∆n (f )) ≤ E k ′ ,ℓ ′ + δ 2 x k ′ ,ℓ ′ + z 2n ,
except on a set of probability not larger than Σe -z . Then, inequality [START_REF] Bartlett | Model selection and error estimation[END_REF] implies

D(f ⋆ , f) ≤ D(f ⋆ , f k,ℓ ) + ∆n (f k,ℓ ) + E k, l + δ 2 x k, l + z 2n -pen( k, l) + pen(k, ℓ) ≤ D(f ⋆ , f k,ℓ ) + ∆n (f k,ℓ ) + E k, l + δ 2 x k, l 2n -pen( k, l) + pen(k, ℓ) + δ 2 z 2n ,
except on a set of probability not larger than Σe -z . Consequently, if for all k ′ ≥ 1 and ℓ ′ ∈ L,

pen(k ′ , ℓ ′ ) ≥ E k ′ ,ℓ ′ + δ 2 x k ′ ,ℓ ′ 2n , then D(f ⋆ , f) ≤ D(f ⋆ , f k,ℓ ) + ∆n (f k,ℓ ) + pen(k, ℓ) + δ 2 z 2n ,
except on a set of probability not larger than Σe -z . Put differently,

P δ -2 √ 2n[D(f ⋆ , f) -D(f ⋆ , f k,ℓ ) + ∆n (f k,ℓ ) + pen(k, ℓ)] ≥ √ z ≤ Σe -z , or, letting z = u 2 , P [δ -2 √ 2n[D(f ⋆ , f ) -D(f ⋆ , f k,ℓ ) + ∆n (f k,ℓ ) + pen(k, ℓ)] ≥ u ≤ Σe -u 2 .
Recalling that ∞ 0 e -u 2 du = √ π 2 and letting g + = max(g, 0), we obtain

E (D(f ⋆ , f) -D(f ⋆ , f k,ℓ ) + ∆n (f k,ℓ ) + pen(k, ℓ)) + ≤ δ 2 Σ 2 3/2 π n .
Hence, as E[ ∆n (f k,ℓ )] = 0,

E[D(f ⋆ , f)] ≤ D(f ⋆ , f k,ℓ ) + pen(k, ℓ) + δ 2 Σ 2 3/2 π n .
Since this is true for all k and ℓ, we finally get

E[D(f ⋆ , f)] ≤ inf k≥1,ℓ∈L D(f ⋆ , F k,ℓ ) + pen(k, ℓ) + δ 2 Σ 2 3/2 π n , where D(f ⋆ , F k,ℓ ) = inf f ∈F k,ℓ D(f ⋆ , f ).
This concludes the proof of Theorem 2.1.

Proof of Proposition 2.1

The first step consists in proving that the quantity

E sup f ∈F k,ℓ (E[∆(f , X)] -∆ n (f ))
may be upper bounded by means of the Rademacher average

E sup f ∈F k,ℓ 1 n n i=1 ε i ∆(f , X i ) ,
where ε 1 , . . . , ε n are independent Rademacher random variables, defined by P {ε i = 1} = P {ε i = -1} = 1/2, independent of X 1 , . . . , X n . Let X ′ 1 , . . . , X ′ n be independent copies of X 1 , . . . , X n , also independent of ε 1 , . . . , ε n . A symmetrization argument yields

E sup f ∈F k,ℓ (E[∆(f , X)] -∆ n (f )) = E sup f ∈F k,ℓ E 1 n n i=1 ∆(f , X ′ i ) X 1 , . . . , X n - 1 n n i=1 ∆(f , X i ) ≤ E sup f ∈F k,ℓ 1 n n i=1 ∆(f , X ′ i ) -∆(f , X i ) = E sup f ∈F k,ℓ 1 n n i=1 ε i ∆(f , X ′ i ) -∆(f , X i ) ≤ E sup f ∈F k,ℓ 1 n n i=1 ε i ∆(f , X ′ i ) + E sup f ∈F k,ℓ 1 n n i=1 (-ε i )∆(f , X i ) = 2E sup f ∈F k,ℓ 1 n n i=1 ε i ∆(f , X i ) .
Next, the Rademacher average

E sup f ∈F k,ℓ 1 n n i=1 ε i ∆(f , X i )
may be bounded by resorting to a Dudley integral. More precisely, let

S k,ℓ = {∆(f , •), f ∈ F k,ℓ }
be a subset of the continuous functions from C to R + , endowed with the sup-norm • ∞ , and denote by N (S k,ℓ , • ∞ , ε) the covering number of S k,ℓ , i.e., the minimal number of closed balls of radius ε needed to cover S k,ℓ . According to Dudley [START_REF] Dudley | Uniform Central Limit Theorems[END_REF], there exists an absolute constant c > 0 such that, for all X 1 , . . . , X n ,

E sup f ∈F k,ℓ 1 n n i=1 ε i ∆(f , X i ) ≤ c √ n ∞ 0 ln N (S k,ℓ , • ∞ , ε)dε.
To evaluate the covering number of S k,ℓ , we may use Lemma 2 in Kégl [START_REF] Kégl | Principal Curves: Learning, Design, and Applications[END_REF], which ensures that

N (S k,ℓ , • ∞ , ε) ≤ 2 ℓδ/ε+3k+1 V k+1 d δ 2 √ d ε + √ d d ℓδ √ d kε + 3 √ d kd ,
where V d denotes the volume of the d-dimensional unit ball. Observe that

ln N (S k,ℓ , • ∞ , ε) ≤ ℓδ ε + 3k + 1 ln 2 + (k + 1) ln V d + d ln δ 2 √ d ε + √ d + kd ln ℓδ √ d kε + 3 √ d = ℓδ ε ln 2 + (3k + 1) ln 2 + (k + 1) ln V d + d(k + 1) ln √ d + d ln δ 2 ε + 1 + kd ln 3 + kd ln ℓδ 3kε + 1 = ℓδ ε ln 2 + d ln δ 2 ε + 1 + kd ln ℓδ 3kε + 1 + kd ln 3 + (3k + 1) ln 2 + (k + 1)(ln V d + d 2 ln d).
Hence, recalling that the support of f is included in a set C with diameter δ, we obtain

∞ 0 ln N (S k,ℓ , • ∞ , ε)dε = δ 2 0 ln N (S k,ℓ , • ∞ , ε)dε ≤ √ ℓδ ln 2I 1 + √ dI 2 + √ kdI 3 + δ 2 A(k, d),
where

I 1 = δ 2 0 1 √ ε dε, I 2 = δ 2 0 ln δ 2 ε + 1 dε, I 3 = δ 2 0 ln ℓδ 3kε + 1 dε, and
A(k, d) = kd ln 3 + (3k + 1) ln 2 + (k + 1)(ln V d + d 2 ln d) 1/2 .
Control of I 1 . Clearly,

I 1 = δ 2 0 1 √ ε dε = 2δ.
Control of I 2 . We have

I 2 ≤ δ 2 0 ln 2δ 2 ε dε = 2δ 2 1/2 0 ln 1 u du ≤ δ 2 ( √ ln 2 + √ π).
Control of I 3 . Let M = max(3k, L/δ). Clearly, for all ℓ ∈ L, δ ≥ ℓ M , and then δ 2 ≥ ℓδ M . Let us cut up the integral I 3 and write

I 3 = δ 2 0 ln ℓδ 3kε + 1 dε = ℓδ/M 0 ln ℓδ 3kε + 1 dε + δ 2 ℓδ/M ln ℓδ 3kε + 1 dε. (7)
Observe, since ε ≤ ℓδ M , that ℓδ 3kε ≥ 1. Consequently,

ℓδ/M 0 ln ℓδ 3kε + 1 dε ≤ ℓδ/M 0 ln 2ℓδ 3kε dε = 2ℓδ 3k 3k/2M 0 ln 1 u du ≤ ℓδ M ln 2M 3k + √ π .
The second integral in equality ( 7) may be bounded using the fact that the integrand is a decreasing function of ε:

δ 2 ℓδ/M ln ℓδ 3kε + 1 dε ≤ δ 2 - ℓδ M ln M 3k + 1 ≤ δ 2 - ℓδ M ln 2M 3k .
As a result, where the nonnegative constants a 0 , . . . , a 2 only depend on the maximal length L, the dimension d and the diameter δ of the convex set C.

I 3 ≤ δ 2 ln 2M 3k + ℓδ M √ π. Thus, ∞ 0 ln N (S k,ℓ , • ∞ , ε)dε ≤ 2δ √ δℓ

Proof of Proposition 3.1

Let S k,κ = {∆(f , •), f ∈ F k,κ } be a subset of the continuous functions from C to R + , endowed with the sup-norm • ∞ . Starting as in the proof of Proposition 2.1, we know that, for all X 1 , . . . , X n ,

E sup f ∈F k,κ E[∆(f , X)] -∆ n (f ) ≤ c √ n ∞ 0 ln N (S k,κ , • ∞ , ε)dε,
for some absolute constant c > 0. Now, according to Lemma 5 in Sandilya and Kulkarni [START_REF] Sandilya | Principal curves with bounded turn[END_REF], we may write, for each ε > 0, Control of I 1 . We clearly have

I 1 = δ 2 0 1 √ ε dε = 2δ.
Control of I 2 . We have

I 2 ≤ δ 2 0 ln 2δ 2 ε dε = 2δ
The first integral is bounded by using the inequality ζ(κ)δ 2 3kε ≥ 1 for all ε ∈ ]0,

≤ δ 2 a 1 √ k + a 2 ζ(κ) + a 3 ζ(κ) √ k 1 { ζ(κ) 3k <1} + a 4 k ln ζ(κ) k 1 { ζ(κ) 3k ≥1} + a 0 ,
where the nonnegative constants a 0 , . . . , a 4 only depend on the dimension d.

Figure 1 :

 1 Figure 1: An example of principal curve.

Figure 3 :

 3 Figure 3: Principal curves fitted with [A] a too small number k of segments, [B] a too large k and [C] an appropriate one.

Figure 4 :

 4 Figure 4: Denoting by -→ s j the vector ---→ v j v j+1 for all j = 1, . . . , k, the angles involved in the definition of the turn are defined by φ j+1 = ( -→ s j , --→ s j+1 ).

Lemma 3 . 1 .

 31 For a proof of this result, we refer the reader to Alexandrov and Reshetnyak [2, Chapter 5]. Let f be a curve with turn κ and let δ be the diameter of C. Then L (f ) ≤ δζ(κ), where the function ζ is defined by

Figure 5 :Proposition 3 . 1 .

 531 Figure 5: Graph of the function ζ.

Theorem 3 . 1 .

 31 Consider a family of nonnegative weights {x k,κ } k≥1,κ∈K such that k≥1,κ∈K e -x k,κ = Σ < ∞, and a penalty function pen : N * × K → R + . Let f = fk ,κ . There exist nonnegative constants c 0 , . . . , c 2 , depending only on the dimension d, such that, if for all

Figure 6 :

 6 Figure 6: 100 observations distributed around a half-circle with radius 1.

Figure 7 :

 7 Figure 7: [A] Contour plot of ∆ n ( fk,ℓ ) as a function of √ k and ℓ for the half-circle data (n=100). [B] Zoom in the zone of large values of k and ℓ.

Figure 8 :

 8 Figure 8: Selected principal curves for the half-circle data (n=100). [A] Method MS: k = 20, l = 3. [B] PL algorithm.

Figure 9 :

 9 Figure 9: Method MS: Examples of principal curves for some prespecified values of k and ℓ (n=100). [A] k = 20, ℓ = 2.5. [B] k = 20, ℓ = 3.1. [C] k = 20, ℓ = 3.4. [D] k = 35, ℓ = 4.

Figure 10 :

 10 Figure 10: [A] Observations sampled around the digit 2 (n=150), 3 (n=150) and 5 (n=250). [B] Principal curves selected by the method MS: k = 27, l = 24; k = 23, l = 23; k = 17, l = 21. [C] PL principal curves.

Figure 11 :

 11 Figure 11: [A] Three NIST database handwritten digits. [B] Principal curves selected by the method MS: k = 23, l = 80; k = 38, l = 82; k = 15, l = 66. [C] PL principal curves.

Figure 13 :

 13 Figure 13: Localization of the two considered seismic zones Z1 (about 60 • S 50 • W to 40 • N 0 • ) and Z2 (about 65 • S 0 • to 25 • S 160 • E).

Figure 14 :

 14 Figure 14: Selected principal curves for the seismic zone Z1 (n=252). [A] Method MS: k = 55, l = 31. [B] PL principal curve.

Figure 15 :

 15 Figure 15: Selected principal curves for the seismic zone Z2 (n=322). [A] Method MS: k = 22, l = 38. [B] PL principal curve.

2 0

 2 S k,ℓ , • ∞ , ε)dε = δ ln N (S k,ℓ , • ∞ , ε)dε ≤ δ ζ(κ) ln 2I 1 + √ dI 2 + √ kdI 3 + δ 2 A(k, d), , d) = δ 2 kd ln 3 + (2k + 1) ln 2 + (k + 1)(ln V d +

  ln N (S k,κ , • ∞ , ε)

	≤ =	ζ(κ)δ 2 ε + kd ln ζ(k)δ 2 ε ln 2 + (2k + 1) ln 2 + (k + 1) ln V d + d(k + 1) ln + 2k + 1 ln 2 + (k + 1) ln V d + d ln δ 2 √ d + ε ζ(κ)δ 2 √ d + 3 √ d kε	√ √ d + d ln d	δ 2 ε	+ 1
		+ kd ln 3 + kd ln	ζ(κ)δ 2 3kε	+ 1
	=	ζ(k)δ 2 ε	ln 2 + d ln		δ 2 ε	+ 1 + kd ln	ζ(κ)δ 2 3kε	+ 1 + kd ln 3 + (2k + 1) ln 2
		+ (k + 1)(ln V d +	d 2	ln d).

  With respect to the second integral in[START_REF] Birgé | From model selection to adaptive estimation[END_REF], we note that the function under the integral is decreasing in ε, so that + lnV d + 2 ln 2) 1/2 + a 0 ,where a 0 is a nonnegative constant. Finally,

	Hence, collecting the different results,
	0	∞	ln N (S k,ℓ , • ∞ , ε)dε
	≤ 2δ 2 ζ(κ) ln 2 +	√	dδ 2 (	√	ln 2 +	√	π) + δ 2 √	kd ln	ζ(κ) 3k	+	√	π +	√	ln 2 1 { ζ(κ) 3k ≥1}
		+ δ 2 √	kd	ζ(κ) 3k	√	π +	√	ln 2 1 { ζ(κ) 3k <1} + δ 2 A(k, d)
	≤ δ 2 2 ζ(κ) ln 2 + + √ k √ d( √ π + √ ζ(κ) 3 √ k ln 2) + d ln 3 + √ πd1 { ζ(κ) 3k <1} + kd ln d 2 ln d	ζ(κ) 3k	1 { ζ(κ) 3k ≥1}
	We obtain											ζ(κ)δ 2 3k ].
				0	ζ(κ)δ 2 /3k		ln		ζ(κ)δ 2 3kε	+ 1 dε ≤	0	ζ(κ)δ 2 /3k	ln	2ζ(κ)δ 2 3kε	dε
															=	2ζ(κ)δ 2 3k	0	1/2	ln	1 u	du
															≤	ζ(κ)δ 2 3k	(	√	ln 2 +	√	π).
					δ 2 ζ(κ)δ 2 /3k	ln	ζ(κ)δ 2 3kε	+ 1 dε ≤ δ 2 -	ζ(κ)δ 2 3k	√	ln 2.
	Thus, we have
					I 3 ≤	           δ 2	ln δ 2 ζ(κ) ζ(κ) 3k 3k √ π + + √ π + √ ln 2 √	ln 2 if if	ζ(κ) 3k ζ(κ) 3k	≥ 1 < 1.

∞ 0 ln N (S k,ℓ , • ∞ , ε)dε

Acknowledgments. The authors thank three referees and the Associate Editor for valuable comments and insightful suggestions.

 [START_REF] Alcorn | Preprocessing of data for character recognition[END_REF]Research partially supported by the French "Agence Nationale pour la Recherche" under grant ANR-09-BLAN-0051-02 "CLARA". 2 Research carried out within the INRIA project "CLASSIC" hosted by Ecole Normale Supérieure and CNRS.


Jean-Daniel Kant
email: jean-daniel.kant@lip6.fr

Alain Lifchitz
email: alain.lifchitz@lip6.fr

Web-R: a tool to record & replay personal web navigation

Keywords: Personal Web navigation, Personal Web archive, Web history search, Web accessibility, Offline browsers, Tools and techniques for Web personalization, Web visualization, Information management tool

This poster presents a useful tool to capture the content of browsing sessions. Web-R saves systematically all the components sufficient and necessary to visualize offline the page, the same way it was displayed online. In order to avoid the saturation of local disk, Web-R integrates a page comparison mechanism to avoid unnecessary redundancy and, in addition, provides a way to manage the maximum storage space. We give some insights than, in typical cases, this storage is compatible with today's casual hard disk capacities and costs, along years. Moreover, since the full content of visited pages is stored, the Web-R user can access and process (e.g. sort, extract) several information including path and visit statistics, therefore giving him/her a global view on his/her personal navigation.

INTRODUCTION

In this work, we want to design a companion recording tool of the browser, thus integrating navigation, browsing and recording, whether they occur online or offline. We propose to save automatically the pages the user already browsed in order to help this user to organize his/her own Web. In order to be effective, the recording tool has to be non-intrusive (Figure 1 below), fast and should provide space management services in order to avoid possible local disk saturation. It should also provide an exact and complete storage of user's web pages and web navigation, by saving all the components sufficient and necessary to visualize offline the page the same way it was displayed online. Most of existing tools (online and/or offline browsers [START_REF] Surfsaver | [END_REF], site rippers [6], recorders [START_REF] Anupam | Automating Web Navigation with the WebVCR[END_REF]) do not fulfill our requirements since they are usually based on URL recording that cannot capture the dynamic and sometimes transient nature of web pages. Moreover, these approaches lack automation and user-centric (no adaptation, works only at request since the user has to ask to record the page each time he/she browses it).

DESIGN

Web-R is an application that will take control of its child running instance of Microsoft Internet Explorer ® (IE) using OLE Automation ® . Only the instances that have been created through Web-R will be monitored.

Page recording

The main architecture is depicted in Figure 2 The Pilot module is the heart of Web-R. It controls the configuration, tools and database modules, and a list of IE instances. Each IE instance in this list is coupled with a Coop process, which is used to communicate between this instance and Pilot. IE instances are independent, so the user could browse several pages simultaneously.

The page recording proceeds as depicted in Figure 3 below. The IE instance retrieves the page from the Web (1), and interprets its code for display. While the page is downloaded by IE, the Coop module monitors this download process, and waits until IE has finished its processing [START_REF] Keepoint | [END_REF]. Coop recuperates the result of this processing (i.e. the processed page) and sends it to the Page Source Processor (PSP) module [START_REF] Milic-Frayling | WebScout: Web Navigation Aid and Personal Web History Explorer[END_REF]. The aim of this module is to perform all the additional processing that is necessary to ensure a perfect storage of the current page. Therefore, the PSP might collect additional files (e.g. images, sounds, Flash Macromedia code, etc.) from the Web (4). Once the page processing is fully completed, the page and its components are sent to the Pilot module [START_REF] Surfsaver | [END_REF] to store all the files on the database (6). The main task of the PSP is to handle properly file name redirection of live Web to local hard drive. It handles the complete storage of frames, images, secured pages, animations (e.g. Flash Macromedia ® ), scripts and dynamic pages (e.g. ASP, PHP, CGI, etc.).

Storage

Page comparison algorithm

In order to avoid redundancy, we implemented a simple surface comparison algorithm to determine if the current page already exists in the database. This algorithm proceeds in 3 steps to determine if a page P needs to be stored: This algorithm is fast and will work for most of the situations encountered during realistic web navigation.

Storage space management

Since all visited pages will be stored on local disk during a Web-R session, one may be concerned about space on local disk. In fact, we believe that the total size of a personal web content will not exceed 5 GB / year for a typical user. Recent studies confirm this estimation ([4]). Anyway, we must provide a tool to manage storage space. Thus, the Web-R user can set a size limit of the storage space that could be none, a fixed size or a percent of the total disk size. When the limit is reached, a purge mechanism is started, which could be either manual or automatic.

Managing the personal web

In the current version, Web-R also offers a basic tool to manage the stored pages. When started, it displays a multicolumn list window that contains all the stored pages. The page entries can be sorted / filtered by the following attributes: URL, title, number of access, date and time of first visit, date and time of last visit, logical / physical size. The managing tool could also be used to visualize a stored page. Its (offline) content will be displayed when the user doubleclick on one item of the page list.

CONCLUSION AND FUTURE DIRECTIONS

In this work, we have shown that a systematic storage of the personal web is technically feasible and realistic since it will not require too much disk space.

Two very recent systems have very similar features to Web-R. Keepoint [START_REF] Keepoint | [END_REF] also provides an automatic and complete storage of web pages browsed by IE but does not integrate a page comparison mechanism hence requiring the user to decide whether a similar page should be overwritten. The Microsoft WebScout [START_REF] Milic-Frayling | WebScout: Web Navigation Aid and Personal Web History Explorer[END_REF] project seems to be very ambitious and its HistoryExplorer component shares the same goals as Web-R. However, we do not have any details on this implementation, so it is difficult to compare it with our proposed solution.

Web-R is currently under development and several technical issues need to be addressed. In the short run, we need to consolidate the page processor in order to solve the cases where the storage was not completed (typically dynamical pages and scripts). Moreover, we will facilitate offline navigation, where links to already stored pages (local) and links to non-stored pages (exterior) will be distinguished. In the latter case, the system may ask the user if he wants to locally store the exterior page or navigate this page online. We also plan to make a better use of the browser resources, like using the cache to avoid multiple -and therefore unnecessary -fetches of components on the live Web. We should also compress the stored files to optimize space storage, and encrypt all the data to ensure a total privacy.

Figure 1 .

 1 Figure 1. Web-R root interface (actual scale)

  below.

Figure 2 .

 2 Figure 2. Overall Web-R architecture

Figure 3 .

 3 Figure 3. Page recording flowchart

1 .

 1 Let SP a set of stored page with the same URL as P. If SP is empty, store P; 2. Otherwise, look for a page P ' SP with the same HTML code as P. If none, store P; 3. Otherwise, compare respective individual size of the components in P ' and P. If they differ, store P; otherwise do not store P.

ACKNOWLEDGEMENTS

We would like to thank the LIP6 for its financial support on this project. We also thank Nicolas Limare for his participation to the development of the preliminary version of the Web-R software.


Philippe Négrel

Emmanuelle Petelet-Giraud

Romain Millot

Eline Malcuit

Source and mobility of Rare Earth Elements in a sedimentary aquifer system: Aquitaine basin (Southern France)

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

Source and mobility of Rare Earth Elements in a sedimentary

aquifer system: Aquitaine basin (Southern France)

Ph. NEGREL 1 , E. PETELET-GIRAUD 1 , R. MILLOT 1 , E. MALCUIT 2 1 BRGM, Metrology Monitoring Analysis Department, Orléans, France p.negrel@brgm.fr, e.petelet@brgm.fr, r.millot@brgm.fr 2 BRGM, Regional Geological Survey Service, Bordeaux, France, e.malcuit@brgm.fr

The study of rare earth elements (REEs) in natural waters initially involved an examination of their occurrence and behavior in seawater and coastal waters such as estuaries. Since the 1990s, REE geochemistry has been applied to continental waters such as rivers and lakes and groundwaters. Rare earth elements) are of great interest because of their unique characteristics and have been used in the study of many geological processes like weathering and water-rock interaction processes, provenance of sediments, etc... With the evolution of analytical techniques like new generation ICP-MS, much attention had been paid towards the water geochemistry of REEs. However, there is a need of more investigations devoted to REEs in large groundwater systems, especially on the understanding of the distribution of REEs and their evolution in such systems. In this frame, large sedimentary aquifer systems often constitute strategic water resources for drinking water supply, agriculture irrigation and industry, but can also represent an energetic resource for geothermal power. Large water abstractions can induce complete modification of the natural functioning of such aquifer systems. These large aquifer systems thus require water management at the basin scale in order to preserve both water quantity and quality.

The large Eocene Sand aquifer system of the Aquitaine sedimentary basin was studied through various hydrological, chemical and isotopic tools. This system extends over 116,000 km² in the South west part of the French territory. The aquifer being artesian in the west of the district and confined with piezometric levels around 250-m depth in the east. The 'Eocene Sands', composed of sandy Tertiary sediments alternating with carbonate deposits, is a multi-layer system with high permeability and a thickness of several tens of metres to a hundred metres. The Eocene Sand aquifer system comprises at least five aquifers: Paleocene, Eocene infra-molassic sands (IMS), early Eocene, middle Eocene, and late Eocene. One important feature, in these confined systems isolated from anthropogenic influence, is the range in salinities by a factor of 10, from 250 mg/L up to 2.5 g/L. The ΣREE, in the range 2-54 ng/L, with a dependence on salinity when expressed in % HCO 3 or SO 4 , reflect the carbonate or evaporite source of REEs. The UCC normalized-REE patterns show a large variability as exemplified by the REE flat patterns-low SREE associated with salinity controlled by HCO 3 .

In the present work, the REEs are investigated in terms of saturation indices, speciation modelling, REE patterns in order to recognize the aquifer type hosting groundwater and decipher the origin of the salinity of the groundwater as some part of the aquifer display in the groundwater concentration of chemical element exceeding the drinking water standard (SO 4 , F…). Such high concentrations of naturally-occurring substances (e.g. unaffected by human activities) can have negative impacts on groundwater thresholds and deciphering their origin by means of geochemical tools like REE is a remaining challenge.


Manuela Campiglio

Rosaria Bufalino

Marco Sandri

Elisa Ferri

Rosa Anna Aiello

Andrea Matteis

Marcella Mottolese

Sabino Placido

Patrizia Querzoli

Andrea De Matteis

Sabino De Placido

Antonio Jirillo

Alberto Bottini

Manuela Fantini

Andrea Bonetti

Fulvia Pedani

Maria Mauri

Annamaria Molino

Antonella Ferro

Serenella M Pupa

Marianna Sasso

Sylvie Ménard
email: sylvie.menard@istitutotumori.mi.it

Andrea Balsari

Elda Tagliabue

Eleonora Taibi

Giuseppe Nuciforo

Gianfranco Di Feo

Giulio Metro

Federica Cuppone

Giuseppe Cancello

Maria Lunardi

Annamaria Scola

Alessandra Bersiga

Simone Bonardi

Lorenzo Gianni

Alberto Ravaioli

Marta Zaninelli

Laura Furini

Ester Berardengo

Leopoldo Costarelli

Eliana Gilioli

Michela Frisinghelli

Alberto Morabito

Mario Brandi

Giovanni Simone

Giovanna Cavazzini

Alberto Bellomi

Giorgio Mustacchi

Fabrizio Zanconati

Alessandro Bertolini

Salvatore Ambrosi

Angelo Di Leo

Claudino Wederson

Anna Sapino

Isabella Castellano

Dario Giuffrida

Lorenzo Memeo

Andrea Decensi

Increased overall survival independent of RECIST response in metastatic breast cancer patients continuing trastuzumab treatment: evidence from a retrospective study

Keywords: HER2, trastuzumab, metastatic breast cancer, RECIST, propensity score matching

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

INTRODUCTION

First line treatment of HER2 positive metastatic breast cancer (MBC) patients with the anti-HER2 monoclonal antibody trastuzumab is now well established [START_REF] Slamon | Use of chemotherapy plus a monoclonal antibody against HER2 for metastatic breast cancer that overexpresses HER2[END_REF][START_REF] Marty | Randomized phase II trial of the efficacy and safety of trastuzumab combined with docetaxel in patients with human epidermal growth factor receptor 2-positive metastatic breast cancer administered as first-line treatment: the M77001 study group[END_REF]. However, patients may develop resistance to trastuzumab and show metastatic disease progression. Several retrospective studies published thus far (reviewed in [START_REF] Fabi | Do HER-2 positive metastatic breast cancer patients benefit from the use of trastuzumab beyond disease progression? A mono-institutional experience and systematic review of observational studies[END_REF][START_REF] Cancello | Continuing trastuzumab beyond disease progression: outcomes analysis in patients with metastatic breast cancer[END_REF]) [START_REF] Extra | Efficacy of trastuzumab in routine clinical practice and after progression for metastatic breast cancer patients: the observational Hermine study[END_REF][START_REF] Montemurro | Retrospective evaluation of clinical outcomes in patients with HER2-positive advanced breast cancer progressing on trastuzumab-based therapy in the pre-lapatinib era[END_REF] have indicated a certain overall response rate, clinical benefit and extended time to progression from the use of a second trastuzumabcontaining regimen in HER2-positive MBC patients who progressed on first-line trastuzumab treatment. Moreover, two recent randomized trials designed to directly evaluate the efficacy of continuing trastuzumab after progression demonstrated the potential clinical utility of continuing trastuzumab even for patients with progressive disease [START_REF] Minckwitz | Trastuzumab beyond progression in human epidermal growth factor receptor 2-positive advanced breast cancer: a german breast group 26/breast international group 03-05 study[END_REF][START_REF] Blackwell | Randomized study of Lapatinib alone or in combination with trastuzumab in women with ErbB2-positive, trastuzumab-refractory metastatic breast cancer[END_REF]. To further evaluate the benefits of trastuzumab-continuation therapy we carried out a retrospective analysis of 272 HER2 positive MBC patients first line trastuzumab treated in the metastatic setting at different oncology centers in Italy, who continued or discontinued trastuzumab at progression. Since RECIST (Response Evaluation Criteria in Solid Tumors) was used for determining efficacy of trastuzumabcontaining regimens by the oncologists participating in this study, we specifically evaluated RECIST responses upon first line trastuzumab treatment, as a potential predictive marker for the therapeutic effects of trastuzumab-continuation beyond metastatic disease progression. The present study is the first report of increased overall survival independent of RECIST response achievement on initial chemotherapy regimens that include trastuzumab.

METHODS

Response evaluation

Clinical response of patients during first line trastuzumab treatment was evaluated according to RECIST 1.0 criteria [START_REF] Therasse | New guidelines to evaluate the response to treatment in solid tumors. European Organization for Research and Treatment of Cancer, National Cancer Institute of the United States, National Cancer Institute of Canada[END_REF] for the target lesions as follows: i) complete response (CR) based on disappearance of all lesions; ii) partial response (PR) based on at least 30% reduction of the sum of maximum lesion diameters; iii) progressive disease (PD) based on at least 20% increase in the sum of maximum lesion diameters; and iv) stable disease (SD) based on a variation of the sum of maximum target lesion diameters between PR and PD. Clinical response of non-target lesions was categorized as: complete response (CR), progression (PD) or noCR/noPD. Response to trastuzumab was evaluated based on ultrasound, total body scan, magnetic resonance imaging, computed axial tomography or other radiological examination. Patients were considered as responsive when they presented a complete or partial response in at least one metastatic lesion.

Independent central response review has been done at Fondazione IRCCS-Istituto Nazionale

Tumori. The institutional review board of each participant institution approved the observational study.

Statistical methods

To reduce potential bias in effects of treatment choices by physicians/patients or of differences in prognosis before treatment, either of which might result in large differences between treated and untreated groups in covariates, we used a propensity score [START_REF] Rosenbaum Pr Rubin | The central role of the propensity score in observational studies for causal effects[END_REF] to assess the conditional probability of treatment according to the individual's covariates and to balance treatment choicerelated variables such that the analysis simulates random assignment. We estimated a propensity score using a logistic regression model where the outcome variable is the oncologists' decision to continue/discontinue trastuzumab and the explanatory variables are the covariates selected according to their variable importance, estimated using random forests [START_REF] Breiman | Random Forests[END_REF]. Our random forest contained 30000 trees and 6 random features to split each node. Mean decrease in prediction accuracy was the estimated measure of variable importance.

Multivariate 1:1 matching with replacement performed by a 0.1 caliper matching on the estimated propensity score, together with an exact matching with site of metastases, generated 238+238 matched treated and controls units; 34 (12.5%) observations were dropped by matching. The degree of covariate imbalance was measured using the standardized (mean and proportion) difference proposed by Austin et al [START_REF] Austin | A comparison of the ability of different propensity score models to balance measured variables between treated and untreated subjects: a Monte Carlo study[END_REF].

Overall Survival (OS) of trastuzumab-treated patients was calculated as the time from the date of first documented progression of disease until death from any cause or until the end of the observation period. The hazard ratios related to trastuzumab treatment suspension after disease progression were estimated on matched data by a Cox proportional hazards model with random effects (frailty), stratified by investigational center [START_REF] Clayton | A model for association in bivariate life tables and its application in epidemiological studies of familial tendency in chronic disease incidence[END_REF]. Predictive covariates for survival were selected according to their variable importance, estimated using survival random forests [START_REF] Ishwaran | Random Survival Forests[END_REF].

The assumption of proportional hazards was verified using the test of Grambsch and Therneau [15], based on Schoenfeld residuals. The proportional hazards assumption was tested globally and for each covariate, with all null hypotheses set at a 5% level of significance.

Differences of treatment effects in the responder and non-responder subgroups were investigated considering in the Cox model the interaction between patient response to trastuzumab-based therapy and continuation/discontinuation of treatment and by testing (using a likelihood-ratio test) whether the coefficient of the interaction term was significantly different from zero. The survival curve of Figure 1 was estimated on unmatched data by the Kaplan-Meier method.

Figure 2 shows the predicted survivor functions estimated on matched data using the Cox proportional hazard model with frailty.

Missing values of the analyzed covariates were imputed by the method described in Iacus and Porro [START_REF] Iacus | Missing data imputation, matching and other applications of random recursive partitioning[END_REF] based on proximity of random forests. Imputed data were used only for variable selection and propensity score estimation. Matching and model estimation were applied to selected variables with missing values. The Cox proportional hazard model was estimated using 442 observations (34 were deleted due to missingness).

Statistical analyses were carried out using the R software for statistical computing, v. 2.8.1 (R Development Core Team, 2008), with the following packages: randomForest, randomSurvivalForest, survival and Matching.

RESULTS

To evaluate RECIST responses upon first line trastuzumab treatment, as a potential predictive marker for any therapeutic effect of trastuzumab-continuation beyond metastatic disease progression, we have examined the benefits of trastuzumab-continuation therapy by carrying out a retrospective analysis of 272 HER2 positive MBC patients under treatment at 22 different oncology centers in Italy during the years of 2000 and 2001 (Fig. 1). Our analysis included 154 (56.6%) progressed patients who continued trastuzumab for an additional 1 month to 5.5 years after progression (median exposure time to trastuzumab of 15.5 months), and 118 patients that were changed to a different chemotherapeutic regimen. Of the 154 patients that continued trastuzumab, 79 (51.3%) showed responses (CR+PR) based on RECIST criteria during first-line trastuzumab treatment. Of the 118 patients that suspended trastuzumab, RECIST responses had been observed in 44 (37.3%). The clinical-pathological characteristics of these patients are shown in Table 1.

Continuation of trastuzumab led to a significant increase in overall survival (p<0.001) in both RECIST responder (Fig. 2a) and non-responder groups (Fig. 2b). Specifically, RECIST responder patients continuing trastuzumab showed a median survival of 27 months (95%CI = 24 to 37) compared to 12 months (95%CI = 7 to 16) for those halting trastuzumab (p<0.001).

Similarly, within RECIST non-responder patients the median survival benefit was 21 (95%CI = 18 to 38) versus 11 months (95%CI = 8 to 15, p<0.001). The strong gain in overall survival in patients non-responsive to first line trastuzumab treatment who continued trastuzumab at disease progression was unexpected and surprising. It is possible that the survival benefit was due to clinical-pathological differences between those patients selected to continue trastuzumab therapy versus those that were switched to a different regimen. The continuation group represented patients of younger age (52 vs. 56 years), smaller size of first metastatic lesion at progression (20 vs. 25 mm), a higher number of resected lymph nodes (18 vs. 15), more taxane (68.8% vs. 39.8%) and gemcitabine (14.3% vs. 0.8%) therapy, fewer visceral-only metastases (8.4% vs.21.2%), more multiple metastases (80.5% vs. 59.3%) and fewer liver-only metastases (6.5% 1.73 to 7.21), with no significant differences in the two estimated HRs (p-value of the likelihoodratio test = 0.690). These results suggest that continued trastuzumab treatment after disease progression has clinically and statistically significant effects in both responder and nonresponder subgroups, even after balancing the differences in the groups who continued or discontinued the treatment. In addition, responders and non-responders did not show significant differences in overall survival (HR = 0.91, 95% CI = 0.46 to 1.67, p = 0.760).

DISCUSSION

Two recent randomized clinical trials have documented the benefit on overall response of continued use of trastuzumab in women with HER2 positive breast cancer who experienced progression during trastuzumab treatment [START_REF] Minckwitz | Trastuzumab beyond progression in human epidermal growth factor receptor 2-positive advanced breast cancer: a german breast group 26/breast international group 03-05 study[END_REF][START_REF] Blackwell | Randomized study of Lapatinib alone or in combination with trastuzumab in women with ErbB2-positive, trastuzumab-refractory metastatic breast cancer[END_REF]. Here we show from a study of clinical practice in Italy that continued trastuzumab treatment after disease progression has clinically and statistically significant effects, remarkably in both RECIST responder and non-responder MBC patients upon first line trastuzumab treatment.

Although we are aware that retrospective studies suffer from limitations, such as the heterogeneity we found in the patients either continuing or halting trastuzumab beyond progression, we reduced bias by building a predicted survival function on matched patients balanced for variables associated with continuation/discontinuation of trastuzumab treatment and for predictors of survival. Finding a significant clinical benefit in both RECIST responder and non-responder patients continuing trastuzumab at progression in Cox proportional hazards analysis of progressed patients, matched using propensity score, strengthened our results. Specifically, in our sample the overall survival approximately doubled in patients with RECIST responses continuing trastuzumab compared to those halting trastuzumab (median overall survival 27 versus 12 months). More interestingly, a similar significant increase (21 versus 11 months) was found also in non-responder patients.

While a number of retrospective analyses have shown a potential advantage for the "continuation policy" (reviewed in [START_REF] Fabi | Do HER-2 positive metastatic breast cancer patients benefit from the use of trastuzumab beyond disease progression? A mono-institutional experience and systematic review of observational studies[END_REF][START_REF] Cancello | Continuing trastuzumab beyond disease progression: outcomes analysis in patients with metastatic breast cancer[END_REF][START_REF] Extra | Efficacy of trastuzumab in routine clinical practice and after progression for metastatic breast cancer patients: the observational Hermine study[END_REF][START_REF] Montemurro | Retrospective evaluation of clinical outcomes in patients with HER2-positive advanced breast cancer progressing on trastuzumab-based therapy in the pre-lapatinib era[END_REF]) in agreement with our results, two retrospective studies indicated no statistically significant benefit from continuing trastuzumab [START_REF] Montemurro | Outcome of patients with HER2-positive advanced breast cancer progressing during trastuzumab-based therapy[END_REF][START_REF] Montemurro | Retrospective evaluation of clinical outcomes in patients with HER2-positive advanced breast cancer progressing on trastuzumab-based therapy in the pre-lapatinib era[END_REF]. In these studies, carried out by the same group, the authors themselves suggest that patients continuing trastuzumab might have been subject to closer clinical evaluation because of the weekly administration schedule of the antibody compared to the control arm of untreated patients, and an earlier detection of disease progression may have partially clouded the benefit of continuation.

To explain the gain in survival beyond progression in RECIST non-responder patients first observed in our series, we speculate that the two major mechanisms of action described for trastuzumab i.e., antibody-dependent cell-mediated cytotoxicity (ADCC) and cytostatic activity through the blockage of HER2 proliferation pathways [START_REF] Spector | Understanding the mechanisms behind trastuzumab therapy for human epidermal growth factor receptor 2-positive breast cancer[END_REF] can explain their effectiveness in responder and non-responder patients, alike. ADCC likely accounts for tumor shrinkage and is thus probably the main mechanism reflected by the RECIST system response evaluation. ADCCmediated tumor destruction requires a functional innate immune system [START_REF] Gennari | Pilot study of the mechanism of action of preoperative trastuzumab in patients with primary operable breast tumors overexpressing HER2[END_REF] presumably compromised in metastatic patients who have received one or more chemotherapy regimens. However, for those patients lacking cellular components that mediate ADCC and that do not account for RECIST responses, trastuzumab therapy might still be beneficial because of its cytostatic activity: inhibition of HER2-mediated proliferation pathways could lead to improved survival even without a decrease in the tumor mass.

With the exclusion of patients with tumors insensitive to HER2 blockage, including those with PTEN or p27 loss [START_REF] Pandolfi | Breast cancer -loss of PTEN predicts resistance to treatment[END_REF][START_REF] Nahta | P27(kip1) down-regulation is associated with trastuzumab resistance in breast cancer cells[END_REF], or expressing a truncated HER2 receptor [START_REF] Scaltriti | Expression of p95HER2, a truncated form of the HER2 receptor, and response to anti-HER2 therapies in breast cancer[END_REF], or displaying enhanced Akt signaling, or signaling through other receptors such as cMet [START_REF] Shattuck | Met receptor contributes to trastuzumab resistance of Her2-overexpressing breast cancer cells[END_REF] and integrins [START_REF] Mocanu | Associations of ErbB2, beta1-integrin and lipid rafts on Herceptin (Trastuzumab) resistant and sensitive tumor cell lines[END_REF], or those producing high amount of growth factors leading to increased activation of HER2 no longer sensitive to antibody blockage [START_REF] Motoyama | The efficacy of ErbB receptor-targeted anticancer therapeutics is influenced by the availability of epidermal growth factor-related peptides[END_REF] continuation of trastuzumab treatment in patients who show no RECIST response may lead to a clinically relevant cytostatic response. Compared with cytotoxic chemotherapy, trastuzumab-based therapy and other similar targeted therapies often result in disease stabilization or delayed disease progression rather than tumor shrinkage. Thus, to predict benefits of trastuzumab continuation at progression calculated in terms of survival beyond progression, indicators other than RECIST should be considered. Promising results obtained in preclinical models have shown that molecular imaging techniques such as PET may be adequate to detect trastuzumab-induced inhibition of HER2+ tumor growth [START_REF] Shah | Imaging biomarkers predict response to anti-HER2 (ErbB2) therapy in preclinical models of breast cancer[END_REF], but clinical explorations in randomized trials are warranted.

In conclusion, our findings suggests that continued trastuzumab treatment after disease progression has clinically and statistically significant effects in both RECIST responder and nonresponder MBC, indicating that continuing trastuzumab beyond disease progression has efficacy also in the latter group. Our results may also suggest that is possible to restart treatment with trastuzumab even in patients relapsed in the adjuvant setting.

Patients

  272 patients with HER2-positive MBC all receiving first line trastuzumab in the metastatic setting during 2000-2001 in 22 Italian oncologic centers and recorded as in progression under trastuzumab treatment during the 2000-2007 observation period were considered for this study. Following progression, 154 patients continued trastuzumab-containing regimens, whereas 118 patients stopped trastuzumab and changed therapy. All primary tumors were scored 3+ by immunohistochemistry (Herceptest) or 2+/FISH-positive prior to first line trastuzumab treatment for MBC.

LEGENDSFigure 1 .

 1 Figure 1. CONSORT Diagram of the Italian retrospective study

Figure 2 .

 2 Figure 2. Overall survival from the first documented progression under trastuzumab therapy for both responder (a) and non-responder groups (b) according to discontinuation (dotted) and continuation (solid line) of trastuzumab treatment.

Figure 3

 3 Figure 3Predicted survivor function from the first documented progression under trastuzumab therapy for matched patients with disease progression, according to discontinuation/continuation of trastuzumab treatment and according to response/no response: discontinuation and response (D/R, solid line), discontinuation and no response (D/NR, dashed), continuation and response (C/R, dotted), and continuation and no response (C/NR, dot-dashed).

n=272HER2-Follow

  Figure&OLFN KHUH WR GRZQORDG )LJXUH )LJXUDUHYLVHGSSW

Table 1 .

 1 Descriptive statistics of the unmatched study sample and of the subgroups with suspension/continuation of Trastuzumab treatment at disease progression. Continuous covariates and counts: median and (interquartile range). Categorical covariates: absolute and (relative) frequency. P-value of the test of independence between covariates and discontinuation/continuation of Trastuzumab treatment at progression. Continuous covariates and counts: Wilcoxon's rank-sum test statistic. Categorical covariates: Fisher's exact test.

		Total sample	Treatment continuation †	Treatment suspension †	p §
	Total	272	154 (56.6)	118 (43.4)
	Patient-related variables
	Age (years)	54 (44 -61)	53 (40 -60)	56 (48 -63)	0.009
	Menopause	184 (69.4)	99 (65.1)	85 (75.2)	0.082
	LVEF (%)	64 (59 -68)	64 (58 -68)	63 (59 -67)	0.846
	Dead patients	131 (48.2)	63 (40.9)	68 (57.6)	0.007
	Tumor-related variables
	ER positive	133 (49.8)	69 (45.1)	64 (56.1)	0.084
	PgR positive	104 (39.2)	52 (34.4)	52 (45.6)	0.076
	Grade (I+II vs. III)	160 (70.8)	85 (68.5)	75 (73.5)	0.464
	Ductal histotype	196 (83.8)	121 (83.4)	75 (84.3)	1.000
	N. of resected lymph nodes	17 (12 -22)	18 (12 -23)	15 (12 -20)	0.037
	N. of positive lymph nodes	4 (1 -10)	4 (1 -9)	3 (1 -10)	0.704
	Relapses and recurrences
	Target/not target metastases	38 (14.0)	13 (8.4)	25 (21.2)	0.001
	Visceral only
	Not Visceral	62 (22.8)	30 (19.5)	32 (27.1)
	Other	172 (63.2)	111 (72.1)	61 (51.7)
	N. of metastatic lesions	2 (1 -3)	2 (2 -3)	2 (1 -3)	0.010
	Multiple metastases	194 (71.3)	124 (80.5)	70 (59.3)	<0.001
	Site of metastasis	30 (11)	10 (6.5)	20 (16.9)	0.022
	Liver
	Bone	12 (4.4)	5 (3.2)	7 (5.9)
	Lung	15 (5.5)	8 (5.2)	7 (5.9)
	Other	215 (79.0)	131 (85.1)	84 (71.2)
	N. of target lesions	1 (1 -2)	2 (1 -2)	1 (1 -2)	0.001
	Size of first metastatic lesion at	22 (15 -35)	20 (11 -30)	25 (20 -40)	0.005
	progression (mm)
	Associated therapies
	Taxane	153 (56.3)	106 (68.8)	47 (39.8)	<0.001
	Gemcitabine	23 (8.5)	22 (14.3)	1 (0.8)	<0.001
	Aromatase inhibitor	34 (12.5)	24 (15.6)	10 (8.5)	0.079
	Other drugs	23 (8.5)	15 (9.7)	8 (6.8)	0.384
	RECIST response to first-line trastuzumab-based therapy	123 (45.2)	79 (51.3)	44 (37.3)	0.021
	§

ACKNOWLEDGEMENTS

This work was supported by AIRC (Associazione Italiana per la Ricerca sul Cancro) and FIRB/MUR (Ministero Università e Ricerca) and Roche spa.

vs. 16.9%), a higher number of target lesions (2 vs. 1), and a higher response to trastuzumabbased therapy (51.3% vs. 37.3%) in the metastatic disease (see also Supplementary Table 1).

To reduce bias due to patient selection, we used the method of propensity score matching estimated by a logistic regression model using continuation/discontinuation of trastuzumab as the outcome variable and the following covariates as predictors: taxane treatment, number of resected lymph nodes, axillary dissection, tumor size, response to trastuzumab-based therapy, multiple metastases, gemcitabine, time to progression, target and non-target metastases, age, site of metastases, aromatase inhibitor, number of positive lymph nodes and number of target lesions. These variables were selected using the graph of the estimated variable importances in the Supplementary Figure 1. The area under the ROC curve of the logistic model was 0.80 (95% CI = 0.75 to 0.86). The estimate values of the model's parameters are given in Supplementary Table 2.

The most predictive variables of survival were: response to trastuzumab therapy, size of first metastatic lesion at progression, menopause, brain metastases, site of metastases, other drugs associated to trastuzumab therapy, time to progression, aromatase inhibitor associated to trastuzumab therapy, and HER2 score (see the estimated variable importance in the Supplementary Figure 2). Propensity score 1:1 matching with replacement generated 238 matched treated and control patients. The two groups reached a "good" balance (standardized mean differences below 10) for all the variables associated with continuation/discontinuation of trastuzumab treatment and for predictors of survival (see Supplementary Table 1). Limiting the analysis to the matched patients (Fig. 3) still showed that discontinuation of trastuzumab at metastatic disease progression was a risk factor that significantly reduced overall survival in both responder (HR = 2.23; 95% CI = 1.03 to 4.82) and non-responder groups (HR = 3.53, 95% CI =