M Kieffer

E Walter

Guaranteed estimation of the parameters of nonlinear continuous-time models: contributions of interval analysis †

Keywords: Bounded-error estimation, guaranteed estimation, interval analysis, parameter bounding, parameter optimization, sensitivity functions

This paper is about guaranteed parameter estimation in two contexts, namely bounded-error and optimal estimation. In bounded-error estimation, one looks for the set of all parameter vectors that are consistent with some prior bounds on the errors deemed acceptable between the model behavior and that of the system. In optimal estimation, one looks for the set of all parameter vectors that minimize some cost function quantifying the discrepancy between the behaviors of the system and its model. In both cases, guaranteed means that proven statements are made about the set of interest. The situation is made much more difficult when the model output is assumed to depend nonlinearly in the parameters to be estimated and when dealing with continuous-time models, as here. Important tools based on interval analysis that contribute to allowing guaranteed estimation in these challenging conditions are presented. Some are absolutely classical in the context of interval analysis but not so well known in the community of parameter estimation at large. Others have been developed recently and were mainly presented in conferences. Some, such as the use of sensitivity functions to reduce more quickly the size of outer approximations of the sets of interest, are new. Challenges for future research on the context of guaranteed nonlinear estimation are mentioned.

INTRODUCTION

Parameter estimation is classically conducted by minimizing some cost function that may be purely heuristic or deduced from information (or hypotheses) about the noise corrupting the data and possibly about the prior distribution for the parameter vector.

Bounded-error parameter estimation, by contrast, does not use a cost function and aims instead at characterizing the set of all values of the parameter vector that are consistent with bounds on the errors one is prepared to accept between the behavior of the model and that of the system to be modeled. It can be traced back to the seminal work of F.C. Schweppe [START_REF] Schweppe | Recursive state estimation: unknown but bounded errors and system inputs[END_REF][START_REF] Schweppe | Uncertain Dynamic Systems[END_REF] on state estimation for linear models. His idea was recursively to compute ellipsoids guaranteed to contain the actual state of the system provided that some bounds on the state perturbations and measurement noise were satisfied. This type of parameter or state bounding based on a description of uncertainty by sets is often referred to as guaranteed estimation, see, e.g., [START_REF] Chernousko | Optimal guaranteed estimates of indeterminacies with the aid of ellipsoids[END_REF] and the books [START_REF] Chernousko | State Estimation for Dynamic Systems[END_REF] and [START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF]. We shall see in this paper that guaranteed can be given a broader meaning in our context. Of course, other types of containers than ellipsoids could and have been used, such as boxes [START_REF] Milanese | Estimation theory and uncertainty intervals evaluation in presence of unknown but bounded errors: Linear families of models and estimators[END_REF], parallelotopes [START_REF] Chisci | Recursive state bounding by parallelotopes[END_REF], zonotopes [START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF] or other limited-complexity polytopes. Matlab toolboxes implementing ellipsoidal or polytopic calculus are readily available [START_REF] Kurzhanski | Ellipsoidal toolbox[END_REF], [START_REF] Veres | Users manual -reference of the geometric bounding toolbox -version 7.3[END_REF].

Parameter or state bounding when the model output depends nonlinearly on the quantities to be estimated is of course much more difficult, and nonlinear ellipsoidal calculus remains a largely open subject. A generalization of the ellipsoidal approach to nonlinear state estimation is briefly considered in [START_REF] Chernousko | State Estimation for Dynamic Systems[END_REF]. It is based on the use of auxiliary (time-varying) linear state equations on which ellipsoidal calculus can be applied. Outer (or inner) approximations of the set of feasible state vectors at any given instant of time can then be computed under the form of ellipsoids using the same methodology as for linear models, provided suitable inclusion relations between the nonlinear and auxiliary models can be established. Along a similar line, one may also try to consider a nonlinear model as a linear one with uncertain regressors, where the uncertainty is described by intervals, and use techniques such as described in [START_REF] Nazin | Interval parameter estimation under model uncertainty[END_REF]. In both cases, the nonlinear problem is transformed into a linear uncertain one to which linear techniques can be applied. Linear models with time-varying coefficients could also be considered, as in [START_REF] Kurzhanski | Ellipsoidal techniques for reachability analysis[END_REF].

Instead of attempting linearization, one may use the tools provided by interval analysis directly to compute an approximate but guaranteed characterisation of sets in the nonlinear case. In the context of bounded-error parameter estimation, this approach made it possible to derive algorithms that could be used to bracket the set of possible values for the parameter vector between computable unions of boxes in parameter space [START_REF] Jaulin | Set inversion via interval analysis for nonlinear bounded-error estimation[END_REF]. At the beginning, these algorithms were limited to models for which a closed-form expression of the model output as a function of the parameters was available, but the concept of guaranteed integration [START_REF] Berz | Verified integration of ODEs and flows using differential algebraic methods on high-order Taylor models[END_REF][START_REF] Lohner | Enclosing the solutions of ordinary initial and boundary value problems[END_REF][START_REF] Nedialkov | Methods for initial value problems for ordinary differential equations[END_REF] made it applicable to models for which no such expression exists [START_REF] Kieffer | Guaranteed parameter estimation for cooperative systems[END_REF][START_REF] Raissi | Set membership state and parameter estimation for systems described by nonlinear differential equations[END_REF][START_REF] Walter | Guaranteed nonlinear parameter estimation in knowledge-based models[END_REF]. This paper will not describe guaranteed integration in any detail, but it may be useful to mention that two approaches at least could be considered. The use of high-order Taylor models is considered in [START_REF] Berz | Verified integration of ODEs and flows using differential algebraic methods on high-order Taylor models[END_REF][START_REF] Lohner | Enclosing the solutions of ordinary initial and boundary value problems[END_REF][START_REF] Nedialkov | Methods for initial value problems for ordinary differential equations[END_REF], while W. Kühn uses a set-valued dynamical system approach based on the description of sets by zonotopes, where again interval analysis plays an important role [START_REF] Kühn | Rigorously computed orbits of dynamical systems without the wrapping effect[END_REF].

The success of interval analysis as a tool to compute guaranteed set estimates in the context of bounded-error estimation seems to have somewhat obscured the fact that the very same tools can also be used to obtain guaranteed estimates in the more traditional approach where a cost function is to be optimized.

Besides presenting our most recent results on the use of interval analysis and guaranteed integration for bounded-error nonlinear parameter estimation in the context of continuoustime models, we would therefore like to show how these tools can also contribute to the more mainstream parameter estimation through optimization. We shall also explain why this is actually more difficult, and point out challenges for future research.

NOTATION AND PROBLEM STATEMENT

We assume that the model consists of a set of ordinary differential equations in state-space form, involving some vector p of parameters to be estimated

x ′ (t) = f (x (t) , p, t) , y m (p, t) = h (x (t) , p) . (1)
x ′ (t) is the derivative of x with respect to time and the dependency of f in t may result from the influence of a known, time-varying vector u of external inputs. For the sake of simplicity, the initial state x (t 0 ) is assumed here to be known, but it could also be parametrized. With such a model, y m (p, t), as computed from (1), is almost always nonlinear in the parameters p, even when it is linear in some input. Let y(t i ) be the vector of the measurements collected on the system at the instants of time t i , i = 1, . . . , n t , and let y be the vector obtained by concatenating all y(t i )s y = y T (t 1 ), . . . , y T (t nt )

T .

(

) 2
Two types of estimation problems will be considered, namely those based on the minimization of some cost function and those involving the characterization of a set of parameter vectors that satisfy error bounds. We shall now state these problems more precisely and explain what we mean by guaranteed in both cases.

Guaranteed estimation via optimization

The most classical statistical estimation methods involve the optimization of a cost function c (p), which may be deduced from probabilistic assumptions on the noise and parameters, see, e.g., [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF][START_REF] Walter | Identification of Parametric Models from Experimental Data[END_REF]. For example, if one assumes that there exists a true value p * of the parameter vector and that the measurements satisfy

y = y m (p * ) + e, (3)
with y m (p) consisting of subvectors y m (p, t i ) and e a zero-mean white Gaussian measurement noise, then the maximum-likelihood estimate of p is the argument of the minimum of the cost function

c (p) = (y -y m (p)) T (y -y m (p)) . (4)
When y m (p) is linear in p, it is easy, at least in principle, to check whether the experimental conditions are such that the global minimizer of ( 4) is unique, and to compute it when this is the case. When y m (p) is nonlinear in p, as is most often the case with models such as (1), optimization becomes much more difficult. The usual local methods (such as quasi-Newton or conjugate gradient algorithms) converge at best to a local minimizer of the cost function.

Global optimization based on random search (for instance simulated annealing or genetic algorithms) cannot be guaranteed to locate global minimizers in finite time, so that, most often, no guarantee at all is provided about the estimate. Often, we do not even know how many global minimizers there is. In this context, guaranteed estimation aims at deriving proven statements about the global minimum of the cost function and the associated set of global minimizers (which is also the set of optimal estimates).

Guaranteed estimation via parameter bounding

An alternative parameter estimation approach is parameter bounding , see [START_REF]Bounding Approaches to System Identification[END_REF][START_REF] Norton | Special Issue on Bounded-Error Estimation: Issue 1[END_REF][START_REF] Norton | Special Issue on Bounded-Error Estimation: Issue 2[END_REF][START_REF] Walter | Special Issue on Parameter Identification with Error Bounds[END_REF] and the references therein. In parameter bounding, one looks for the set of all parameter vectors that are consistent with the experimental data, the structure of the model, and the bounds on the errors one is prepared to accept. With each experimental data vector y (t i ), we assume that we can associate some known vector of intervals (or box) [e i , e i ] of acceptable errors. The parameter vector p, assumed to belong to some known prior box [p] 0 , is deemed acceptable if e i y (t i )y m (t i , p) e i for all i = 1, . . . , n t .

(

) 5
Parameter estimation then amounts to characterizing the set

S = {p ∈ [p] 0 | p is acceptable} = {p ∈ [p] 0 | y m (p) ∈ [y]} , (6) with
In this context, guaranteed estimation aims at deriving proven statements about S, such as S is void or S is included in this computed set. When y m (p) is linear in p, a number of methods have been proposed to compute sets guaranteed to contain S, which is then a convex polytope. When the model is nonlinear in p, as assumed here, S may consist of disconnected and nonconvex subsets, which makes the situation much more difficult to handle.

Aims of the paper

This paper is aimed at presenting important tools based on interval analysis and guaranteed integration that contribute to allowing guaranteed estimation in these challenging conditions. Some are absolutely classical in the community of interval analysts but not so well known in the control and statistics communities. Others have been developed recently and were mainly presented in conferences. Some, such as the use of sensitivity functions to reduce more quickly the size of outer approximations of the sets of interest, are new. A seemingly simple but nevertheless challenging example of nonlinear parameter estimation problem will be treated with the optimal and bounding approaches, and will be used to demonstrate the interest of the new approach for set reduction based on the use of sensitivity functions.

TOOLS FOR GUARANTEED ESTIMATION

In this section, we briefly recall the three main tools to be used to make estimation guaranteed in either of the senses defined above. Interval analysis (IA), which is at the core of the other two, is presented first. We then proceed to guaranteed global optimization, and guaranteed set characterization.

Specific tools, such as guaranteed integration and sensitivity functions, which are useful to consider the estimation of systems described by continuous-time state-space models, will be introduced in Section 4.

Interval analysis

IA [START_REF] Jaulin | Applied Interval Analysis[END_REF][START_REF] Moore | Methods and Applications of Interval Analysis[END_REF][START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF] provides tools to compute with vectors of intervals (more concisely, boxes)

[x] = [x, x] ⊂ R n . For any function f : D ⊂ R n -→ R m
defined by combining arithmetical operators and elementary functions, IA makes it possible to build inclusion functions

[f ] satisfying ∀ [x] ⊂ D, f ([x]) ⊂ [f ] ([x]) , (9)
where

[f ] ([x]
) is a box. Natural inclusion functions, for example, are obtained by replacing all real variables, operators and elementary functions by their interval counterparts. It thus becomes possible to enclose the set f ([x]) , which most often cannot be computed, in a computable box

[f ] ([x]). If f is differentiable over [x], the mean-value theorem implies that ∀x ∈ [x] , ∃z ∈ [x] such that f (x) = f (m) + J (z) • (x -m) ,
where J is the Jacobian matrix of f and m = mid([x]). Thus,

∀x ∈ [x] , f (x) ∈ f (m) + [J] ([x]) • (x -m) ,
where [J] is an inclusion function for J. So

f ([x]) ⊂ f (m) + [J] ([x]) • ([x] -m) .
Therefore, the interval function

[f c ] ([x]) = f (m) + [J] ([x]) • ([x] -m) (10)
is another inclusion function for f , called its centered inclusion function. This function becomes especially interesting when the width of [x] is small, because the overestimation resulting from the interval evaluation of [J]([x]) is reduced by the product with [x]m, which is a small box centered on zero. The centered inclusion function can be significantly improved at the cost of a slightly more complicated formulation [START_REF] Hansen | On solving systems of equations using interval arithmetic[END_REF].

Guaranteed global optimization

Let P ⊂ R dim p be a (possibly very large but compact) domain of interest in parameter space.

For simplicity, here P will be a box [p] 0 . We wish to characterize the set S of all the global minimizers of the cost function that are inside

P S = arg min p∈P c (p) . (11)
These minimizers are not expected to lie on the boundary of P, so this is unconstrained minimization. We now present the principles of a simplified version of Hansen's algorithm, see [START_REF] Hansen | Global Optimization Using Interval Analysis[END_REF][START_REF] Jaulin | Applied Interval Analysis[END_REF] for more details.

Let L be a list of boxes, the union of which contains all global minimizers belonging to P. Initially, L contains only [p] 0 . The basic structure of the algorithm is as follows, where ε > 0 is a tolerance parameter to be tuned by the user.

[p] is such that 0 / ∈ [g]([p]), (12)
then no p in [p] can be an unconstrained minimizer so [p] can be eliminated. If a contractor is available, see Section 3.4, Step 4 uses it to reduce the size of [p] without bisecting it, which makes Step 4 an essential ingredient in the struggle against the curse of dimensionality. A number of potentially useful contractors are presented in [START_REF] Jaulin | Applied Interval Analysis[END_REF].

Since the tolerance parameter ε is strictly positive, the algorithm stops after a finite number of steps and the set S of all global minimizers of c in P is included in the union of the boxes in L, each of which has now a width smaller than ε.

Set Inversion Via Interval Analysis

Let f be a possibly nonlinear function from R np to R ny and Y be a subpaving of R ny , i.e., a union of nonoverlapping boxes. Set inversion is the characterization of the reciprocal image of Y

S = {p ∈ R np | f (p) ∈ Y} = f -1 (Y). ( 13
)
Using an inclusion function [f ] for f , two subpavings S and S such that S ⊂ S ⊂ S can be obtained with the algorithm Sivia, for Set Inverter Via Interval Analysis [START_REF] Jaulin | Set inversion via interval analysis for nonlinear bounded-error estimation[END_REF]. A (possibly very large) search box [p] 0 to which S is guaranteed to belong must be provided by the user. It will be used as an initial outer approximation to grow inner and outer approximations of S under the form of subpavings by successive bisections and selections. Consider a box [p] obtained by bisection of some box belonging to the current outer subpaving. Several cases may be encountered. The resulting recursive algorithm is given in Table I, where the subpavings S and S should be initialized as empty before the first call.

Algorithm Sivia(in: f , Y, [p] , ε; inout: S, S) 1 [p] := C ([p]); 2 if [p] = ∅ or [f ] ([p]) ∩ Y = ∅ return; 3 if [f ] ([p]) ⊂ Y then 3 {S := S ∪ [p] ; S := S ∪ [p ′ ]; return;}; 4 if w ([p]) < ε then {S := S ∪ [p]; return;}; 5 Sivia(f , Y, L [p] , ε, S, S); Sivia(f , Y, R [p] , ε, S, S). Table I. Sivia
All the boxes of S that are not in S form an uncertainty layer ∆S. All of them have a width smaller than ε, so the accuracy of the description of S may be increased by decreasing ε, at the cost of more computations.

Contractors

Assume that the parameters to be estimated must satisfy the vector inequality constraint k (p) 0, to be understood componentwise. Assume moreover that p is known to belong to

[p]. A contractor C k for p is then an algorithm to compute a box C k ([p]) such that C k ([p]) ⊂ [p] , {p ∈ [p] |k (p) 0} ⊂ C k ([p]) . (14)
The first relation in [START_REF] Kieffer | Interval analysis for guaranteed nonlinear parameter estimation[END_REF] ensures that p is contracted, while the second guarantees that no value of p satisfying the constraints is lost. Contractors can be similarly defined in the case of equality constraints.

When a closed-form expression of the cost function is available, this constraint may be used to build a contractor for p via interval constraint propagation, see [START_REF] Jaulin | Applied Interval Analysis[END_REF]. When this is not the case, as usual when considering state-space models such as those described by [START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF], one may use a centered inclusion function of the cost function over

[p] = [p 1 ] × • • • × p np to get from (15) c (m) + np i=1 ([p i ] -m i ) ∂c ∂p i ([p]) ⊂ [0, c] , (16)
since the cost function is positive. Thus each [p i ], i = 1, . . . , n p , has to satisfy

[p i ] ⊂   [0, c] -c (m) - j =i ([p j ] -m j ) ∂c ∂p j ([p])   / ∂c ∂p i ([p]) + m i , (17)
leading to a contracted box [p] new = C c ([p]
), the components of which are defined as

[p i ] new = [p i ] ∩     [0, c] -c (m) - j =i ([p j ] -m j ) ∂c ∂p j ([p])   / ∂c ∂p i ([p]) + m i   (18)
for i = 1, . . . , n p . The main advantage of this contractor is that it only requires the computation of an inclusion function for the gradient of the cost function, which is easily obtained using sensitivity functions, see Section 4.

Using [START_REF] Kurzhanski | Ellipsoidal Calculus for Estimation and Control[END_REF] new = C k ([p]), the components of which are defined as

(m) + np j=1 ([p j ] -m j ) ∂y m k ∂pj ([p]) -y k 0, -y m k (m) - np j=1 ([p j ] -m j )
[p i ] new = [p i ] ∩     y k , y k -y m k (m) - j =i ([p j ] -m j ) ∂y m k ∂p j ([p])   / ∂y m k ∂p i ([p]) + m i   , (20)
for i = 1, . . . , n p . A similar procedure may be applied for all components of y m (p). Again, only an inclusion function for the Jacobian matrix of y m (p) is needed to apply the contractor C k . For state-space models such as those described by [START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF] , this inclusion function may be obtained using sensitivity functions, see Section 4.2.

With standard interval arithmetic, contraction may only be obtained in [START_REF] Kurzhanski | Ellipsoidal toolbox[END_REF]

when 0 / ∈ [∂y m k /∂p i ] ([p]
). The conditions leading to 0 ∈ [∂y m k /∂p i ] ([p]) depend mainly on the structure of the model and the overestimation introduced by the inclusion function for ∂y m k /∂p i . They are not necessarily linked to whether [p] is in the solution set S. Extended intervals could readily be used to improve the efficiency of [START_REF] Kurzhanski | Ellipsoidal toolbox[END_REF]. With the Sivia algorithm, boxes tend to accumulate on the boundary of S [START_REF] Jaulin | Set inversion via interval analysis for nonlinear bounded-error estimation[END_REF]. For boxes that are very small and across the boundary of S, no significant contraction may be obtained. However, for the first iterations of Sivia, very large boxes may be efficiently eliminated or their size may be strongly reduced using C k .

Many other contractors may be considered. For example the Newton contractor, which can be seen as an interval counterpart to the Newton algorithm, can be used to reduce the search domain in global optimization. Assuming that the cost function is sufficiently differentiable, it exploits the constraint g (p) = 0, that must be satisfied by any p to qualify as an unconstrained global minimizer. This, however, requires an inclusion function for the Hessian of the cost function, which significantly complicates implementation and execution. Other interesting contractors include the Krawczyk contractor and contractors based on constraint propagation. For more details, see [START_REF] Jaulin | Applied Interval Analysis[END_REF].

DEALING WITH CONTINUOUS-TIME MODELS

Apart from the evaluation of an inclusion function of the cost function c, the simplified version of Hansen's algorithm presented in Section 3.2 requires the evaluation of an inclusion function of its gradient g. When a closed-form expression of the model output is available, a formal expression of the gradient is easily obtained by formal differentiation and the standard methods for obtaining inclusion functions apply. In general, however, it is not possible to obtain closedform solutions for models such as those described by [START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF], which makes the obtention of inclusion functions for c and g particularly difficult.

To get an inclusion function for c, the naive approach would be to use one of the guaranteed ODE solvers provided by IA, e.g., AWA [START_REF] Lohner | Computation of guaranteed enclosures for the solutions of ordinary initial and boundary value-problem[END_REF], COSY [START_REF] Hoefkens | Efficient high-order methods for ODEs and DAEs[END_REF] or VNODE [START_REF] Nedialkov | Methods for initial value problems for ordinary differential equations[END_REF]. Unfortunately, these solvers prove unable to provide accurate enclosures for the solutions when there are uncertain parameters, as here. The main idea to bypass this difficulty is to bound the solutions of uncertain dynamical systems using deterministic dynamical systems. Efficient guaranteed numerical integrators can then be used to compute the corresponding bounding solutions. Preliminary results were presented in [START_REF] Kieffer | Interval analysis for guaranteed nonlinear parameter and state estimation[END_REF][START_REF] Walter | Interval analysis for guaranteed nonlinear estimation[END_REF] in the context of cooperative dynamical models, i.e., models such as [START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF] for which the off-diagonal terms of the Jacobian matrix of f are positive. These results were inspired by the interval observer proposed in [START_REF] Gouzé | Interval observers for uncertain biological systems[END_REF]. Müller's theorems [START_REF] Müller | Über das Fundamentaltheorem in der Theorie der gewöhnlichen Differentialgleichungen[END_REF], which have recently been used in the context of guaranteed simulation [START_REF] Gennat | Simulation of uncertain systems with guaranteed bounds[END_REF], make it possible to bound the solutions of more general dynamical models.

Müller's theorem

In order to apply one of the theorems in [START_REF] Müller | Über das Fundamentaltheorem in der Theorie der gewöhnlichen Differentialgleichungen[END_REF] to the bounding of the solutions of uncertain dynamical models such as [START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF] in the presence of uncertain parameters p ∈ [p 0 , p 0 ] and uncertain initial conditions x (t 0 ) ∈ [x 0 , x 0 ], we reformulate it as follows.

Theorem 1. Assume that f (x, p, t) in (1) is continuous on

T :    ω (t) x Ω (t) p 0 p p 0 t 0 t t 1
where ω i (t) and Ω i (t) , i = 1 . . . n x , are continuous on [t 0 , t 1 ] and such that 1. ω (t 0 ) = x (t 0 ) and Ω (t 0 ) = x (t 0 ), 2. the left derivatives D -ω i (t) and D -Ω i (t) and right derivatives D + ω i (t) and

D + Ω i (t) of ω i (t)
and Ω i (t) satisfy, for i = 1 . . . n x and t ∈ [t 0 , t 1 ] ,

D ± ω i (t) min T i (t) f i (x, p, t) and D ± Ω i (t) max Ti(t) f i (x, p, t) ,
where T i (t) is the subsets of T defined by

T i (τ ) :        x i = ω i (τ ) , ω j (τ ) x j Ω j (τ ) , j = i, p 0 p p 0 , t = τ,
and T i (t) is the subset of T defined by

T i (τ ) :        x i = Ω i (τ ) ω j (τ ) x j Ω j (τ ) , j = i, p 0 p p 0 , t = τ.
Then, for any given x (t 0 ) ∈ [x 0 , x 0 ] and p ∈ [p 0 , p 0 ], a solution to (1) exists such that

ω (t) x(t) Ω (t) ∀t ∈ [t 0 , t N ],
and this solution is equal to x (t 0 ) at t = t 0 . Moreover, if for any p ∈ [p 0 , p 0 ] and t ∈ [t 0 , t 1 ] , f (x, p, t) is Lipschitz with respect to x, then for any given x (t 0 ) ∈ [x 0 , x 0 ] and p ∈ [p 0 , p 0 ] this solution is unique. ♦

This theorem allows the evaluation of lower and upper bounds for the solution of (1) provided that two functions ω (t) and Ω (t) are available. The interval function [Φ] (t) = [ω (t) , Ω (t)] can then be seen as an inclusion function for all solutions x (t) of the state equation in [START_REF] Alamo | Guaranteed state estimation by zonotopes[END_REF].

The construction of ω (t) and Ω (t) is usually easy on a case-by-case basis, as illustrated in Section 5.1.

An inclusion functions for g may similarly be obtained using the sensitivity functions of the output with respect to the parameters, see Section 4.2.

Sensitivity functions

Denote the first-order sensitivity of x j with respect to p k by s jk

s jk (p, t) = ∂x j ∂p k (p, t) . (21)
Assume, for the sake of simplicity, that the model output is linear in the state and given by

h (x (t) , p) = Mx (t) , (22)
where M is a known matrix. The Jacobian matrix of h (x (t) , p) is then given by

J h (p, t) = M ∂x (p, t) ∂p , (23)
with ∂x (p, t) ∂p = (s jk (p, t)) , j = 1, . . . , dim x, k = 1, . . . , dim p. (24)
When considering all measurements between t 1 and t nt , as in ( 8), one gets

∂y m (p) ∂p =     M ∂x(p,t1) ∂p . . . M ∂x(p,tn t ) ∂p     . (25)
The gradient of c may then be obtained as

dc dp (p) = -2 (y -y m (p)) T dy m (p) dp ,
where dym(p) dp may be evaluated using first-order sensitivity functions. To compute s jk , differentiate the jth row of (1) to get

s ′ jk = ∂f j (x, p) ∂x j s jk + ∂f j (x, p) ∂p k . (26)
Since x (t 0 ) is assumed to be known, the initial conditions are

s jk (t 0 ) = ∂x (t 0 ) ∂p k = 0.
The sensitivity function may then be obtained by considering the extended state-space model consisting of the dynamical part of (1), and of all differential equations (26) satisfied by the sensitivity functions. Müller's theorem will then turn out to be especially useful, as the extended state-space model is seldom cooperative, even if this is the case of the initial state space model, see Section 5.

EXAMPLES

Consider the seemingly simple compartmental model [START_REF] Kieffer | Interval analysis for guaranteed nonlinear parameter estimation[END_REF] described by Figure 1. Although it is linear from an input-output point of view, its output is highly nonlinear in its parameters. Its state equation is readily obtained from conservation law as

x ′ = f (x, p, u) , (27)
where p = (k 21 , k 12 , k 01 ) T and

f (x, p, u) = -(p 1 + p 3 ) x 1 + p 2 x 2 + u p1x 1 -p 2 x 2 . (28)
The quantity of material in Compartment 2 is assumed to be observed, so

y m (t i , p) = x 2 (t i , p) , i = 1, ..., n t .
Assume that there is no input (u ≡ 0) and that the initial condition is known to be x 0 = (1, 0) T .

Differentiating [START_REF] Nazin | Interval parameter estimation under model uncertainty[END_REF] with respect to each of the parameters in turn, one gets the differential equations

               s ′ 11 = -(p 1 + p 3 ) s 11 + p 2 s 21 -x 1 , s ′ 21 = p 1 s 11 -p 2 s 21 + x 1 , s ′ 12 = -(p 1 + p 3 ) s 12 + p 2 s 22 + x 2 , s ′ 22 = p 1 s 12 -p 2 s 22 -x 2 , s ′ 13 = -(p 1 + p 3 ) s 13 + p 2 s 23 -x 1 , s ′ 23 = p 1 s 13 -p 2 s 23 . (29)
Since x 0 does not depend on p, the initial conditions for all of these sensitivity equations are zero. Note that the 8th-order extended state equation to be solved for computing x and the s jk s is not cooperative, which makes the Müller theorem recalled in Section 4.1 particularly useful. Since all p k s and x j s are positive, it is easy, using Müller's theorem, to build bounding dynamical systems to evaluate guaranteed enclosures of the model output y m (p), of the cost c and of its gradient g with respect to p. One may, for instance, compute guaranteed enclosures of x, s 11 and s 21 via the simulation of the two 6th-order deterministic ODEs

                 x ′ 1 = -(p 1 + p 3 )x 1 + p 2 x 2 , x ′ 2 = p 1 x 1 -p 2 x 2 , x ′ 1 = -(p 1 + p 3 )x 1 + p 2 x 2 , x ′ 2 = p 1 x 1 -p 2 x 2 , s ′ 11 = -(p 1 + p 3 )s 11 + p 2 s 21 -x 1 , s ′ 21 = p 1 s 11 -p 2 s 21 + x 1 (30)
and

                 x ′ 1 = -(p 1 + p 3 )x 1 + p 2 x 2 , x ′ 2 = p 1 x 1 -p 2 x 2 , x ′ 1 = -(p 1 + p 3 )x 1 + p 2 x 2 , x ′ 2 = p 1 x 1 -p 2 x 2 , s ′ 11 = -(p 1 + p 3 )s 11 + p 2 s 21 -x 1 , s ′ 21 = p 1 s 11 -p 2 s 21 + x 1 . (31)
The same procedure can be applied for the sensitivities with respect to p 2 and p 3 , and one could, at least in principle, build a single 16th-order deterministic ODE, the solution of which would provide all required bounding functions. This is not to be recommended, however, as the complexity of guaranteed ODE solvers depends heavily on the state dimension, so one is better off dealing with six 6th-order deterministic ODEs.

To generate artificial data, a two-compartment model with a "true" value of the parameter vector p * = (0.6, 0.15, 0.35)

T has been simulated. The data were then obtained by rounding the value of x 2 (t i ) to the nearest two-digit number for t i = i∆t, with ∆t = 1 s and i = 1, . . . , 15. The initial search domain is [p] 0 = [0.01, 1] ×3 .

Global Optimization

The cost function is written as

c (p) = nt i=1 (x 2 (p, t i ) -y (t i )) 2 .
For a precision parameter ε = 0.001, the deterministic global optimization algorithm described in Section 3.2 with the contractor C c presented in Section 3.4 provides a list of boxes whose projections onto the (p 1 , p 2 ) plane and (p 2 , p 3 ) plane are shown in Figure 2. Only boxes for which it was not possible to prove that they do not contain any global minimizer are represented, so this is an outer approximation. This result has been obtained in 3 h on a Pentium IV at 2 GHz. The symmetry of the projection on the (p 2 , p 3 )-plane of the solution boxes suggests an identifiability problem that could be confirmed by a theoretical study [START_REF] Walter | Identification of Parametric Models from Experimental Data[END_REF]: the values of the parameters p 2 and p 3 may be exchanged without modifying the model output. This identifiability study was not used to get the results presented here.

Bounded-error estimation

The Sivia algorithm described in Section 3.3 is now employed. To be tolerated, the error must belong to [e] = [-0.005, 0.005] for all measurements. Three versions have been considered. In Using the same data as in the optimization case, we obtain the results summarized in Tables II andIII. The guaranteed numerical integration toolbox VNODE [START_REF] Nedialkov | Methods for initial value problems for ordinary differential equations[END_REF] has been used to obtain guaranteed solutions to [START_REF] Neumaier | Interval Methods for Systems of Equations[END_REF], [START_REF] Norton | Special Issue on Bounded-Error Estimation: Issue 1[END_REF] , and to similar systems of ODEs used for evaluating the sensitivities with respect to p 2 and p For a given value of ε, NIF is the most efficient in terms of computational complexity, see Table II, since only cooperativity is used, and no sensitivity function is involved. The price to be paid is, by far, the worst performance in terms of accuracy, measured by the volume of the outer approximation to the solution set, see Table III. When considering a centered form, since the accuracy of the inclusion function is much better, much less undetermined boxed have to be considered. The volume of the solution set is then much smaller for a given ε. Nevertheless, the centered form is much more costly to obtain, since additional ODEs have to be solved to get the sensitivity functions. Using the contractor in conjunction with the centered inclusion function allows to reduce computing time significantly, at least for large values of ε. The volume of the solution set is also significantly reduced, when compared to the other two approaches. For larger values of ε, CF-CP behaves similarly to CF in terms of computing time, but the quality of the solution is still better for CF-CP. The advantage brought by the contractors reduces when ε decreases, since many small boxes accumulates at the boundary of the solution set. The boxes within the solution set and those at the boundary cannot be reduced by contraction. This explains why the number of simulations become similar. precision parameter ε has to be very small (less than 0.001) to make it possible to see that the actual solution set S consists of two subsets, suggesting an identifiability problem. This conclusion is already reached with ε = 0.005 for CF and with ε = 0.01 for CF-CP.

The results of Tables II and III suggest combining CF-CP and NIF within a single search algorithm. For each box [p], CF-CP should only be used when the contractor C k applied on [p] may produce a contraction significant enough to compensate for the additional computing time required to obtain the sensitivity functions needed by C k . If obtaining a NIF takes α% of the time required to obtain a CF-CP, one may stop using CF-CP as soon as a the ratio of the volume of the contracted box to that of the initial box is larger than α%. Many other strategies may be considered.

CONCLUSIONS AND PERSPECTIVES

Interval analysis is an attractive option for anyone who feels disturbed by the fact that the usual techniques for nonlinear parameter estimation are unable to deliver any guarantee about the estimates that they produce. In this paper, we have presented some already available tools that allow global and guaranteed results to be obtained in nontrivial cases involving continuoustime state-space models. Guaranteed simulation, in combination with Müller's theorems, allow one to deal with ODEs for which no closed-form solution is available, in a more general setting than the notion of cooperativity previously used. It should be clear, however, that guaranteed estimation is far from being applicable to all nonlinear estimation problems of interest, In the general case, the characterization of the set of all optimal or acceptable parameter vectors is an NP-complete problem, and the struggle against the curse of dimensionality is absolutely crucial. The use of contractors, which reduces the size of boxes of interest without bisecting them or losing any solution, is an important ally in this endeavor. We have shown how sensitivity functions, combined with guaranteed integration and one of Müller's theorems, can be used to build such contractors.

Guaranteed parameter optimization turns out to be more difficult than guaranteed parameter bounding. The presence of a sum of terms depending on the uncertain parameters in the cost function entails overestimation in the evaluation of the associated inclusion function. This is because all interval quantities are treated as independent, even if they actually correspond to the same box in parameter space. More efficient contractors than those employed here could be considered to address this difficulty. A first option, which we plan to investigate, is to build contractors for each term in the sum, taking into account the fact that each of them must comply with the constraint [START_REF] Kieffer | Guaranteed parameter estimation for cooperative systems[END_REF]. Convexity tests based on second order sensitivity functions could also be considered.

Moreover, significant improvement in computing times can be expected by an appropriate management of the switching between NIF and CP-CF, which may differ from box to box and from time to time.

1 . 2 .

 12 While L contains at least one box with width greater than ε, do { (a) Pop a box [p] out of L, (b) Update the upper bound ub( c) of the value of the global minimum c over [p] 0 . (c) Try to eliminate [p] via the midpoint and monotonicity tests. (d) Should this fails, try to eliminate or at least reduce [p] using a contractor. (e) If [p] resists elimination and the width of [p] is greater than ε, then bisect [p] and push the two resulting boxes into L. Otherwise push [p] back as is into L. } At last eliminate all boxes of L that no longer pass the midpoint test. A good strategy for selecting the box to be considered at Step 1 is to take one of the boxes with the largest potential, i.e., a box [p] such that the lower bound of [c]([p]) is minimal, with [c] an inclusion function for the cost function c. To update ub( c) as required by Step 2, it suffices to compute the value of the cost at the center of [p], which is of course an upper bound of the optimal value of the cost, and to compare it with the best upper bound available so far. Two methods for eliminating the box [p] under consideration are tried at Step 3. The midpoint test is based on the evaluation of [c] ([p]). If the lower bound of [c]([p] is greater than ub( c), then no p in [p] can be a global minimizer of c over P, so the entire box [p] can be eliminated. The monotonicity test takes advantage of the fact that if c is once differentiable then any unconstrained minimizer corresponds to a stationary point at which the gradient g of c takes the value zero. Thus if an inclusion function [g] is available for g and if

  Consider first guaranteed global optimization. If an upper bound c for the global minimum c of c (p) over [p] 0 is available, then any global minimizer p ∈ [p] has to satisfy c ( p) c.

  2. A division by [∂c/∂p i ] ([p]) is performed in[START_REF] Kühn | Rigorously computed orbits of dynamical systems without the wrapping effect[END_REF]. If standard interval arithmetic is used, a division by an interval containing zero leads to the whole real line, so if 0 ∈ [∂c/∂p i ] ([p]), no contraction takes place along the component of [p] being considered. This situation occurs as soon as [p] contains an argument of a minimum of the cost function, but also when [p] is large, which may lead to a large overestimation of the range of ∂c/∂p i over [p]. Applying extended interval arithmetic [17] when 0 ∈ [∂c/∂p i ] ([p]) would allow one to get up to two intervals of the form ]-∞, a] and [b, +∞[ when performing the division by [∂c/∂p i ] ([p]). After intersection with [p i ], the ith component of C c ([p]) would then be empty or consist of one or two subintervals. Consider now bounded-error parameter estimation. All values of the parameter vector p ∈ S satisfy y m (p) ∈ [y] = y, y , which leads to y m (p)y 0 -y m (p) + y 0 .

y k 0 ,

 0 Copyright c 2009 John Wiley & Sons, Ltd.Int. J. Adapt. Control Signal Process. 2009; 00:1-6 Prepared using acsauth.cls for k = 1, . . . , dim y m (p). The kth output measurement thus leads to a contracted domain[p]

Figure 1 .

 1 Figure 1. Two-compartment model

Figure 2 .

 2 Figure 2. Projection onto the (p1, p2) plane (left) and (p2, p3) plane (right) of a guaranteed outer approximation of the set of all global minimizers (ε = 0.001)

Figure 3 .Figure 3

 33 Figure 3. Volume of the solution set as a function of computing time (in seconds)

Figures 4 and 5

 5 Figures4 and 5show planar projections of the guaranteed outer approximations S ε to the solution set S obtained by the three versions of the Sivia algorithm for four values of the precision parameter. Presentation is in the same order as in TablesII and III. With NIF, the

Figure 4 .

 4 Figure 4. Projection on the (p1, p2)-plane of outer-approximations of the solution set obtained for various values of the precision parameter ε (from left to right, ε = 0.01, ε = 0.005, ε = 0.001, and ε = 0.0005), and for three versions of Sivia (from top to bottom, natural inclusion function (NIF), centered form (CF), and centered form with constraint propagation (CF-CP)); note the scale of the subfigures of the two right-most columns.

Figure 5 .

 5 Figure 5. Projection on the (p2, p3)-plane of outer-approximations of the solution set obtained for various values of the precision parameter ε (from left to right, ε = 0.01, ε = 0.005, ε = 0.001, and ε = 0.0005), and for three versions of Sivia (from top to bottom, natural inclusion function (NIF), centered form (CF), and centered form with constraint propagation (CF-CP)); note the scale of the subfigures of the two right-most columns.

  1. First, one may use some contractor C on [p] to make [p] smaller by eliminating parts of it that do not satisfy f (p) ∈ Y, see Section 3.4.2. If [p] has become empty after contraction, or if [f ] ([p]) has an empty intersection withY, then [p] has an empty intersection with S and should be discarded.3. If [f ] ([p]) is included in Y, then [p] is included in S, and should be stored in S and in S.4. If [f ] ([p]) has a nonempty intersection with Y, but is not entirely in Y, then [p] may contain a part of S without being included in S and is said to be undetermined. If it has a width greater than a prespecified precision parameter ε, then this undetermined box should be bisected into the boxes L [p] and R [p], each of which should be investigated by the algorithm.

  , and provided that a closed-form expression for y m (p) is available, one may use interval constraint propagation to get a contracted domain for[p]. When such expression is not available, one may again resort to a centered inclusion function, now for the model output. For the kth component y m

		translates into
	 	y m k
	

k (p) of y m (p) , for all p ∈ S ⊂ [p] and m ∈ [p] , (19)

Table II .

 II [START_REF] Chernousko | Optimal guaranteed estimates of indeterminacies with the aid of ellipsoids[END_REF] . Computing times in seconds for three versions of the Sivia algorithm

	ε	0.01 0.005 0.001 0.0005
	NIF	21	48	294	823
	CF	650 1265 7720 28328
	CF-CP 180	336	5818 27876
	ε	0.01 0.005 0.001 0.0005
	NIF	1.03	0.31	0.032	0.014
	CF	0.23 0.057 0.093 0.0051
	CF-CP 0.049 0.016 0.0050 0.0040

Table III .

 III Volumes (×10 3 ) of guaranteed outer approximations of the solution set for three versions of the Sivia algorithm

Copyright c 2009 John Wiley & Sons, Ltd. Int. J. Adapt. Control Signal Process. 2009; 00:1-6 Prepared using acsauth.cls

† This work has been partly supported by the ANR CPP and Digiteo-DIM PASO projects


Kyung In Kim

Etienne Roquain

Mark A Van De Wiel

Spatial clustering of array CGH features in combination with hierarchical multiple testing

We propose a new approach for clustering DNA features using array CGH data from multiple tumor samples. We distinguish data-collapsing: joining contiguous DNA clones or probes with extremely similar data into regions, from clustering: joining contiguous, correlated regions based on a maximum likelihood principle. The model-based clustering algorithm accounts for the apparent spatial patterns in the data. We evaluate the randomness of the clustering result by a cluster stability score in combination with crossvalidation. Moreover, we argue that the clustering really captures spatial genomic dependency by showing that coincidental clustering of independent regions is very unlikely. Using the region and cluster information, we combine testing of these for association with a clinical variable in an hierarchical multiple testing approach. This allows for interpreting the significance of both regions and clusters while controlling the Family-Wise Error Rate simultaneously. We prove that in the context of permutation tests and permutationinvariant clusters it is allowed to perform clustering and testing on the same data set. Our procedures are illustrated on two cancer data sets.

Introduction

Array Comparative Genomic Hybridization (array CGH) was designed as a high-resolution measurement device for copy number aberrations, which are known to be involved in the cancer development process [START_REF] Kallioniemi | CGH microarrays and cancer[END_REF]. Through hybridization, fluorescent labeled DNA is extracted from test and reference samples resulting in log-ratios of intensities of the two types of samples. When plotted against chromosomal position, the log-ratio data appears as segments at various amplitudes. One is particularly interested in those chromosomal segments which show levels of loss or gain, because those segments may possibly harbor oncogenic genes, cancer progression markers or other relevant features. Array CGH data is often encoded as loss (deletion), normal and gain by assigning discrete copy number states to chromosomal segments. Generally, two processing steps are performed for this purpose: segmentation and calling.

Segmentation denotes a process to find breakpoints from single copy number profiles. Many segmentation algorithms have been proposed. [START_REF] Willenbrock | A comparison study: applying segmentation to array CGH data for downstream analyses[END_REF]; [START_REF] Lai | Comparative analysis of algorithms for identifying amplifications and deletions in array CGH data[END_REF] discuss and compare several of those, while Dìaz- [START_REF] Dìaz-Uriarte | AdaCGH: A parallelized web-based application and R package for the analysis of aCGH data[END_REF] provide software to simultaneously apply several popular algorithms to the same data set. Recently, new algorithms have been introduced to deal with the computational burden caused by the increasing resolution of the arrays (f.e. [START_REF] Pique-Regi | Sparse representation and Bayesian detection of genome copy number alterations from microarray data[END_REF]. After segmentation, assigning discrete states to the segmented data is referred to as calling [START_REF] Van De Wiel | CGHcall: an algorithm to call aberrations for multiple array CGH tumor profiles[END_REF].

Here, we will assume that the data has been properly segmented and discretized. While the measurement unit is clearly defined from the array platform used, there is no a priori biologically most relevant unit of interest in DNA copy number experiments. This contrasts the majority of mRNA experiments, for which the genes as coding regions for proteins are a natural relevant unit. By nature, the crucial DNA copy number events, aberrations, arise when a piece of DNA is either deleted or gained. Such a piece can be an entire chromosomal arm, but also just 1/10 6 th of an arm. Therefore, we aim to learn the relevant units from the data. Such units are particularly important for interpretation and power when testing association with clinical data. For simplifying terminology, we refer to the array CGH features as probes although this could also reflect clones or cDNAs. Then, our proposed two-step dimension reduction approach enables one to consider the data at two levels with decreasing resolution: regions and clusters. The regions are the result of collapsing neighboring probes of which the discretized values are highly repetitious and redundant within (almost) all profiles (Van de [START_REF] Van De Wiel | CGHregions: dimension reduction for array CGH data with minimal information loss[END_REF], while the clusters contain contiguous similarly behaving regions.

After collapsing, our study proceeds in two directions. First, we develop a model-based clustering algorithm that considers spatial patterns. Second, we perform simultaneous regionwise and cluster-wise multiple testing. The idea of multiple testing based on clustered results is not new. [START_REF] Benjamini | False discovery rates for spatial signals[END_REF] adapt the False Discovery Rate (FDR) to allow for cluster-level multiple testing in the analysis of functional Magnetic Resonance Imaging (fMRI) data. Moreover, comparing to voxel-wise (feature level) multiple testing, they asserted that their method improved SNR (signal to noise ratio) and gained statistical power. Even though both array CGH data and fMRI data have strong spatial correlation patterns, we cannot apply the methodology of [START_REF] Benjamini | False discovery rates for spatial signals[END_REF] directly. Using fMRI it is often feasible to obtain the clusters of voxels from independent, preparatory scans. In array CGH data, however, it is hardly possible to obtain such information a priori. One solution is splitting the set of samples in two parts, where one is used for clustering and the other for testing, which guarantees independence between clustering and testing. However, in that case, we should accept significant power loss due to splitting. Perone [START_REF] Perone Pacifico | False discovery control for random fields[END_REF] offer an alternative solution, also in the context of imaging data and FDR: clustering on the basis of p-values. Their theory is based on smooth Gaussian random fields, which is not a realistic model for array CGH data. Moreover, this would render meaningful clusters from a testing perspective only.

We show that, in the context of permutation testing and family-wise error rate (FWER), it is possible to apply both clustering and testing to the same data set if the cluster result is permutation invariant. We consider selection of a suitable test statistic on both levels (regions and clusters) and use an hierarchical testing procedure to control the FWER. The entire clustering plus testing procedure is illustrated on two array CGH cancer data sets.

High-dimensionality

High-dimensionality is a common theme for the analysis of data produced by high-throughput genomic technologies. It denotes the situation that the number of features is much larger than the number of samples. Due to the large dimension, traditional multiple testing procedures such as Bonferroni yield very conservative results.

While the term 'high-dimensionality' is used neutrally in terms of the number of features, some authors note that technological developments in data generating processes can allow a more specific approach for the data. [START_REF] Benjamini | False discovery rates for spatial signals[END_REF] suggest that the assumption of infinitely many elements in a cluster is not unrealistic for fMRI data in an environment of improved technical resolution. A similar development has occurred for our type of genomic data. Here, the resolution increased from chromosome arms (∼ 40) to Bacterial Artificial Clones (∼ 3000), to oligos (< 400, 000), to next generation arrays (> 1, 000, 000) [START_REF] Park | Experimental design and data analysis for array comparative genomic hybridization[END_REF]. The increasing number of features does not necessarily imply an increasing number of 'distinct' relevant units, because these new features may appear as repetitions of similar copies on a coarser scale.

We briefly illustrate the above idea. Consider the extreme case in which a certain region of the genome spanning r probes never contains a genomic breakpoint for the type of samples under study. Doubling the resolution in this genomic region to 2r will only lead to measuring more of the same (due the absence of a breakpoint) and hence for both resolutions the relevant unit can be collapsed to one data point per sample for this genomic region.

Discretized array CGH data usually appears as large numbers of distinct blocks which consist of probe vectors with (almost) the same discretized status (across samples) for as many as hundreds or thousands of probes. So, regions in the data correspond to distinct blocks. Those regions may be considered as biologically relevant units for the type of samples under study.

Collapsing repetitious probes is useful, since it reduces unnecessary dimension of data and enhances computational convenience. Conventional clustering methods using correlational association among probes seem not appropriate in this stage, because these may yield mingled clusters: repetitiously copied probes mixed with similarly behaving probes. Hence this may weaken the performance and interpretation of the clustering methods usually performed in the next stage. From this point of view, it is desirable to separate two procedures: collapsing and clustering. The first step handles the technical resolution and is used for dimension reduction. The latter step is applied for combining similar behaving probe regions so that one may obtain hints of collective behaviors within certain genomic neighborhoods.

We do not suggest a new algorithm for the collapsing process, but instead use the one proposed by [START_REF] Van De Wiel | CGHregions: dimension reduction for array CGH data with minimal information loss[END_REF]. Here, the relative amount of information lost by collapsing can be controlled; we use 0.5% as an upper bound. Note that this algorithm maintains the high resolution where desired: small genomic segments that differ consequently from their neighbors are kept. In construction of the clustering models, we incorporate the collapsing information into our model in two ways: either via base-pair distance between two regions or via the number of collapsed probes between two regions.

Methods

Mathematical setting

In this paper, the overall statistical model corresponds to n i.i.d. copies of a pair of random variables denoted by (X 1 , Y 1 ), ..., (X n , Y n ), where, for each individual i = 1, ..., n, X i ∈ {0, -1, 1} m is the state vector along all the (ordered) regions and Y i ∈ {0, 1} is a label determining which group the individual belongs to. In this model X = (X 1 , ..., X n ) is a sample of i.i.d. variables, in which the distribution of X i can be seen as a mixture between the two groups Y i = 0 and Y i = 1. In particular, note that the distribution of the X i 's does not involve the relationship between the Y i 's and the X i 's.

In a nutshell, our approach firstly clusters the rows of X that are similar from a modelbased point of view and secondly detects the regions and clusters that are significantly associated to the label vector Y = (Y 1 , ..., Y n ):

• The clustering phase uses the variables X, by assuming a specific parametric model for the distribution of the X i 's . Since X 1 , ..., X n are i.i.d., the clustering result A(X) is invariant under permutations of the columns of X.

• The testing phase tests the independence between regions X 1:n j = (X i j ) 1≤i≤n and Y and between clusters (X 1:n j ) j∈A and Y , conditional on the clustering A(X). The (conditional) p-values are computed by performing permutations of the observed labels, which is valid because A(X) is permutation invariant. Then, these p-values are integrated in an hierarchical multiple testing procedure controlling the family-wise error rate both on the region and cluster level.

Clustering model

Our motivation for spatial clustering comes from Figure 1. In Figure 1, we illustrate Kendall correlations for pairs of all regions for two discretized array CGH data sets. Strong correlations are concentrated in the diagonal parts, so spatial dependency patterns are notably blockwise structures. Such correlations are caused by a high likelihood of the same aberration to span multiple consecutive regions (as opposed to off-diagonal correlations between regions on different chromosomes, which necessarily represent different aberrations). We incorporate the spatial correlation in our modeling approach. Log-linear models are a natural candidate for describing discrete data. Hence, we apply such models for our purpose, while adapting these to take specific data characteristics such as spatial correlations, but also physical base pair distance between regions, into account. Moreover, we consider practical concerns for efficient optimization of the clustering.

Let m be the number of regions and n be the number of i.i.d. samples under consideration. Moreover, let each region be indexed by its spatial location on the DNA, that is, the (j -1)-th and (j + 1)-th regions are the neighbors of the j-th region. The goal is to identify clusters so that regions in one cluster behave similarly.

We first consider the quadratic exponential model [START_REF] Cox | A note on the quadratic exponential binary distribution[END_REF] which approximates log-linear models for binary data up to 2nd order interactions. In our case we consider trinary valued variables. Because the samples are initially assumed i.i.d., we suppress the sample index when describing the model, so X j = x j refers to the observed state of the jth region. Hence, for observations x A = (x j ) j∈A , with x j ∈ {0, -1, 1} and A ⊂ {1, . . . , m}, the marginal distribution of the corresponding random variable X A , p A , is modeled as

log p A (x A ) = α A + j∈A β A,j x j + j<k, j,k∈A γ A,jk f (x j , x k ), where f (x j , x k ) = -1, if x j = x k and f (x j , x k ) = x j x k , otherwise.
Note that α A is not a free parameter, because it is determined by the 'summing-to-one' condition. The discrete values '0', '-1' and '1' denote normal (= 2 copies), loss (< 2 copies) and gain (> 2 copies), respectively. We do not consider amplification (high copy number gain) and double deletion (0 copies) as separate states because of computational complexity of the clustering algorithm, but alternative solutions to accommodate such different states are discussed in the Discussion.

Since our purpose is not to specify every 2nd order interaction γ A,jk , but to bind features which are in the same cluster, it may be better to simplify the above model by assigning one average 2nd order interaction per cluster. Since in actual array CGH data we observe physical base-pair distance between two regions, we use this information in the model. We assume the dependence of two regions becomes weaker for more distant pairs of regions (similar to p.430 of [START_REF] Cressie | Statistics for spatial data[END_REF]). Hence the interaction parameter γ A,jk is represented as γ A djk and djk = d -q jk / max{d -q jk } if j and k are in the same chromosome, 0 otherwise.

where d jk is the base-pair distance between two regions j and k and the maximum is computed over pairs in the same chromosome. Alternatively, the distance function can be used to incorporate the collapsing information; for example, the distance between two regions j and k can be the actual number of probes (measurement unit) between them. Here, q is the decreasing rate of spatial dependence. Common choices for q are 0, 1, 2 [START_REF] Cressie | Statistics for spatial data[END_REF]. Let

θ ℓ = (α A ℓ , β A ℓ ,1 , . . . , β A ℓ ,m ℓ , γ A ℓ )
, where A ℓ is a set of contiguous regions and m ℓ is the number of regions in A ℓ . Then, our full model for a partition A = {A ℓ } K A ℓ=1 of {1, . . . , m}, where K A denotes the number of clusters, is

p A (x) = K A ℓ=1 p A ℓ (x A ℓ ; θ ℓ ), (1)
with log p A ℓ (x A ℓ ; θ ℓ ) = α A ℓ + j∈A ℓ β A ℓ ,j x j + γ A ℓ j<k,j,k∈A ℓ djk f (x j , x k ). (2)
From (1), we know that X A ℓ and X A ℓ ′ are assumed independent for different ℓ and ℓ ′ . Denote the value of x A ℓ for sample i = 1, . . . , n by x i A ℓ . Then, we search a partition and parameter setting that optimizes the following log-likelihood:

( Â, θ) = arg max A∈C,θ n i=1 A ℓ ∈A log p A ℓ (x i A ℓ ; θ ℓ ) (3)
where C is the set of contiguous partitions and θ = (θ ℓ ) K A ℓ=1 . Note that the model ( 1) contains an intrinsic trade-off for dividing a cluster into two subclusters: potentially higher likelihood due to an extra parameter γ (α is not a free parameter) with respect to one cluster versus potentially lower likelihood due to the forced independency for two regions in a different subcluster. Which effect prevails depends on how similar the regions between two subclusters are. As a consequence, the number of clusters K A is not defined a priori, but instead computed through the optimization step.

Generally, searching the space C for large m is computationally intractable. Therefore, we restrict the maximum size of a cluster. The components of the product in the right hand side of (3) are computed in two steps; first, we estimate the parameters of each cluster smaller than the maximum size. From these parameters, the likelihood of each candidate contiguous partition can be determined. Second, given the likelihoods, we search the partition space and, due to the specific partition structure [START_REF] Cormen | Introduction to Algorithms[END_REF], we achieve the global optimum for the equation ( 3) by applying Dijkstra's shortest path algorithm. Note that, from a biological perspective, a cluster should not cover more than one chromosome, so we apply the algorithm per chromosome separately.

One could question the usefulness of this cluster algorithm when the results would be sensitive to similar region-wise differences between, say, two clinical groups. However, observe from (2) that the likelihood for a given partition A ℓ is to a large extent determined by withinsample similarity. This means that, as opposed to a p-value type clustering method, our method is unlikely to cluster two independent consecutive regions that coincidentally have a similar group difference, because the independence implies that these two differences are likely to be caused by (partly) different samples. Moreover, in terms of causality, the perception is that the DNA copy number data (which are early onset markers) potentially affects the group labels, rather than the other way around. Nevertheless, in rare cases, two independent consecutive regions may have realizations that coincide for many samples, in particular when the marginal distributions within groups would (almost) degenerate (e.g. all 0's in one group and all 1's in the other). This may cause coincidental clustering. Therefore, we investigate how likely such coincidental clustering is for a given data set using a shuffling argument.

Denote the event that a random pair of consecutive regions clusters by E, and the event that these two regions are dependent due to their genomic proximity by D. Its complement, lack of genomic dependency, is denoted by D ′ . If our cluster method performs properly P (D|E) should be very high. Applying Bayes' rule, we have

P (D|E) = P (E|D)P (D) P (E) = P (E|D)P (D) P (E|D)P (D) + P (E|D ′ )P (D ′ ) . ( 4
)
From the data we easily estimate P (E) by counting consecutive pairs of regions that cluster. Then, if P (E|D ′ ) is very small and P (E) is fairly large, we necessarily have: P (E|D)P (D) ≫ P (E|D ′ )P (D ′ ), which implies P (D|E) ≈ 1. To show that P (E|D ′ ) is small, we break the genomic dependency structure by shuffling regions such that two regions of the same chromosome are not allowed to be neighbors in the shuffled data set. Then, we apply our clustering method to the shuffled data set. If very few regions cluster, then, necessarily, P (E|D ′ ) is small. The shuffling is repeated several times to account for the arbitrariness of the shuffling.

Clustering stability

Since a clustering result depends on the data used, especially on the sample size, there exists uncertainty on the clustering result. We suggest to use the adjusted Rand index [START_REF] Hubert | Comparing partitions[END_REF] to evaluate stability of the clustering results. If a small change of the data results in large change of the clustering result, then the clustering method is instable. Therefore, we check the stability using a cross validation framework. Suppose regions are numbered from 1 to m and a clustering result is represented as a contiguous partition of {1, . . . , m} as in the previous section.

Let A = {A ℓ } K A ℓ=1 and A ′ = {A ′ ℓ ′ } K A ′
ℓ ′ =1 be two such partitions which we aim to compare and let |A| be the cardinality of index subset A. Then the adjusted Rand index (ARI) for the two partitions is defined as

ARI(A, A ′ ) = r -E(R) max(R) -E(R) , r = K A ℓ=1 K A ′ ℓ ′ =1 |A ℓ ∩ A ′ ℓ ′ | 2 ,
where r is the observed value of random variable R. The expectation of R is computed conditional on the fixed margins,

|A ℓ |, |A ℓ ′ |, ℓ = 1, . . . , K A , ℓ ′ = 1, . . . , K A ′
and the maximum is computed as the maximum of R regardless of the fixed margins. The adjusted Rand index takes numbers between -1 to 1. After obtaining V clustering results from V -fold cross validation, we compute the average adjusted Rand index, aARI, which serves as our stability measure:

aARI = 1 V V v=1 ARI(A, A (v) ) ( 5
)
where A is the original clustering result from the full data and A (v) is the v-th cross-validated clustering result.

Hierarchical multiple testing for clusters and regions

Once we have the clusters A = {A ℓ } K A ℓ=1 , our focus turns to testing association of the region and cluster-wise data with clinical information. Let us first shortly consider association tests between groups (reference and test) and states (loss, normal and gain) for a 2 × 3 contingency

For a cluster A ℓ , we simultaneously test

H ℓ 0 : j∈A ℓ H 0j . (7)
Besides the usual scala of test statistics for a particular testing problem, two classes of approaches may be distinguished: unconditional and conditional tests. In the first case, usually only the sample sizes (row margins) are fixed while column margins are variable. In the second case, both margins are fixed. Permutation tests fall in the latter category, since these only permute the labels (reference or test) for all samples. Permutation tests are useful to approximate (summaries of) multivariate distributions, which is exactly what we need for testing association on the level of clusters. Therefore, we use permutation tests in combination with the popular Pearson X 2 statistic for determining the significance of regions and clusters.

A test statistic for testing ( 7) is M ℓ = min j∈A ℓ p j , where p j is the p-value obtained for region j in cluster A ℓ . Using the p-values to define the cluster test statistic standardizes the regions within a cluster. Before further motivating the use of M ℓ , we first introduce the hierarchical testing approach.

The hierarchical multiple testing procedure we propose to use for the cluster-region data is based on that by [START_REF] Meinshausen | Hierarchical testing of variable importance[END_REF]. Such a procedure allows to test both clusters and regions within one multiple testing framework. The procedure controls Family-Wise Error Rate (FWER) for hierarchical hypotheses. Why FWER and not the more popular False Discovery Rate (FDR)? Firstly, on the level of regions (usually a few hundreds) FDR may not reflect what it should, because a large 'cluster' of highly correlated regions could have a disproportional large share in the number of discoveries, which highly impacts the estimated FDR for other regions. This is a consequence of the fact that FDR does not provide control on subsets [START_REF] Finner | On the false discovery rate and expected type I errors[END_REF]. Secondly, our constructed clusters are more independent than regions, but their number is much lower (often below 100). Then, the power enhancement of an FDR procedure w.r.t. FWER is usually quite subtle, which may not outbalance the stronger conclusion one is allowed to draw with FWER control.

The procedure by [START_REF] Meinshausen | Hierarchical testing of variable importance[END_REF] assumes an a priori known clustering. Such a setting can be relevant for these data as well, but often it is preferable to use the same data for both clustering and testing (see Discussion). In a clustering plus permutation-testing framework, control of the FWER conditional to the data-based clustering is feasible when the clustering is permutation invariant: when testing H ℓ 0 using permutations, we permute the clinical responses (group labels Y in our case), while keeping the regions X 1:n j fixed for j ∈ A ℓ . If the result of the cluster algorithm is column-wise permutation invariant, the clusters may, in the testing phase, be assumed to be known when the null-hypotheses are formulated conditionally to the clusters. More mathematical details on the validity of permutation tests in this setting are presented in Appendix II. A similar argument for combining permutation testing with clustering is given by [START_REF] Goeman | The inheritance procedure: multiple testing of tree-structured hypotheses[END_REF].

The cluster method introduced in this paper is permutation-invariant, because it does not use the group labels. We use the sequential hierarchical testing procedure proposed by [START_REF] Goeman | The inheritance procedure: multiple testing of tree-structured hypotheses[END_REF], with critical values that depend on the rejection set R, which contains both the rejected clusters and regions. This procedure is a slightly more powerful alternative to the one by [START_REF] Meinshausen | Hierarchical testing of variable importance[END_REF]. It applies the so-called inheritance principle in combination with the [START_REF] Shaffer | Modified sequentially rejective multiple test procedures[END_REF] improvement in an hierarchical testing context. The inheritance principle is a variation of the fall-back principle [START_REF] Wiens | The fallback procedure for evaluating a single family of hypotheses[END_REF] that allows to test an hypothesis less strictly (applying a larger p-value threshold) when neighboring hypotheses are rejected. Likewise, the Shaffer improvement allows use of a larger p-value threshold for regions in a cluster, because it uses the connection between the cluster and region hypotheses: if H ℓ 0 is false, at least one H 0j should be false too for j ∈ A ℓ (see Appendix II). We emphasize that this procedure guarantees strong control of FWER [START_REF] Goeman | The inheritance procedure: multiple testing of tree-structured hypotheses[END_REF]. It relies only on the individual permutation p-values of clusters and regions. It is a Holm-type procedure and hence subset pivotality [START_REF] Dudoit | Multiple hypothesis testing in microarray experiments[END_REF] is not required.

For a given rejection set R (both for clusters in A and regions in {1, ..., m}), critical values for clusters are defined by α R = α/(K A -D R ), where D R equals the number of clusters in R for which all regions are members of R as well. Note that we opt, as opposed to [START_REF] Meinshausen | Hierarchical testing of variable importance[END_REF], to weigh clusters equally, because for our application small clusters may be as relevant as large ones. This does not affect control of FWER (see Appendix III for a proof). Critical values for regions in cluster A ℓ are denoted by α R,ℓ = α R /(|A ℓ |max(1, D R,ℓ )), where D R,ℓ equals the number of regions in cluster A ℓ that are members of R.

Then, the hierarchical testing procedure, which is initiated by R = ∅, is as follows.

1. Reject H ℓ 0 for cluster

A ℓ if P 0 (M ℓ ≤ m ℓ ) ≤ α R , where m ℓ is the realization of M ℓ . 2. If H ℓ 0 is rejected, reject H 0j for region j in cluster A ℓ if p j ≤ α R,ℓ .
3. Update the rejection set R and critical values α R,ℓ .

4. Repeat steps 2 and 3 for the non-rejected regions until no more regions are rejected.

5. Update cluster critical values α R .

6. Repeat steps 1 to 5 for the non-rejected clusters and regions therein.

7. Stop when no hypothesis is rejected anymore.

One could argue that for testing clusters M ℓ = min j∈A ℓ p j may have less power than a statistic that focuses more on 'average behavior' (such as a median p-value or a sum of standardized region-wise test statistics). This is true when small effects add up to one larger cluster-wise effect, which is quite common in mRNA gene expression studies. However, we believe the following scenarios to be more relevant for the clustered array CGH data: 1. a cluster is homogenous (large γ), and hence high positive dependencies between regions within the cluster are present; 2. a cluster is heterogenous, and for only a few regions the association exists. It is clear that in the first scenario little power is lost when using min p with respect to f.e. median p, while in the second scenario median p has less power than min p. Using the hierarchical testing approach, the latter could, after rejection of H ℓ 0 , still identify significant regions in the heterogenous cluster.

Results

We analyze two data sets from [START_REF] Chin | Genomic and transcriptional aberrations linked to breast cancer pathophysiologies[END_REF] (Data1) and from Douglas et al. (2004) (Data2). Both data sets have been discretized to ternary values, 0, -1 and 1 ( [START_REF] Van De Wiel | CGHcall: an algorithm to call aberrations for multiple array CGH tumor profiles[END_REF]. After collapsing, Data1 contains 383 regions in rows and 96 and 49 samples in columns, which correspond to ER-positive and ER-negative breast cancer samples, respectively. Data2 contains 436 regions in rows and 7 and 30 samples in columns involved in colorectal cancer, representing two group states: microsatellite instable and chromosomal instable, respectively.

We first concentrate on the clustering results. For the clustering algorithm, the maximum number of regions per cluster is constrained to 9, which we generally found to be sufficiently large. Parameter q, the decreasing rate of spatial dependence is set to 1. Finally, γ is re-scaled by γ = (e γ -1)/(e γ + 1) so that it lies between -1 to 1 as the nominal correlation coefficient (equation (6.5.10) in p.430 of [START_REF] Cressie | Statistics for spatial data[END_REF]).

Figures 2 and3 illustrate the clustering results for Data1 and Data2, respectively. The number of clusters is 63 and 69, respectively. The 10-fold cross-validation stability score, the adjusted average Rand index (5), is high for both data sets: 0.969 and 0.963, respectively. These indicate that the cluster results do not strongly depend on the in-or exclusion of 10% of the samples.

For both data sets we also investigate the probability on coincidental clustering. More precisely, we show that for both data sets the probability that two (consecutive) regions in a cluster are really dependent is high. Following the argument outlined in the Methods section (see equation ( 4)) we first need to show that the probability that two consecutive regions cluster, P (E), is fairly large. Given the relatively small number of clusters in both data sets, this is the case. Next, we need to show that probability that two independent regions cluster, P (E|D ′ ), is low by considering shuffled data sets. Indeed, using 25 shuffled data sets, we observe that for Data 1 (383 regions), on average 382.2 (range: 381-383) clusters are formed, while for Data 2 (436 regions) on average 434.5 (range: 432-436) clusters are formed. Hence, we are confident that clusters in the original data sets are almost always created because of genomic spatial dependency.

For testing association of the class labels with the group labels, 20.000 permutations of the group labels are performed. The above hierarchical testing procedure was applied to the clusters and regions. Tables 2 and3 contain the results of two significant clusters for both data sets, whereas Figures 4 and5 combine the overall collapsing, clustering and testing results. Note that for both data sets the two groups are rather unbalanced. Hence, when a peak in the differential proportion of either gains or losses is caused by a large proportion of that aberration in the smaller group, these tend to be less significant (f.e. second cluster on chromosome 16, Figure 4).

To elucidate the potential benefit of the clusters in a testing setting we compared the results of the hierarchical testing procedure with a simple Holm step-down procedure applied to the regions alone. Two FWER cut-offs were considered: α = 0.05 and α = 0.1. Our findings can be summarized as follows. For both data sets and both cut-offs the Holm procedure only identifies regions that are part of a cluster identified as significant by the hierarchical procedure. For Data1, α = 0.05, the hierarchical procedure identifies 8 clusters, two of which contain no regions that are identified by the Holm procedure. Holm identifies 13 regions, hierarchical two extra, 15. For Data1, α = 0.1, the hierarchical procedure identifies 9 clusters, for which each contains at least one region that is identified by the Holm procedure. Both procedures identify the same 18 regions. For Data2, α = 0.05, the hierarchical procedure identifies 1 significant cluster of which none of the regions are identified by Holm. No significant regions are identified by both. For Data2, α = 0.1, the hierarchical procedure identifies 3 significant clusters. Two clusters contain not a single significant region according to Holm. Both Holm and hierarchical identify (the same) one region. In summary, both procedures are comparable in terms of the number of identified regions, but the hierarchical procedure has a clear advantage: it is able to identify clusters of which none of the regions are identified by Holm. [START_REF] Meinshausen | Hierarchical testing of variable importance[END_REF] shows similar results using an hierarchical Bonferroni-type procedure for other data types.

We also compared the hierarchical procedure with a Holm step-down procedure on clusters. As expected, the number of detections using these procedures differs very little, because clusters are the highest level in the hierarchy and the proportion of differential clusters is small.

Discussion

We introduced a conceptual idea for spatial high-dimensional data: dimension reduction in two steps, collapsing and clustering. These two steps are performed separately: fast collapsing on all array features and more rigorous, model-based clustering on the resulting regions, much fewer in number.

Conditional versus unconditional null-hypotheses

The examples presented here require only one data set to be available for the purpose of both clustering and testing. This implies that the null-hypotheses are conditional on the clusters. However, the methodology developed here also applies to the unconditional setting, where the inference is performed on samples independent of those used for clustering. First, note that since the clustering algorithm always separates chromosomes, the conditional and unconditional nulls are likely to be very accordant on the chromosome level. 0.0000 0.0000 0.0000 0.7999 1 0.0034 1.0000 0.2245 16 0.0000 0.0000 0.0000 0.7999 2 0.0009 0.3060 0.1139 16 0.0000 0.0000 0.0000 0.7999 3 0.0000 0.0000 0.0000 16 0.0000 0.0000 0.0000 0.7999 4 0.0000 0.0000 0.0000 16 0.0000 0.0000 0.0000 0.7999 5 0.0001 0.0374 0.0268 16 0.0000 0.0000 0.0000 0.7999 6 0.0001 0.0374 0.0268 16 0.0000 0.0000 0.0000 0.7999 7 0.0001 0.0189 0.0201 16 0.0000 0.0000 0.0000 0.7999 8 0.0001 0.0189 0.0201 In the conditional setting, we emphasize that the FWER on clusters estimated from the same data is well-defined, albeit from conditional, and hence random, null hypotheses (because the clusters are random). This point of view was investigated by Perone [START_REF] Perone Pacifico | False discovery control for random fields[END_REF] in the context of a clustered version of FDR. However, as opposed to approach of Perone Pacifico et al. ( 2004) using p-value-based clustering, we argue that the conditional approach introduced here clearly delineates the information in the data which is used to make clusters from the information used to perform testing on the clusters, because the clusters are permutation invariant, so that we obtain an interpretable procedure on the permutation probability space as well.

We acknowledge that the conditional setting is not fully in accordance with the traditional hypothesis testing setting, which is why one might prefer the unconditional setting in some cases. Reasons why we emphasize the conditional setting are the following. Tumor array CGH is still very expensive, and, more importantly, for many clinical studies it is not practically feasible to obtain good quality tumor material for a large amount of samples. Moreover, cancer is a very heterogeneous disease, and hence the power to detect differences between groups for particular genomic locations is likely to be small for small sample sizes. For example, some aberrations occur in only 5% of the entire population, but it may still be relevant to detect such an aberration when (most of) it 5% would belong to one group. Therefore, one often prefers to use as many samples as possible in the testing phase, and not 'sacrifice' samples for the purpose of clustering. However, the unconditional setting may be particularly attractive when aCGH data (X) is externally available, but the response (Y ) is not. In such a case, the external data can anyhow not be used in the testing phase. Our software easily applies to such a setting as well.

Finally, we studied to what extent the results of the conditional and unconditional approaches differ for Data1. We repeatedly split this data set in two parts. The conditional approach uses the second part for both clustering and testing, whereas the unconditional approach uses the first part for clustering and the second for testing. We registered which regions are significant and which regions are a member of a significant cluster, at α = 0.2. Summarizing, on average 94.3% (99.9%) of the regions rejected (non-rejected) by the conditional approach were confirmed by the unconditional approach, whereas on average 93.2% (99.8%) of the regions that are member of a rejected (non-rejected) cluster as determined by the conditional approach were confirmed by the unconditional approach. Hence, it is comfortable to notice that the stability of the clustering implies a high agreement between the unconditional and conditional approach, whereas the latter uses only half of the data.

Other issues and conclusion

A byproduct of our clustering method is the strength measure of association (γ) within a cluster. This helps in interpreting the results of the hierarchical multiple testing on clusters and regions: if γ is relatively small for a significant cluster, one expects that the significance is driven by a few regions within the cluster, which should be reflected in the adjusted regionwise p-values. For large values of γ one expects rather similar adjusted region-wise p-values. Hence, to some extent γ and the region-wise p-values may help in distinguishing genomic regions that are potentially causally related to the response and those that are just correlated with neighboring, more strongly associated regions. In any case, further biological validation, also at other molecular levels, such as mRNA or protein, is needed to decide which genomic DNA regions are really causally related to the response. Some prefer the use of undiscretized rather than discretized array CGH data for testing association with clinical information. We believe our cluster-testing approach to be useful in this setting as well. The clustering would then group features that share the underlying discrete characteristics of the data, while the corresponding undiscretized data would be used to achieve (supposedly) more power in the permutation testing procedure.

For clusters that contain regions that are amplified (high copy number gain) rather than gained, one could apply our algorithm on a 4-digit code (-1,0,1,2) but this may be very timeconsuming. When few losses (-1) are present in such a cluster (which is not unlikely given the presence of amplifications), one may locally re-define the 4-digit code to a 3-digit one where the losses and normals would be joined in one class. Another alternative solution is to keep the 4-digit code, but instead of using an extra parameter in model (1) use the same parameter as for x i = 1 (gains), but now with x i = 2. This would give more weight to two consecutively amplified regions than for two consecutively gained regions. Double deletions may be dealt with in a similar way.

Our clustering algorithm suggests a development for DNA probe design in low dimensional platforms such as MLPA [START_REF] Schouten | Relative quantification of 40 nucleic acid sequences by multiplex ligation-dependent probe amplification[END_REF]. If strongly associated clusters could be verified in further studies, then one may consider to use only one MLPA probe per cluster. Such low dimensional platforms are usually less costly, more reproducible and more straightforward to analyze. Another application of our genomic clusters is the clustering of samples and prediction of clinical outcome. [START_REF] Van Wieringen | Weighted clustering of called aCGH data[END_REF] show that use of genomic regions rather than individual probes to cluster samples based on array CGH data enhances clustering performance. In prediction, the smaller number of clusters may result in simpler and more robust feature selection.

In summary, we developed a dedicated clustering algorithm which, in combination with permutation testing, allows a multiple resolution perspective on association of array CGH data with clinical information. With the introduction of extremely high-resolution data, f.e. obtained by massive parallel sequencing, the need for such methods will only increase.

Figure 1 :

 1 Figure 1: Two correlation (Kendall's τ ) heatmaps for Chin et al. (2006) (left) and Douglas et al. (2004) (right) data sets. Data sets contain 383 and 436 regions, respectively. Regions are plotted according to chromosomal position. Colors represent correlations from -1 (cyan) to 1 (pink).

Figure 2 :

 2 Figure 2: Clustering results for Data1. Losses are plotted in red, normals in black, gains in green. Clusters are order according to chromosomal position from bottom to top and depicted alternately in yellow and blue. Sample labels are plotted on the bottom axis, "+" indicates an ER-positive sample, "-" an ER-negative one.

Figure 3 :

 3 Figure 3: Clustering results for Data2. See caption Figure 2.

Figure 4 :

 4 Figure 4: Clustering plus testing, Data1. Chromosomes on the bottom axis, separated by solid lines. Alternating light-grey and white bars demarcate clusters. Crosses show the association coefficients γ (scale on left axis). Dark-grey bars: significant clusters, tick marks at the top axis: significant regions (α = 0.1). Bottom plot: absolute difference between proportions of gains in the two groups (green) and between proportions of losses (red) (scale on right axis).

Figure 5 :

 5 Figure 5: Clustering plus testing, Data2. See caption Figure 4.

Table 1 :

 1 table n = (n 11 , n 12 , n 13 , n 21 , n 22 , n 23 ) as Table1. Let π gs be the population proportion for Contingency table representation of copy number aberrations between two groups. group g and state s. Then, for each region j we test H 0j : for any g, s, π gs,j = π g,j π s,j .

	group \ state loss normal gain
	reference	n 11	n 12	n 13 n 1+
	test	n 21	n 22	n 23 n 2+
		n +1	n +2	n +3 n ++

Table 2 :

 2 Comparison of p-values for two clusters in Data1. Columns 2-4 and 7-9 denote the raw (unadjusted), Holm and hierarchical p-values for clusters and regions, respectively. Fifth column contains the association coefficient γ.

	Clust. p clust. p clust. p clust.	γ	Reg. p reg. p reg. p reg.
		raw	Holm	Hier.		raw	Holm Hier.
	6	0.0006 0.0366 0.0402 0.4653 1 0.0058 1.0000 1.0000
	6	0.0006 0.0366 0.0402 0.4653 2 0.0076 1.0000 1.0000
	6	0.0006 0.0366 0.0402 0.4653 3 0.0023 0.8004 0.9246
	6	0.0006 0.0366 0.0402 0.4653 4 0.0007 0.2541 0.3283
	6	0.0006 0.0366 0.0402 0.4653 5 0.0002 0.0738 0.0938
	6	0.0006 0.0366 0.0402 0.4653 6 0.0403 1.0000 1.0000
	6	0.0006 0.0366 0.0402 0.4653 7 0.2804 1.0000 1.0000
	6	0.0006 0.0366 0.0402 0.4653 8 0.1696 1.0000 1.0000
	16

Table 3 :

 3 Comparison of p-values for two clusters in Data2. See caption Table 2. Clust. p clust. p clust. p clust.

				γ	Reg. p reg. p reg. p reg.
		raw	Holm	Hier.	raw	Holm Hier.
	66	0.0013 0.0863 0.0875 0.7076 1 0.0006 0.2371 0.1925
	66	0.0013 0.0863 0.0875 0.7076 2 0.0065 1.0000 0.9030
	66	0.0013 0.0863 0.0875 0.7076 3 0.0017 0.7174 0.3570
	66	0.0013 0.0863 0.0875 0.7076 4 0.0012 0.4888 0.3220
	66	0.0013 0.0863 0.0875 0.7076 5 0.0007 0.3010 0.2450
	66	0.0013 0.0863 0.0875 0.7076 6 0.0112 1.0000 0.9030
	68	0.0007 0.0455 0.0455 0.5499 1 0.0336 1.0000 1.0000
	68	0.0007 0.0455 0.0455 0.5499 2 0.0387 1.0000 1.0000
	68	0.0007 0.0455 0.0455 0.5499 3 0.0132 1.0000 1.0000
	68	0.0007 0.0455 0.0455 0.5499 4 0.0119 1.0000 1.0000
	68	0.0007 0.0455 0.0455 0.5499 5 0.0160 1.0000 1.0000
	68	0.0007 0.0455 0.0455 0.5499 6 0.0052 1.0000 1.0000
	68	0.0007 0.0455 0.0455 0.5499 7 0.0002 0.0654 0.0630

Acknowledgements

We thank Jelle Goeman and the Associate Editor for discussing the testing procedures with us and Wessel van Wieringen for providing code to create heatmaps. The second author was partly supported by the French Agence Nationale de la Recherche (ANR) (references ANR-PARCIMONIE, ANR-09-JCJC-0027-01, ANR-09-JCJC-0101-01).

Appendix I: Software

The methods discussed here, both for clustering and hierarchical testing, are implemented in the R-package dnaCplusT. The package contains example data sets, documentation on functions and parameters and the actual R-code. For the implementation of the procedures, we used R-Bioconductor software packages RBGL, partitions and multtest. The package is available from the last author's web site: http://www.few.vu.nl/∼mavdwiel/dnaCplusT.html.

Appendix II: Validity of permutation test

Here, we prove the validity of the permutation tests in our setting. In this section, let us write X j = X 1:n j for short. Under the null H 0,j : Y and X j are independent conditionally on A(X) = A, the following proves that the distribution of (Y, X j ) is equal to the distribution of (Y σ , X j ), conditionally on A(X) = A, for any deterministic permutation σ of {1, ..., n} (where the superscript "σ" codes for a σ-column-wise permutation of the vector or of the matrix):

Similarly, we may prove that under the null

the distribution of (Y, (X j ) j∈A ℓ ) is equal to the distribution of (Y σ , (X j ) j∈A ℓ ), conditionally on A(X) = A, for any deterministic permutation σ of {1, ..., n}.

Null-hypothesis ( 8) is sufficient for applying permutation and allows for FWER control using the inheritance principle. However, in order the sharpen FWER control using [START_REF] Shaffer | Modified sequentially rejective multiple test procedures[END_REF] improvement, we need to assume that (8) is implied by

because the intersection hypothesis contains the logical relationship between the region-wise null-hypotheses and the cluster hypothesis needed to apply this improvement. This assumption means that the dependency between Y and the vector (X j ) j∈A ℓ is fully described by dependencies between Y and X j 's, j

Note the similarity between this assumption and one often made in step-wise regression (with random covariates), where an interaction term X 1 X 2 (which models dependency between Y and (X 1 , X 2 )) is only considered once at least one of the main terms (X 1 or X 2 ) is present in the model.

Appendix III: Proof of FWER control

Here, we provide a proof that the hierarchical multiple testing R controls the FWER both on the cluster and on the region levels. We consider the case where we include the Shaffer correction (the other case is similar). The proof was provided in [START_REF] Goeman | The inheritance procedure: multiple testing of tree-structured hypotheses[END_REF], but for a different type of hypothesis weighting. Our arguments are very similar.

The proof uses the so-called "sequential rejection principle, presented by [START_REF] Goeman | The sequential rejection principle of familywise error control[END_REF] (see also [START_REF] Arlot | Some nonasymptotic results on resampling in high dimension. II. Multiple tests[END_REF]). Our hierarchical rejection procedure R, rejecting both clusters A ∈ A and regions j ∈ A, can be expressed as the following sequential rejective procedure: R = R k where R 0 = ∅ and for all i ≥ 0,

and where k is the first i ≥ 0 for which R i+1 = R i . In the above recursion relation, p A and p j are the p-values for cluster A and region j, respectively, while the threshold on clusters

In the latter thresholds, we use the convention 1/0 = +∞ and 0/0 = 0.

We aim to establish that our procedure controls the hierarchical FWER at a pre-specified level α. From the sequential rejection principle, the latter is true, as soon as both a "monotonicity condition" and a "single step condition" hold [START_REF] Goeman | The sequential rejection principle of familywise error control[END_REF]. The first condition is satisfied because α R and α R,A are nondecreasing in R. The second condition is satisfied if we have

where we denoted by H 0 the set of the true clusters and regions and by H c 0 its complementary. We now prove (9): we have

. Now, the sum appearing in the right hand side of the above relation can be taken only over

Moreover, we may use the logical relation between the hierarchical hypotheses saying that if a cluster is false then at least one of its regions is false:

Combining these two facts, we obtain

Thus, we have

which proves (9) and the required FWER control.